Development, calibration, and validation of a novel human ventricular myocyte model in health, disease, and drug block

*For correspondence: [email protected] (JT); [email protected] (BR) Competing interests: The authors declare that no competing interests exist. Funding: See page 21 Received: 29 May 2019 Accepted: 18 December 2019 Published: 24 December 2019 Reviewing editor: Jose´ D Faraldo-Go´ mez, National Heart, Lung and Blood Institute, National Institutes of Health, United States Copyright Tomek et al. This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited. Development, calibration, and validation of a novel human ventricular myocyte model in health, disease, and drug block Jakub Tomek 1 *, Alfonso Bueno-Orovio 1 , Elisa Passini 1 , Xin Zhou 1 , Ana Minchole 1 , Oliver Britton 1 , Chiara Bartolucci 2 , Stefano Severi 2 , Alvin Shrier 3 , Laszlo Virag 4 , Andras Varro 4 , Blanca Rodriguez 1 * 1 Department of Computer Science, British Heart Foundation Centre of Research Excellence, University of Oxford, Oxford, United Kingdom; 2 Department of Electrical, Electronic, and Information Engineering "Guglielmo Marconi", University of Bologna, Bologna, Italy; 3 Department of Physiology, McGill University, Montreal, Canada; 4 Department of Pharmacology and Pharmacotherapy, Faculty of Medicine, University of Szeged, Szeged, Hungary Abstract Human-based modelling and simulations are becoming ubiquitous in biomedical science due to their ability to augment experimental and clinical investigations. Cardiac electrophysiology is one of the most advanced areas, with cardiac modelling and simulation being considered for virtual testing of pharmacological therapies and medical devices. Current models present inconsistencies with experimental data, which limit further progress. In this study, we present the design, development, calibration and independent validation of a human-based ventricular model (ToR-ORd) for simulations of electrophysiology and excitation-contraction coupling, from ionic to whole-organ dynamics, including the electrocardiogram. Validation based on substantial multiscale simulations supports the credibility of the ToR-ORd model under healthy and key disease conditions, as well as drug blockade. In addition, the process uncovers new theoretical insights into the biophysical properties of the L-type calcium current, which are critical for sodium and calcium dynamics. These insights enable the reformulation of L-type calcium current, as well as replacement of the hERG current model. Introduction Human-based computer modelling and simulation are a fundamental asset of biomedical research. They augment experimental and clinical research through enabling detailed mechanistic and systematic investigations. Owing to a large body of research across biomedicine, their credibility has expanded beyond academia, with vigorous activity also in regulatory and industrial settings. Thus, human in silico clinical trials are now becoming a central paradigm, for example, in the development of medical therapies (Pappalardo et al., 2018). They exploit mature human-based modelling and simulation technology to perform virtual testing of pharmacological therapies or devices. Human cardiac electrophysiology is one of the most advanced areas in physiological modelling and simulation. Current human models of cardiac electrophysiology include detailed information on the ionic processes underlying the action potential such as the sodium, potassium and calcium ionic currents, exchangers such as the Na/Ca exchanger and pumps such as the Na/K pump. They also include representation of the excitation-contraction coupling system in the sarcoplasmic reticulum, an important modulator of the calcium transient, through the calcium-induced calcium-release mechanisms and the SERCA pump. Several human models have been proposed for ventricular electrophysiology, and amongst them the ORd model (O’Hara et al., 2011). Its key strengths are the Tomek et al. eLife 2019;8:e48890. DOI: https://doi.org/10.7554/eLife.48890 1 of 48 RESEARCH ARTICLE representation of CaMKII signalling, capability to manifest arrhythmia precursors such as alternans and early afterdepolarisation, and good response to simulated drug block and disease remodelling (Dutta et al., 2016;Dutta et al., 2017a;Passini et al., 2016;Tomek et al., 2017). Consequently, ORd was selected by a panel of experts as the model best suited for regulatory purposes (Dutta et al., 2017a). Most of the ORd model development has focused on repolarisation properties such as its response to drug block, repolarisation abnormalities and its rate dependence. However, a more holistic comparison of ORd-based simulations with human ventricular experimental data reveals important inconsistencies. Firstly, the plateau of the action potential (AP) is significantly higher in the ORd model than in experimental data used for ORd model construction (O’Hara et al., 2011; Britton et al., 2017) and in data from additional studies using human cardiomyocytes (Coppini et al., 2013;Jost et al., 2013). Secondly, the dynamics of accommodation of the AP duration (APD) to heart rate acceleration, which are known to be modulated by sodium dynamics, show only limited agreement with a comparable experimental dataset (Franz et al., 1988;O’Hara et al., 2011). Thirdly, we identify that simulations of the sodium current block has an inotropic effect in the ORd model, increasing the amplitude of the calcium transient, in disagreement with its established negatively inotropic effect in experimental/clinical data (encainide, flecainide, and TTX) (Gottlieb et al., 1990;Tucker et al., 1982;Legrand et al., 1983;Bhattacharyya and Vassalle, 1982). All those properties, namely AP plateau potential, APD adaptation and response to sodium current block, have strong dependencies on sodium and calcium dynamics. We therefore hypothesise that ionic balances during repolarisation require further research. We specifically focus on an indepth re-evaluation of the L-type calcium current (I CaL ) formulation, given its fundamental role in determining the AP, the calcium transient and sodium homeostasis through the Na/Ca exchanger. The second main focus is the re-assessment of the rapid delayed rectifier current (I Kr ), the dominant repolarisation current in human ventricle, under conditions that reflect experimental data-driven plateau potentials. Using a development strategy based on strictly separated model calibration and validation, we sought to design, develop, calibrate and validate a novel model of human ventricular electrophysiology and excitation contraction coupling, the ToR-ORd model (for Tomek, Rodriguez – following ORd). Our aim for simulations using the ToR-ORd model is to be able to reproduce all key depolarisation, repolarisation and calcium dynamics properties in healthy ventricular cardiomyocytes, under drug block, and in key diseased conditions such as hyperkalemia (central to acute myocardial ischemia), and hypertrophic cardiomyopathy. eLife digest Decades of intensive experimental and clinical research have revealed much about how the human heart works. Though incomplete, this knowledge has been used to construct computer models that represent the activity of this organ as a whole, and of its individual chambers (the atria and ventricles), tissues and cells. Such models have been used to better understand lifethreatening irregular heartbeats; they are also beginning to be used to guide decisions about the treatment of patients and the development of new drugs by the pharmaceutical industry. Yet existing computer models of the electrical activity of the human heart are sometimes inconsistent with experimental data. This problem led Tomek et al. to try to create a new model that was consistent with established biophysical knowledge and experimental data for a wide range of conditions including disease and drug action. Tomek et al. designed a strategy that explicitly separated the construction and validation of a model that could recreate the electrical activity of the ventricles in a human heart. This model was able to integrate and explain a wide range of properties of both healthy and diseased hearts, including their response to different drugs. The development of the model also uncovered and resolved theoretical inconsistencies that have been present in almost all models of the heart from the last 25 years. Tomek et al. hope that their new human heart model will enable more basic, translational and clinical research into a range of heart diseases and accelerate the development of new therapies. Tomek et al. eLife 2019;8:e48890. DOI: https://doi.org/10.7554/eLife.48890 2 of 48 Research article Cell Biology Computational and Systems Biology Materials and methods Strategy for construction, calibration and validation of the ToR-ORd model Table 1 lists the properties (left column) and key references (right column) of experimental and clinical datasets considered for the calibration (top) and independent validation (bottom) of the ToRORd model. This represents a comprehensive list of properties, known to characterize human ventricular electrophysiology under multiple stimulation rates, and also drug action and disease. The recordings in were obtained in human ventricular preparations primarily using measurements with microelectrode recordings, unipolar electrograms, and monophasic APs, therefore avoiding photon scattering effects or potential dye artefacts present in optical mapping experiments. In addition, the ToR-ORd model was calibrated to manifest depolarisation of resting membrane potential in response to an I K1 block, based on evidence in a range of studies summarised in Dhamoon and Jalife (2005). The calibration criteria are chosen to be fundamental properties of ionic currents, action potential and single-cell pro-arrhythmic phenomena (described in more detail in Appendix 11). The validation criteria include response to rate changes, drug action and disease, to explore the predictive power of the model under clinically-relevant conditions. We initially performed the evaluation of the ORd model (O’Hara et al., 2011) by conducting simulations for each of the calibration criteria in Table 1. Further details are described throughout the Materials and methods section and Appendix 1-15.1. Simulations with the existing versions of the ORd model failed to fulfil key criteria such as AP morphology, calcium transient duration, several properties of the L-type calcium current, negative inotropic effect of sodium blockers, or the depolarising effect of I K1 block. The results are later demonstrated in Figures 2 and 3, and Methods:Calibration of I K1 block and resting membrane potential. Secondly, we attempted parameter optimisation using a multiobjective genetic algorithm (Torres et al., 2012). However, simulations with the ORd-based models were unable to fulfil key criteria such as AP and Ca morphology, and the effect of sodium and calcium block on calcium transient amplitude and APD, respectively. We then proceeded to reevaluate the ionic current formulations based on experimental data and biophysical knowledge. Key currents included I CaL and specifically its driving force and activation, as Table 1. Criteria and human-based studies used in ToR-ORd calibration and validation. Calibration Action potential morphology (Britton et al., 2017;Coppini et al., 2013;Jost et al., 2013) Calcium transient time to peak, duration, and amplitude (Coppini et al., 2013) I-V relationship and steady-state inactivation of L-type calcium current (Magyar et al., 2000) Sodium blockade is negatively inotropic (Gottlieb et al., 1990;Tucker et al., 1982;Legrand et al., 1983;Bhattacharyya and Vassalle, 1982). L-type calcium current blockade shortens the action potential (O’Hara et al., 2011) Early depolarisation formation under hERG block (Guo et al., 2011) Alternans formation at rapid pacing (Koller et al., 2005) Conduction velocity of ca. 65 m/s (Taggart et al., 2000) Validation Action potential accommodation (Franz et al., 1988) S1-S2 restitution (O’Hara et al., 2011) Drug blocks and action potential duration (Dutta et al., 2017a;O’Hara et al., 2011) Hyperkalemia promotes postrepolarisation refractoriness (Coronel et al., 2012) Hypertrophic cardiomyopathy phenotype (Coppini et al., 2013) Drug safety prediction using populations of models (Passini et al., 2017) Physiological QRS and QT intervals in ECG (Engblom et al., 2005;van Oosterom et al., 2000;Bousseljot et al., 1995; Goldberger et al., 2000) Tomek et al. eLife 2019;8:e48890. DOI: https://doi.org/10.7554/eLife.48890 3 of 48 Research article Cell Biology Computational and Systems Biology well as the I Na , I Kr , I K1 and chloride currents. The multiobjective genetic algorithm optimisation was repeated several times, throughout the introduction of structural changes to the model. Once simulations with an optimised model fulfilled all calibration criteria, validation was conducted through evaluation against additional experimental recordings for drug block, disease, tissue and whole-ventricular simulations. Details concerning the simulations are given in Appendix 1-15.1, namely the description of simulation protocols and ionic concentrations used (Appendix 1-15.1.1), representation of heart disease (Appendix 1-15.1.2), 1D fibre simulations (Appendix 1-15.1.3), population-of-models and drug safety assessment (Appendix 1-15.1.4), transmurality and whole-heart simulations with ECG extraction (Appendix 1-15.1.5), and a technical note on the update to the Matlab ODE solver which facilitates efficient simulation of the multiobjective GA (Appendix 1-15.1.6). Unless specified otherwise, the baseline ORd model (O’Hara et al., 2011) was used for comparison with the ToR-ORd model. ToR-ORd model structure The ToR-ORd model follows the general ORd structure (Figure 1A). The cardiomyocyte is subdivided into several compartments: main cytosolic space, junctional subspace, and the sarcoplasmic reticulum (SR, further subdivided into junctional and network SR). Within these compartments are placed ionic currents and fluxes described by Hodgkin-Huxley equations or Markov models. The main ionic current formulations altered compared to ORd are highlighted in orange in Figure 1A. In-depth revision of the L-type calcium current The I CaL current was deeply revisited, particularly with respect to its driving force, based on biophysical principles. This reformulation is of relevance to almost all models of cardiac electrophysiology. The I CaL formulation in the ORd model is based on Hodgkin-Huxley equations, with the total current being a product of three components: 1) Open channel permeability, 2) A set of gating variables determining the fraction of channels being open, 3) The electrochemical driving force which acts on ions to move through the open channel based on the membrane potential and ionic concentrations on both sides of the membrane (more details in Appendix 1-5). In most Hodgkin-Huxley models of cardia currents, the driving force is computed as (V-E ion ), that is, the membrane potential minus equilibrium potential, either computed from the Nernst equation, or measured experimentally. However, starting with the Luo-Rudy model (LRd) of 1994 (Luo and Rudy, 1994), the driving force of ions via I CaL in cardiac models is modelled based on the Goldman-Hodgkin-Katz (GHK) flux equation. The driving force based on the GHK equation is: ’CaL ¼z2VF2 RTS½ iezVF RTS½ o ezVF RT1; where zis the charge of the given ion, Vis the membrane potential, F,R,Tare conventional thermodynamic constants, and [S] i ,[S] o are intracellular and extracellular activities of the given ionic specie. S½ ¼gm, where gis the ionic activity coefficient and mthe concentration (in either the intracellular or extracellular space, yielding S½ ior S½ o). Determining ionic activity coefficients In order to compute the ionic driving force via the GHK equation, it is necessary to know the ionic activity coefficients of the intracellular (g i ) and extracellular (g o ) space. The Luo-Rudy model and other Rudy-family models use g o = 0.341 for extracellular space and g i = 1 for the intracellular space. Models based on the Shannon model (Shannon et al., 2004) use 0.341 for both intracellular and extracellular space, but we were unable to find the motivation for this change. The Debye-Hu ¨ckel theory is commonly used to compute the activity coefficients. We used the Davies equation, which extends the basic Debye-Hu ¨ckel equation to be accurate for ionic concentrations found in living cells (Mortimer, 2008): loggi¼ Az2 iffiffiI p 1þffiffiI p0:3I  ; Tomek et al. eLife 2019;8:e48890. DOI: https://doi.org/10.7554/eLife.48890 4 of 48 Research article Cell Biology Computational and Systems Biology Figure 1. Model structure. (A) A schematic of the novel human ventricular myocyte model for electrophysiology and calcium handling. Orange indicates components, substituted, or added, compared to the original ORd model. ‘SS’ indicates junctional subspace compartment, where calcium influx via L-type calcium current occurs and where calcium is released from the sarcoplasmic reticulum. ‘JSR’ and ‘NSR’ are junctional and network sarcoplasmic reticulum compartments, respectively. ‘Main cytosolic pool’ is the remaining intracellular space. Transmembrane currents are indicated with an ‘I’ in their name, with fluxes indicated as ‘J’. Components with a green underscore are modulated by CaMKII signalling. (B) The structure of the Lu-Vandenberg (Lu et al., 2001) Markov model used for the rapidly activating delayed rectifier repolarisation current (I Kr ). The transition rates are given in Appendix 1-15.3.5. Tomek et al. eLife 2019;8:e48890. DOI: https://doi.org/10.7554/eLife.48890 5 of 48 Research article Cell Biology Computational and Systems Biology where Ais a constant (~0.5 for water at 25˚C, ~0.5238 at 37˚C), z i is the charge of the respective ion, and Iis the ionic strength of the solution. The ionic strength is defined as: I¼0:5i Xmiz2 i; where m i is the concentration of the i-th ionic specie present. For concentrations in a study measuring properties of I CaL (Magyar et al., 2000), Iis ca. 0.15-0.17. This warrants the use of Davies equation, which was shown to be accurate for Iup to 0.5, unlike the basic Debye-Hu ¨ckel equation, which is accurate for Iup to 0.01 only (Mortimer, 2008). We implemented the computation of ionic coefficients based on the Davies equation dynamically, so that the activity coefficients are estimated at every simulation step. This allows accurate representation of the driving force when ionic concentrations are disturbed, such as at varying pacing rates, or during homeostatic imbalance. The dynamic computation is also used to estimate ionic activity coefficients for potassium and sodium flowing through the calcium channels, taking into account their different charge. Throughout our simulations, both intracellular and extracellular activity coefficients generally lie between 0.61 and 0.66. Importantly, this estimate shows that the intracellular and extracellular activity coefficients are relatively similar (corresponding to the broadly similar total concentration of charged molecules), in contrast with the original values. Particularly, the origin of the intracellular activity coefficient g i = 1 in the Luo-Rudy model is unclear, as by the Davies (or by any Debye-Hu ¨ckel variant) equation, Iwould have to be zero, which is possible only when there are no ions present. Activation curve extraction An additional improvement in the I CaL formulation is the estimation of its activation curve. In brief, we implement a consistent use of the GHK equation for the extraction of the activation curve and for the I CaL formulation in the ToR-ORd model. The activation curve is obtained via dividing the experimentally measured I-V relationship of the current by the expected driving force for each pulse potential (see Appendix 1-3 for a graphical overview of the process). However, we identified a theoretical inconsistency in previous cardiac models across species (e.g. Luo and Rudy, 1994; Hund et al., 2008;O’Hara et al., 2011;Shannon et al., 2004;Grandi et al., 2010;Carro et al., 2011): whereas the Nernstian driving force of (V-E Ca ) is used to derive the activation curve, the GHK driving force is then used to calculate I CaL . Indeed, experimental studies reporting the activation curve of I CaL generally use the Nernstian driving force of (V-E Ca ) with E Ca being the experimentally measured reversal potential of approximately 60 mV. This is explicitly stated in Linz and Meyer (2000), and also the activation curve by Magyar et al. (2000) used in the ORd model is consistent with dividing the IV relationship with (V-60). In this study, we propose that, for consistency, the same equation needs to be applied both to obtain the activation curve from the I-V curve and to represent the driving force in the current formulation. Thus, in the ToR-ORd model, the activation curve for I CaL was obtained by dividing the I-V curve from Magyar et al. (2000) by the GHK-based driving force, computed using ionic activity coefficients based on the Davies equation (as explained in the previous Section) and intracellular and extracellular ionic concentrations as in Magyar et al. (2000). The following capped Gompertz function (a flexible sigmoid) was found to be the best fit to the resulting steady-state activation curve: d¥¼1:0763 e1:007e0:0829Vfor V31:4978 j1otherwise  ; where Vis the membrane potential. Other I CaL changes 20% of I CaL was placed in the main cytosolic space, consistent with the literature (Scriven et al., 2010). This increases the plateau-supporting capability of I CaL , given that the myoplasmic I CaL is subject to a weaker calcium-dependent inactivation than I CaL in the junctional subspace. Other minor changes are given in Appendix 1-15.3.3. Tomek et al. eLife 2019;8:e48890. DOI: https://doi.org/10.7554/eLife.48890 6 of 48 Research article Cell Biology Computational and Systems Biology I Kr replacement The calibration of the ToR-ORd model’s AP morphology to experimental data resulted in problematic response to calcium blockade during an early phase of the model development when the original I Kr formulation was used (further details in Appendix 1-12). I CaL block is known to shorten APD experimentally (O’Hara et al., 2011) but resulted in a major APD prolongation in simulations instead. This discrepancy could not be resolved through parameter optimisation. A mechanistic analysis revealed that this follows from the lack of ORd I Kr activation, which is however not consistent with relevant experimental data (Lu et al., 2001). We therefore considered alternative I Kr formulations and specifically the Lu-Vandenberg (Lu et al., 2001) Markov model (Figure 1B). The Lu-Vandenberg I Kr model is based on extensive experimental data allowing the dissection of activation and recovery from inactivation and provided the best agreement with experimental data, specifically when considering the AP plateau potentials reported experimentally. In Appendix 1-12, we: (1) provide a detailed explanation of origins of AP prolongation following I CaL block in a model which manifests experimental data-like plateau potentials and which contains the ORd I Kr formulation; (2) explain why this phenomenon occurs only in a model with experimental data-like plateau potentials, but not in the original high-plateau ORd model; (3) compare the ORd and Lu-Vandenberg I Kr formulations with experimental data, demonstrating the good agreement with experimental data of the Lu-Vandenberg formulation but not the ORd. Following the inclusion of the Lu-Vandenberg I Kr formulation, all models generated during model calibration exhibited APD shortening in response to I CaL block. Changes in I Na, I (Ca)Cl, I Clb and I K1 The I Na current formulation was replaced by an alternative human-based formulation (Grandi et al., 2010), given established limitations of the original model with regards to conduction velocity and excitability (O’Hara et al., 2011), comment on article from 05 Oct 2012). The Grandi I Na model was updated to account for CaMKII phosphorylation (Appendix 1-15.3.1). Also from the Grandi model, we added the calcium-sensitive chloride current I (Ca)Cl and background chloride current I Clb formulation (Grandi et al., 2010). Neither model was changed compared to the original formulations, but the intracellular concentration of Cl - was slightly increased (Appendix 1-15.1.1). In accordance with recent observations, I (Ca)Cl was placed in the junctional subspace (Magyar et al., 2017). The motivation to add these currents was to facilitate the shaping of post-peak AP morphology (via I (Ca)Cl ), with I Clb playing a dual role stemming from its reversal potential of ca. 50 mV. It slightly reduces plateau potentials during the action potential, but during the diastole, it depolarises the cell slightly, improving the reaction to I K1 block as explained in the next subsection. The I K1 model was replaced with the human-based formulation by Carro et al. (2011), as it was shown to be key for simulations of hyperkalemic conditions. The I K1 replacement was done before hyperkalemia simulation, not violating the classification of hyperkalemia criterion as a validation step. Extracellular potassium concentration in a healthy cell was reduced from 5.4 to 5 mM to fall within the physiological range (Zacchia et al., 2016). Calibration of I K1 block and resting membrane potential When evaluating the baseline ORd model against the selected criteria, we observed that a reduction in I K1 results in hyperpolarisation of the cell (from 88 to 88.16 mV at 1 Hz pacing). However, it is established that I K1 reduction depolarises cells experimentally (Dhamoon and Jalife, 2005). Changes made during ToR-ORd calibration (predominantly the altered balance of currents during diastole and the inclusion of background chloride current) result in ToR-ORd manifesting depolarization in response to I K1 block, consistent with experimental data. Multiobjective genetic algorithm We applied a multiobjective genetic algorithm (MGA, @gamultiobj function in Matlab, Deb, 2001) to automatically re-fit various model parameters. Based on preliminary experimentation, we used a two-dimensional fitness. We used MGA rather than an ordinary genetic algorithm or particle swarm optimisation, given that MGA optimises towards a Pareto front rather than a single optimum, implicitly maintaining population diversity. The Pareto front is the set of all creatures which are not Tomek et al. eLife 2019;8:e48890. DOI: https://doi.org/10.7554/eLife.48890 7 of 48 Research article Cell Biology Computational and Systems Biology dominated by any other creature in the population, that is creatures for which there is no other creature better in all fitness dimensions. Therefore, a subpopulation of diverse solutions is maintained, and the optimiser consequently has less of a tendency to converge to a single local optimum compared to single-number fitness approaches. In addition, the crossover operator of GA is well suited for a task where multiple criteria are optimised, given that creatures in the population may efficiently share partial solutions to various subcriteria. The fitness used in this study is described in greater detail in Appendix 1-1. Evaluation pipeline and code To facilitate the model validation and future work, we also provide an automated ‘single-click’ evaluation pipeline. It runs automatic simulations to extract and visualise single-cell biomarkers including those related to AP morphology, effect of key channel blockers, early afterdepolarisations (EAD), and alternans measurement. The pipeline generates a single HTML report containing all the results; see Appendix 1-15.2 for a visualisation. The code for our model (Matlab and CellML), the validation pipeline, and the experimental data on human AP morphology are available at https://github.com/ jtmff/torord (Tomek, 2019; copy archived at https://github.com/elifesciences-publications/ torord). An informal blog giving further insight into the choices we made, as well as general thoughts on the development of ToR-ORd and computer models in general, is available at https://underlid. blogspot.com/. We designed the Matlab code used to simulate our model so that the simulation core is structured into functions computing currents, making the high-level organisation of code clear, and facilitating inclusion of alternative current formulations. In addition, a CellML file encoding our model is also provided. This makes the model readily runnable in several simulators in addition to Matlab (e.g. Chaste [Pitt-Francis et al., 2009] and OpenCOR [Garny and Hunter, 2015]). Furthermore, the Myokit library (Clerx et al., 2016) enables conversion of the CellML file to other languages (such as C or Python). Results Calibration based on AP, calcium transient, and L-type calcium current properties The AP morphology of the ToR-ORd is within or at the border of the interquartile range of the Szeged-ORd experimental data (Figure 2A). This is a major improvement compared to the original ORd morphology, which overestimates plateau potentials, particularly during early plateau (Figure 2A). The fact that the early plateau potential is around 20–23 mV is clearly apparent from experimental recordings and is further corroborated by additional studies in human tissue samples (Jost et al., 2013, Figure 6) and isolated human cardiomyocytes (Coppini et al., 2013). We note that compared to the Szeged-ORd dataset (Britton et al., 2017), our model manifests a slightly increased peak membrane potential in the single-cell form, similar to single-cell experimental data (Coppini et al., 2013). This is a design choice related to the fact that the Szeged-ORd dataset contains recordings of small tissue samples, which are expected to manifest a reduced peak potential compared to single-cell. When coupled in a fibre, ToR-ORd manifests conduction velocity of 65 cm/ s, which is consistent with clinical data (Taggart et al., 2000). Both time to peak calcium and duration of calcium transient at 90% recovery obtained with the ToR-ORd model are within the standard deviation of experimental data in isolated human myocytes (Coppini et al., 2013), whereas ORd slightly overestimated the calcium transient duration (Figure 2B). The calcium transient amplitude of ToR-ORd also matches the Coppini et al. data after accounting for the different APD (Appendix 1-8). As described in Materials and methods, the ToR-ORd I CaL activation curve was extracted from experimental data, using the Goldman-Hodgkin-Katz formulation of ionic driving force, ensuring theoretical consistency, unlike the ORd I CaL formulation (Figure 2C). This considerably improves the results of simulated protocols to obtain IV relationship (Figure 2D), validating the theory-driven changes (see Appendix 1-4 for the demonstration of how the updated activation curve underlies the improvement). The simulation of the protocol measuring steady-state inactivation also reveals improved agreement of ToR-ORd with experimental data compared to ORd (Figure 2E). The Tomek et al. eLife 2019;8:e48890. DOI: https://doi.org/10.7554/eLife.48890 8 of 48 Research article Cell Biology Computational and Systems Biology difference between measured ORd steady state inactivation and the experimental data (ca. two times stronger inactivation at around 15 mV, which is relevant for EAD formation) is initially surprising, given that the equation of ORd I CaL steady-state inactivation curve provides a good fit to the same experimental data. This difference follows from the formulation of calcium-dependent inactivation of I CaL (see Appendix 1-5 for details). We observed that in cases of elevated I CaL (e.g. in midmyocardial cells), ORd reverses current direction towards positive values, which is an unexpected behaviour given its reversal potential of 60 Figure 2. Action potential, calcium transient, and I CaL in ToR-ORd. Action potential (A) and calcium transient (B) at 1 Hz obtained with the ToR-ORd model following calibration, compared to those obtained with the ORd model and experimental data from O’Hara et al. (2011) and Coppini et al. (2013), respectively. The purple and green zones in (B) stand for mean ±standard deviation. The duration of calcium transient at 90% recovery was extracted from figures in Coppini et al. (2013), adding the time to peak and time from peak to 90% recovery. (C) Activation curves used in the ToRORd and ORd models (blue and red lines, respectively). The points correspond to the IV relationship measured in Magyar et al. (2000) normalised by the Nernstian driving force (d.f.) assuming reversal potential of +60 mV (blue points) and the GHK driving force (red points). (D, E) I-V relationship and steady-state inactivation as measured in ToR-ORd (red line) versus ORd model (blue line) versus experimental data from Magyar et al. (2000) (black points with line). I CaL,tot is the sum of currents corresponding to all ions (Ca 2+ , Na + , and K + ) passing through the L-type calcium current channels. (F) L-type calcium current of a midmyocardial cell, showing current reversal in ORd, but not in ToR-ORd. Only the calcium component of I CaL,tot is shown to demonstrate that the current reversal is not due to other ions. We note that the difference in total amplitude of I CaL in Figure 2F follows predominantly from different action potential shape in ToR-ORd vs ORd, consistent with the I-V relationship. Tomek et al. eLife 2019;8:e48890. DOI: https://doi.org/10.7554/eLife.48890 9 of 48 Research article Cell Biology Computational and Systems Biology Figure 7. Populations of models and drug safety prediction. (A) Percentile-based summary of AP and calcium transient traces for the two populations of human ToR-ORd models (left side) and distribution of ionic current conductances among models in the population (right side), in the ranges [50150]% (top row) and [0–200]% (bottom row) of the baseline values. (B) A comparison of 0–200% populations based on ToR-ORd (left) and ORd (right, based on Passini et al., 2017) in response to high-dose mexiletine (100-fold effective therapeutic dose). Traces classified as EADs are plotted in red (manifesting only in the ORd population). (C) TdP score obtained for simulations of the 62 reference compounds, based on the occurrence of druginduced repolarisation abnormalities at all tested concentrations in the [0–200]% population of ToR-ORd models. The colours associated with drugs signify their established torsadogenicity as specified in Appendix 1-15.1.4. The logarithmic scale was considered to maximise the visual separation between safe and risky drugs. The classification based on the TdP score is summarised as a confusion matrix on the right, and also compared with the corresponding results obtained in a population of models based on the original ORd model (Passini et al., 2017). Tomek et al. eLife 2019;8:e48890. DOI: https://doi.org/10.7554/eLife.48890 16 of 48 Research article Cell Biology Computational and Systems Biology mechanical remodelling, the disease predisposes the hearts to arrhythmia formation, increasing the vulnerability to early afterdepolarisations. HCM induces complex multifactorial remodelling of cell electrophysiology and calcium handling, making it a challenging validation problem for a computer model. We applied the available human experimental data on HCM remodelling (based predominantly on Coppini et al., 2013) to our baseline model using an approach similar to Passini et al. (2016), observing that the dominant features of the remodelling observed by Coppini et al. are captured. The HCM variant of the computer model corresponds to experimental data in the AP morphology, manifesting a significantly higher plateau potential and an overall APD prolongation (Figure 8C). The calcium transient amplitude of the HCM model is slightly reduced, has longer time to peak, and a noticeably longer duration at 90% recovery (Figure 8D), also consistent with the data by Coppini et al. (2013). Ultimately, the HCM variant of our model is more prone to the formation of EADs (Figure 8E), as was shown experimentally (Coppini et al., 2013). This difference is in line with postulated key role of I CaL and NCX in EAD formation (Luo and Rudy, 1994;Weiss et al., 2010), both of which are markedly increased in HCM. Excessive prolongation of APD due to a strong increase in late sodium current in HCM also contributes to the EAD formation as well, as shown by Coppini et al. (2013). Figure 8. Simulation of hypertrophic cardiomyopathy (HCM) and hyperkalemia. (A) The effect of hyperkalemia on AP morphology; measured in the centre of a simulated fibre. (B) APD90 and effective refractory period (ERP) at varying extracellular potassium concentration. For extracellular potassium higher than 9 mM, full AP did not develop, but low-amplitude activation propagated through the fibre (Appendix 1-14). Membrane potential (C) and calcium transient (D) at 1 Hz pacing compared between a single healthy and HCM cell. (E) 50% I Kr block induces EADs in HCM cell, but not in a healthy one. Tomek et al. eLife 2019;8:e48890. DOI: https://doi.org/10.7554/eLife.48890 17 of 48 Research article Cell Biology Computational and Systems Biology Validation: human whole-ventricular simulations - from ionic currents to ECG We conducted 3D electrophysiological simulations using the ToR-ORd model, representing the membrane kinetics of endocardial, epicardial and mid-myocardial cells to investigate their ability to simulate the ECG (see Appendix 1-15.1.5). Transmural and apex-to-base spatial heterogeneities as well as fibre orientations based on the Streeter rule were incorporated into a human ventricular anatomical model derived from cardiac magnetic resonance (Lyon et al., 2018). Figure 9A shows the resulting electrocardiogram computed based on virtual electrodes positioned on a torso model shown in Figure 9B. The ECG manifests a QRS duration of 80 ms (normal range 78 ±8 ms), and a QT interval of 350 ms (healthy:<430 ms); all of these quantitative measurements are in the range of ECGs of healthy persons (Engblom et al., 2005;van Oosterom et al., 2000). ECG morphology also showed normal features, such as R wave progression in the precordial leads from V1 to V6, isoelectric ST segment, and upright T waves in leads V2 to V6, with inverted T wave in aVR. Figure 9C shows the activation sequence is in agreement with Durrer et al. (1970). The APD map shows longer APD in the endocardium and the base, and shorter APDs in the epicardium and the apex, respectively (Figure 9D). Discussion In this study, we present a new model of human ventricular electrophysiology and excitation contraction coupling, which is able to replicate key features of human ventricular depolarisation, repolarisation and calcium transient dynamics. The ToR-ORd model was developed using a defined set of calibration criteria and subsequently validated on features not considered during calibration to demonstrate its predictive power. This article also unravels several important theoretical findings with implications for computational electrophysiology reaching beyond the ToR-ORd model and cardiac electrophysiology: firstly, the reformulation of the L-type calcium current, which is broadly relevant and generally applicable to human and other species, and secondly, the mechanistically guided replacement of I Kr . Discovering the necessity to carry out these theoretical reformulations was enabled by the comprehensive set of calibration criteria and the use of a genetic algorithm to fulfil them. Finally, to enable reproducibility, we openly provide an automated model evaluation pipeline, which provides a rapid assessment of a comprehensive set of calibration or validation criteria. The AP morphology of ToR-ORd is in agreement with the Szeged endocardial myocyte dataset used to construct the state-of-the art ORd model (O’Hara et al., 2011). The agreement is considerably better than that of ORd itself, which has important implications for multiple aspects studied in this work. The calcium transient also recapitulates key features of human myocyte measurements (Coppini et al., 2013). The validation of ToR-ORd shows that the model responds well to drug block with regards to APD (Dutta et al., 2017a). Good APD accommodation (reaction to abrupt, but persisting changes in pacing frequency) indicates a good balance between ionic currents (Franz et al., 1988;Pueyo et al., 2010). Replication of arrhythmia precursors such as early afterdepolarisations (Guo et al., 2011) and alternans (Koller et al., 2005) makes the model useful for simulations and understanding of arrhythmogenesis. This is particularly important in the context of heart disease, where ToR-ORd is shown to replicate key features of hyperkalemia (Coronel et al., 2012) and hypertrophic cardiomyopathy (Coppini et al., 2013). The model is also shown to be promising in drug safety testing, and whole-heart simulations demonstrate physiological conduction velocity (Taggart et al., 2000) and produce a plausible ECG signal. Among the improved behaviours compared to the state-of-the-art ORd model (O’Hara et al., 2011), the good response of the ToR-ORd model to sodium blockade is particularly noteworthy. ToR-ORd predicts the negative inotropic effect of sodium blockade, consistent with data (Gottlieb et al., 1990;Tucker et al., 1982; Legrand et al., 1983;Bhattacharyya and Vassalle, 1982), unlike ORd, which suggests a strong pro-inotropic effect. The improvement in ToR-ORd follows from the relatively complex interplay of the theoretically driven reformulation of the L-type calcium current and data-driven changes to the AP morphology. This result is of great importance in the context of pharmacological sodium blockers, but it also plays a crucial role in disease modelling, where both fast (Pu and Boyden, 1997) and late (Coppini et al., 2013) sodium current are altered. An important feature of a model is its predictive power, and validation of a model using data not employed in model calibration is a central aspect of model credibility (Pathmanathan and Gray, Tomek et al. eLife 2019;8:e48890. DOI: https://doi.org/10.7554/eLife.48890 18 of 48 Research article Cell Biology Computational and Systems Biology Figure 9. Simulated and clinical 12-lead electrocardiogram. (A) 12-lead ECGs at 1 Hz: simulation using the ToR-ORd model in an MRI-based human torso-ventricular model (top) and a healthy patient ECG record (bottom, https://physionet.org, PTB database, subject 122; Bousseljot et al., 1995; Goldberger et al., 2000). (B) Electrode positions on the simulated torso. (C) Activation time map. (D) APD map. Tomek et al. eLife 2019;8:e48890. DOI: https://doi.org/10.7554/eLife.48890 19 of 48 Research article Cell Biology Computational and Systems Biology 2018;Carusi et al., 2012). With this in mind, we designed our study to first calibrate the developed model using a set of given criteria, with subsequent validation of the model using separate data that were not optimised for during development. The fact that ToR-ORd manifests a wide range of behaviours consistent with experimental studies, even though it was not optimised for these purposes, suggests its generality and a large degree of credibility. To facilitate future model development, we also created an automated ‘single-click’ pipeline, which evaluates a wide range of calibration and validation criteria and creates a comprehensive HTML report. New follow-up models can thus be immediately tested against criteria presented here, making it clear which features of the model are improved and/or deteriorated by any changes made. The greatest theoretical contribution of this work is the theory-driven reformulation of the L-type calcium current, namely the ionic activity coefficients and activation curve extraction. Activation curve of the current in previous cardiac models was based on the use of Nernst driving force in experimental studies, but the models then used Goldman-Hodgkin-Katz driving force to compute the current. This yields a theoretical inconsistency present in existing influential models of guinea pig, rabbit, dog, or human, for example (Luo and Rudy, 1994;Hund et al., 2008;O’Hara et al., 2011;Shannon et al., 2004;Grandi et al., 2010;Carro et al., 2011). We propose and demonstrate that in order to obtain consistent behaviour, the experimental I-V relationship measurements are to be normalised using the Goldman-Hodgkin-Katz driving force instead. Updated ionic activity coefficients and activation of the L-type calcium current improve key features of the current observed in the study underlying the ORd L-type calcium current model (Magyar et al., 2000), and strongly contribute to the improved reaction of the model to sodium blockade. The changes made are relevant in development of future models which use the Goldman-Hodgkin-Katz equation for L-type calcium current or other currents. A second major contribution of this work reaching beyond the model itself is the set of observations on modelling of I Kr , the dominant repolarising current in human ventricle. We noticed limitations of the ORd I Kr model, which may be a result of the single-pulse voltage clamp protocol to characterise the current behaviour. Approaches enabling the dissection of activation and recovery from inactivation based on more comprehensive experimental data, such as Lu et al. (2001) used in our work, may yield a more general and plausible model. In this study, this change was important predominantly for the response of the ventricular cell to calcium block, but our observations are highly relevant also for models of cells with naturally low plateau, such as Purkinje fibres or atrial myocytes. We anticipate that the main future development of the presented model will focus on the ryanodine receptor and the respective release from sarcoplasmic reticulum. Similarly to most existing cardiac models, the equations governing the release depend directly on the L-type calcium current, rather than on the calcium concentration adjacent to the ryanodine receptors, which is the case in cardiomyocytes. Future development of the ryanodine receptor model and calcium handling will extend the applicability of the model to other calcium-driven modes of arrhythmogenesis, such as delayed afterdepolarisations. Also, while the model represents to a certain degree the locality of I CaL calcium influx and calcium release via the utilization of the junctional calcium subspace, a more direct representation of local control (Stern, 1992;Hinch et al., 2004), realistic spatially distributed calcium handling (Colman et al., 2017), or representation of stochasticity, may improve the insights the model can give into calcium-driven arrhythmogenesis. However, we note that such changes (particularly the detailed distributed calcium handling) will increase computational cost of the model’s simulation. In addition, further research on the mechanisms regulating AP dependence on extracellular calcium concentration is needed to update this feature, not currently reproduced by most current human models (Passini and Severi, 2014). Acknowledgements We are grateful to Prof. Yoram Rudy, Prof. Derek Terrar, and Dr. Derek Leishman for very useful discussions. Tomek et al. eLife 2019;8:e48890. DOI: https://doi.org/10.7554/eLife.48890 20 of 48 Research article Cell Biology Computational and Systems Biology Additional information Funding Funder Grant reference number Author Wellcome 100246/Z/12/Z Blanca Rodriguez Wellcome 214290/Z/18/Z Blanca Rodriguez British Heart Foundation FS/17/22/32644 Alfonso Bueno-Orovio European Commission 675451 Blanca Rodriguez National Centre for the Replacement, Refinement and Reduction of Animals in Research NC/P001076/1 Blanca Rodriguez European Federation of Pharmaceutical Industries and Associations TransQST project (Innovative Medicines Initiative 2 Joint Undertaking 116030) Blanca Rodriguez BHF Centre of Research Excellence, Oxford RE/13/1/30181 Blanca Rodriguez UK National Supercomputing Archer RAP award (322 00180) Blanca Rodriguez Partnership for Advanced Computing in Europe AISBL 2017174226 Blanca Rodriguez Amazon Web Services Machine learning research award Blanca Rodriguez Horizon 2020 TransQST project (Innovative Medicines Initiative 2 Joint Undertaking 116030) Blanca Rodriguez The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication. Author contributions Jakub Tomek, Conceptualization, Software, Formal analysis, Investigation, Visualization, Methodology, Project administration; Alfonso Bueno-Orovio, Resources, Software, Formal analysis, Validation, Investigation, Methodology; Elisa Passini, Ana Minchole, Formal analysis, Investigation, Visualization; Xin Zhou, Software, Formal analysis, Investigation, Visualization; Oliver Britton, Laszlo Virag, Andras Varro, Data curation; Chiara Bartolucci, Investigation; Stefano Severi, Supervision, Investigation; Alvin Shrier, Methodology; Blanca Rodriguez, Conceptualization, Resources, Supervision, Funding acquisition, Methodology, Project administration Author ORCIDs Jakub Tomek https://orcid.org/0000-0002-0157-4386 Stefano Severi http://orcid.org/0000-0003-4306-8294 Decision letter and Author response Decision letter https://doi.org/10.7554/eLife.48890.sa1 Author response https://doi.org/10.7554/eLife.48890.sa2 Additional files Data availability No new experimental data were created. However, codes for simulations are available at https:// github.com/jtmff/torord (copy archived at https://github.com/elifesciences-publications/torord). Tomek et al. eLife 2019;8:e48890. DOI: https://doi.org/10.7554/eLife.48890 21 of 48 Research article Cell Biology Computational and Systems Biology References Beattie KA, Hill AP, Bardenet R, Cui Y, Vandenberg JI, Gavaghan DJ, de Boer TP, Mirams GR. 2018. Sinusoidal voltage protocols for rapid characterisation of ion channel kinetics. The Journal of Physiology 596:1813–1828. DOI: https://doi.org/10.1113/JP275733,PMID: 29573276 Bhattacharyya ML, Vassalle M. 1982. Effects of tetrodotoxin on electrical and mechanical activity of cardiac purkinje fibers. Journal of Electrocardiology 15:351–360. DOI: https://doi.org/10.1016/S0022-0736(82)81008-X, PMID: 6292321 Boukens BJ, Sulkin MS, Gloschat CR, Ng FS, Vigmond EJ, Efimov IR. 2015. Transmural APD gradient synchronizes repolarization in the human left ventricular wall. Cardiovascular Research 108:188–196. DOI: https://doi.org/10.1093/cvr/cvv202,PMID: 26209251 Bousseljot R, Kreiseler D, Schnabel A. 1995. Nutzung der EKG-Signaldatenbank CARDIODAT der PTB u ¨ber das internet. Biomedizinische Technik 40:317–318. DOI: https://doi.org/10.1515/bmte.1995.40.s1.317 Bradley E, Webb TI, Hollywood MA, Sergeant GP, McHale NG, Thornbury KD. 2013. The cardiac sodium current na(v )1.5 is functionally expressed in rabbit bronchial smooth muscle cells. American Journal of Physiology. Cell Physiology 305:C427–C435. DOI: https://doi.org/10.1152/ajpcell.00034.2013,PMID: 23784541 Britton OJ, Bueno-Orovio A, Van Ammel K, Lu HR, Towart R, Gallacher DJ, Rodriguez B. 2013. Experimentally calibrated population of models predicts and explains intersubject variability in cardiac cellular electrophysiology. PNAS 110:E2098–E2105. DOI: https://doi.org/10.1073/pnas.1304382110,PMID: 23690584 Britton OJ, Bueno-Orovio A, Vira´g L, Varro´ A, Rodriguez B. 2017. The electrogenic na + /K + Pump Is a Key Determinant of Repolarization Abnormality Susceptibility in Human Ventricular Cardiomyocytes: A PopulationBased Simulation Study. Frontiers in Physiology 8:278. DOI: https://doi.org/10.3389/fphys.2017.00278,PMID: 2 8529489 Bueno-Orovio A, Cherry EM, Fenton FH. 2008. Minimal model for human ventricular action potentials in tissue. Journal of Theoretical Biology 253:544–560. DOI: https://doi.org/10.1016/j.jtbi.2008.03.029,PMID: 18495166 Bueno-Orovio A, Hanson BM, Gill JS, Taggart P, Rodriguez B. 2012. In vivo human left-to-right ventricular differences in rate adaptation transiently increase pro-arrhythmic risk following rate acceleration. PLOS ONE 7: e52234. DOI: https://doi.org/10.1371/journal.pone.0052234,PMID: 23284948 Cardone-Noott L, Bueno-Orovio A, Minchole´ A, Zemzemi N, Rodriguez B. 2016. Human ventricular activation sequence and the simulation of the electrocardiographic QRS complex and its variability in healthy and intraventricular block conditions. EP Europace 18:iv4–iv15. DOI: https://doi.org/10.1093/europace/euw346 Carro J, Rodrı ´guez JF, Laguna P, Pueyo E. 2011. A human ventricular cell model for investigation of cardiac arrhythmias under hyperkalaemic conditions. Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 369:4205–4232. DOI: https://doi.org/10.1098/rsta.2011.0127 Carusi A, Burrage K, Rodrı ´guez B. 2012. Bridging experiments, models and simulations: an integrative approach to validation in computational cardiac electrophysiology. American Journal of Physiology-Heart and Circulatory Physiology 303:H144–H155. DOI: https://doi.org/10.1152/ajpheart.01151.2011,PMID: 22582088 Clerx M, Collins P, de Lange E, Volders PG. 2016. Myokit: a simple interface to cardiac cellular electrophysiology. Progress in Biophysics and Molecular Biology 120:100–114. DOI: https://doi.org/10.1016/j.pbiomolbio.2015. 12.008,PMID: 26721671 Colman MA, Pinali C, Trafford AW, Zhang H, Kitmitto A. 2017. A computational model of spatio-temporal cardiac intracellular calcium handling with realistic structure and spatial flux distribution from sarcoplasmic reticulum and t-tubule reconstructions. PLOS Computational Biology 13:e1005714. DOI: https://doi.org/10. 1371/journal.pcbi.1005714,PMID: 28859079 Coppini R, Ferrantini C, Yao L, Fan P, Del Lungo M, Stillitano F, Sartiani L, Tosi B, Suffredini S, Tesi C, Yacoub M, Olivotto I, Belardinelli L, Poggesi C, Cerbai E, Mugelli A. 2013. Late sodium current inhibition reverses electromechanical dysfunction in human hypertrophic cardiomyopathy. Circulation 127:575–584. DOI: https:// doi.org/10.1161/CIRCULATIONAHA.112.134932,PMID: 23271797 Coronel R, Janse MJ, Opthof T, Wilde AA, Taggart P. 2012. Postrepolarization refractoriness in acute ischemia and after antiarrhythmic drug administration: action potential duration is not always an index of the refractory period. Heart Rhythm 9:977–982. DOI: https://doi.org/10.1016/j.hrthm.2012.01.021,PMID: 22293142 Crumb WJ, Vicente J, Johannesen L, Strauss DG. 2016. An evaluation of 30 clinical drugs against the comprehensive in vitro proarrhythmia assay (CiPA) proposed ion channel panel. Journal of Pharmacological and Toxicological Methods 81:251–262. DOI: https://doi.org/10.1016/j.vascn.2016.03.009 Deb K. 2001. Multi-Objective Optimization Using Evolutionary Algorithms: An Introduction (Wiley-Interscience Series in Systems and Optimization. Wiley. Dhamoon AS, Jalife J. 2005. The inward rectifier current (IK1) controls cardiac excitability and is involved in arrhythmogenesis. Heart Rhythm 2:316–324. DOI: https://doi.org/10.1016/j.hrthm.2004.11.012,PMID: 15 851327 Drouin E, Charpentier F, Gauthier C, Laurent K, Le Marec H. 1995. Electrophysiologic characteristics of cells spanning the left ventricular wall of human heart: evidence for presence of M cells. Journal of the American College of Cardiology 26:185–192. DOI: https://doi.org/10.1016/0735-1097(95)00167-X,PMID: 7797750 Durrer D, van Dam RT, Freud GE, Janse MJ, Meijler FL, Arzbaecher RC. 1970. Total excitation of the isolated human heart. Circulation 41:899–912. DOI: https://doi.org/10.1161/01.CIR.41.6.899,PMID: 5482907 Dutta S, Minchole´ A, Zacur E, Quinn TA, Taggart P, Rodriguez B. 2016. Early afterdepolarizations promote transmural reentry in ischemic human ventricles with reduced repolarization reserve. Progress in Biophysics and Molecular Biology 120:236–248. DOI: https://doi.org/10.1016/j.pbiomolbio.2016.01.008,PMID: 26850675 Tomek et al. eLife 2019;8:e48890. DOI: https://doi.org/10.7554/eLife.48890 22 of 48 Research article Cell Biology Computational and Systems Biology Dutta S, Chang KC, Beattie KA, Sheng J, Tran PN, Wu WW, Wu M, Strauss DG, Colatsky T, Li Z. 2017a. Optimization of an in silico Cardiac Cell Model for Proarrhythmia Risk Assessment. Frontiers in Physiology 8: 616. DOI: https://doi.org/10.3389/fphys.2017.00616,PMID: 28878692 Dutta S, Minchole´ A, Quinn TA, Rodriguez B. 2017b. Electrophysiological properties of computational human ventricular cell action potential models under acute ischemic conditions. Progress in Biophysics and Molecular Biology 129:40–52. DOI: https://doi.org/10.1016/j.pbiomolbio.2017.02.007,PMID: 28223156 Engblom H, Foster JE, Martin TN, Groenning B, Pahlm O, Dargie HJ, Wagner GS, Arheden H. 2005. The relationship between electrical Axis by 12-lead electrocardiogram and anatomical Axis of the heart by cardiac magnetic resonance in healthy subjects. American Heart Journal 150:507–512. DOI: https://doi.org/10.1016/j. ahj.2004.10.041,PMID: 16169332 Franz MR, Swerdlow CD, Liem LB, 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. Journal of Clinical Investigation 82:972–979. DOI: https://doi.org/10.1172/ JCI113706,PMID: 3417875 Fu ¨lo ¨p L, Ba´nya´sz T, Magyar J, Szentandra´ssy N, Varro´ A, Na´na´ si PP. 2004. Reopening of L-type calcium channels in human ventricular myocytes during applied epicardial action potentials. Acta Physiologica Scandinavica 180: 39–47. DOI: https://doi.org/10.1046/j.0001-6772.2003.01223.x,PMID: 14706111 Garny A, Hunter PJ. 2015. OpenCOR: a modular and interoperable approach to computational biology. Frontiers in Physiology 6:26. DOI: https://doi.org/10.3389/fphys.2015.00026,PMID: 25705192 Gibor G, Yakubovich D, Peretz A, Attali B. 2004. External barium affects the gating of KCNQ1 potassium channels and produces a pore block via two discrete sites. The Journal of General Physiology 124:83–102. DOI: https://doi.org/10.1085/jgp.200409068,PMID: 15226366 Gima K, Rudy Y. 2002. Ionic current basis of electrocardiographic waveforms: a model study. Circulation Research 90:889–896. DOI: https://doi.org/10.1161/01.RES.0000016960.61087.86,PMID: 11988490 Goldberger AL, Amaral LAN, Glass L, Hausdorff JM, Ivanov PC, Mark RG, Mietus JE, Moody GB, Peng C-K, Stanley HE. 2000. PhysioBank, PhysioToolkit, and PhysioNet. Circulation 101:215–220. DOI: https://doi.org/10. 1161/01.CIR.101.23.e215 Gottlieb SS, Kukin ML, Medina N, Yushak M, Packer M. 1990. Comparative hemodynamic effects of procainamide, Tocainide, and encainide in severe chronic heart failure. Circulation 81:860–864. DOI: https:// doi.org/10.1161/01.CIR.81.3.860,PMID: 2106401 Grandi E, Pasqualini FS, Bers DM. 2010. A novel computational model of the human ventricular action potential and ca transient. Journal of Molecular and Cellular Cardiology 48:112–121. DOI: https://doi.org/10.1016/j. yjmcc.2009.09.019,PMID: 19835882 Guo D, Liu Q, Liu T, Elliott G, Gingras M, Kowey PR, Yan GX. 2011. Electrophysiological properties of HBI-3000: a new antiarrhythmic agent with multiple-channel blocking properties in human ventricular myocytes. Journal of Cardiovascular Pharmacology 57:79–85. DOI: https://doi.org/10.1097/FJC.0b013e3181ffe8b3,PMID: 20980921 Hinch R, Greenstein JL, Tanskanen AJ, Xu L, Winslow RL. 2004. A simplified local control model of calciuminduced calcium release in cardiac ventricular myocytes. Biophysical Journal 87:3723–3736. DOI: https://doi. org/10.1529/biophysj.104.049973,PMID: 15465866 Horvath B, Banyasz T, Jian Z, Hegyi B, Kistamas K, Nanasi PP, Izu LT, Chen-Izu Y. 2013. Dynamics of the late na (+) current during cardiac action potential and its contribution to afterdepolarizations. Journal of Molecular and Cellular Cardiology 64:59–68. DOI: https://doi.org/10.1016/j.yjmcc.2013.08.010,PMID: 24012538 Hund TJ, Decker KF, Kanter E, Mohler PJ, Boyden PA, Schuessler RB, Yamada KA, Rudy Y. 2008. Role of activated CaMKII in abnormal calcium homeostasis and I(Na) remodeling after myocardial infarction: insights from mathematical modeling. Journal of Molecular and Cellular Cardiology 45:420–428. DOI: https://doi.org/ 10.1016/j.yjmcc.2008.06.007,PMID: 18639555 Jost N, Vira´g L, Comtois P, Ordo ¨g B, Szuts V, Sepre´ nyi G, Bitay M, Kohajda Z, Koncz I, Nagy N, Sze´ l T, Magyar J, Kova´ cs M, Puska´s LG, Lengyel C, Wettwer E, Ravens U, Na´ na´si PP, Papp JG, Varro´ A, et al. 2013. Ionic mechanisms limiting cardiac repolarization reserve in humans compared to dogs. The Journal of Physiology 591:4189–4206. DOI: https://doi.org/10.1113/jphysiol.2013.261198,PMID: 23878377 Koller ML, Maier SK, Gelzer AR, Bauer WR, Meesmann M, Gilmour RF. 2005. Altered dynamics of action potential restitution and alternans in humans with structural heart disease. Circulation 112:1542–1548. DOI: https://doi.org/10.1161/CIRCULATIONAHA.104.502831,PMID: 16157783 Kristo ´f A, Husti Z, Koncz I, Kohajda Z, Sze´ l T, Juha´ sz V, Biliczki P, Jost N, Baczko´ I, Papp JG, Varro´ A, Vira´ g L. 2012. Diclofenac prolongs repolarization in ventricular muscle with impaired repolarization reserve. PLOS ONE 7:e53255. DOI: https://doi.org/10.1371/journal.pone.0053255,PMID: 23300901 Lancaster MC, Sobie EA. 2016. Improved prediction of Drug-Induced torsades de pointes through simulations of dynamics and machine learning algorithms. Clinical Pharmacology & Therapeutics 100:371–379. DOI: https:// doi.org/10.1002/cpt.367,PMID: 26950176 Legrand V, Vandormael M, Collignon P, Kulbertus HE. 1983. Hemodynamic effects of a new antiarrhythmic agent, flecainide (R-818), in coronary heart disease. The American Journal of Cardiology 51:422–426. DOI: https://doi.org/10.1016/S0002-9149(83)80073-3,PMID: 6823856 Linz KW, Meyer R. 2000. Profile and kinetics of L-type calcium current during the cardiac ventricular action potential compared in guinea-pigs, rats and rabbits. PflGers Archiv European Journal of Physiology 439:588– 599. DOI: https://doi.org/10.1007/s004240050982 Tomek et al. eLife 2019;8:e48890. DOI: https://doi.org/10.7554/eLife.48890 23 of 48 Research article Cell Biology Computational and Systems Biology Lu Y, Mahaut-Smith MP, Varghese A, Huang CL, Kemp PR, Vandenberg JI. 2001. Effects of premature stimulation on HERG K(+) channels. The Journal of Physiology 537:843–851. DOI: https://doi.org/10.1113/ jphysiol.2001.012690,PMID: 11744759 Luo CH, Rudy Y. 1994. A dynamic model of the cardiac ventricular action potential. II. afterdepolarizations, triggered activity, and potentiation. Circulation Research 74:1097–1113. DOI: https://doi.org/10.1161/01.RES. 74.6.1097,PMID: 7514510 Lyon A, Ariga R, Minchole´ A, Mahmod M, Ormondroyd E, Laguna P, de Freitas N, Neubauer S, Watkins H, Rodriguez B. 2018. Distinct ECG phenotypes identified in hypertrophic cardiomyopathy using machine learning associate with arrhythmic risk markers. Frontiers in Physiology 9:213. DOI: https://doi.org/10.3389/fphys.2018. 00213,PMID: 29593570 Magyar J, Iost N, Ko ¨rtve´ ly A ´., Ba´nya´ sz T, Vira´g L, Szigligeti P, Varro´ A, Opincariu M, Sze´ csi J, Papp JG, Na´na´si PP. 2000. Effects of endothelin-1 on calcium and potassium currents in undiseased human ventricular myocytes. Pflu ¨gers Archiv 441:144–149. DOI: https://doi.org/10.1007/s004240000400 Magyar J, Horva´th B, Va´czi K, Hegyi B, Go ¨nczi M, Dienes B, Kistama´ s K, Ba´ nya´ sz T, Baczko´ I, Varro´ A, Sepre´ nyi G, Csernoch L, Na´na´si PP, Szentandra´ssy N. 2017. Calcium activated chloride current in mammalian ventricular myocytes. Biophysical Journal 112:36a. DOI: https://doi.org/10.1016/j.bpj.2016.11.227 Mahajan A, Shiferaw Y, Sato D, Baher A, Olcese R, Xie LH, Yang MJ, Chen PS, Restrepo JG, Karma A, Garfinkel A, Qu Z, Weiss JN. 2008. A rabbit ventricular action potential model replicating cardiac dynamics at rapid heart rates. Biophysical Journal 94:392–410. DOI: https://doi.org/10.1529/biophysj.106.98160,PMID: 18160660 Makielski JC. 2016. Late sodium current: a mechanism for angina, heart failure, and arrhythmia. Trends in Cardiovascular Medicine 26:115–122. DOI: https://doi.org/10.1016/j.tcm.2015.05.006 Minchole ´A, Zacur E, Ariga R, Grau V, Rodriguez B. 2019. MRI-Based computational torso/Biventricular multiscale models to investigate the impact of anatomical variability on the ECG QRS complex. Frontiers in Physiology 10:1103. DOI: https://doi.org/10.3389/fphys.2019.01103,PMID: 31507458 Mortimer RG. 2008. The Activities of Nonvolatile Solutes. In: Physical Chemistry.274 3rd Edition. Elsevier. Muszkiewicz A, Britton OJ, Gemmell P, Passini E, Sa´nchez C, Zhou X, Carusi A, Quinn TA, Burrage K, BuenoOrovio A, Rodriguez B. 2016. Variability in cardiac electrophysiology: using experimentally-calibrated populations of models to move beyond the single virtual physiological human paradigm. Progress in Biophysics and Molecular Biology 120:115–127. DOI: https://doi.org/10.1016/j.pbiomolbio.2015.12.002,PMID: 26701222 Nagatomo T, January CT, Makielski JC. 2000. Preferential block of late sodium current in the LQT3 DeltaKPQ mutant by the class I(C) Antiarrhythmic flecainide. Molecular Pharmacology 57:101–107. PMID: 10617684 Nolasco JB, Dahlen RW. 1968. A graphic method for the study of alternation in cardiac action potentials. Journal of Applied Physiology 25:191–196. DOI: https://doi.org/10.1152/jappl.1968.25.2.191,PMID: 5666097 O’Hara T, Vira´g L, Varro´ A, Rudy Y. 2011. Simulation of the undiseased human cardiac ventricular action potential: model formulation and experimental validation. PLOS Computational Biology 7:e1002061. DOI: https://doi.org/10.1371/journal.pcbi.1002061 Okada J, Washio T, Maehara A, Momomura S, Sugiura S, Hisada T. 2011. Transmural and apicobasal gradients in repolarization contribute to T-wave genesis in human surface ECG. American Journal of Physiology. Heart and Circulatory Physiology 301:200–208. DOI: https://doi.org/10.1152/ajpheart.01241.2010,PMID: 21460196 Pappalardo F, Russo G, Tshinanu FM, Viceconti M. 2018. In silico clinical trials: concepts and early adoptions. Briefings in Bioinformatics 20:1699–1708. DOI: https://doi.org/10.1093/bib/bby043,PMID: 29868882 Passini E, Minchole´ A, Coppini R, Cerbai E, Rodriguez B, Severi S, Bueno-Orovio A. 2016. Mechanisms of proarrhythmic abnormalities in ventricular repolarisation and anti-arrhythmic therapies in human hypertrophic cardiomyopathy. Journal of Molecular and Cellular Cardiology 96:72–81. DOI: https://doi.org/10.1016/j.yjmcc. 2015.09.003,PMID: 26385634 Passini E, Britton OJ, Lu HR, Rohrbacher J, Hermans AN, Gallacher DJ, Greig RJH, Bueno-Orovio A, Rodriguez B. 2017. Human in Silico Drug Trials Demonstrate Higher Accuracy than Animal Models in Predicting Clinical Pro-Arrhythmic Cardiotoxicity. Frontiers in Physiology 8:668. DOI: https://doi.org/10.3389/fphys.2017.00668, PMID: 28955244 Passini E, Severi S. 2014. Computational analysis of extracellular calcium effect on action potential duration. Biophysical Journal 106:721a. DOI: https://doi.org/10.1016/j.bpj.2013.11.3985 Pathmanathan P, Gray RA. 2018. Validation and trustworthiness of multiscale models of cardiac electrophysiology. Frontiers in Physiology 9:6. DOI: https://doi.org/10.3389/fphys.2018.00106 Pitt-Francis J, Pathmanathan P, Bernabeu MO, Bordas R, Cooper J, Fletcher AG, Mirams GR, Murray P, Osborne JM, Walter A, Chapman SJ, Garny A, van Leeuwen IMM, Maini PK, Rodrı ´guez B, Waters SL, Whiteley JP, Byrne HM, Gavaghan DJ. 2009. Chaste: a test-driven approach to software development for biological modelling. Computer Physics Communications 180:2452–2471. DOI: https://doi.org/10.1016/j.cpc.2009.07.019 Pruvot EJ, Katra RP, Rosenbaum DS, Laurita KR. 2004. Role of calcium cycling versus restitution in the mechanism of repolarization alternans. Circulation Research 94:1083–1090. DOI: https://doi.org/10.1161/01. RES.0000125629.72053.95,PMID: 15016735 Pu J, Boyden PA. 1997. Alterations of na+ currents in myocytes from epicardial border zone of the infarcted heart. A possible ionic mechanism for reduced excitability and postrepolarization refractoriness. Circulation Research 81:110–119. DOI: https://doi.org/10.1161/01.res.81.1.110,PMID: 9201034 Pueyo E, Smetana P, Caminal P, deLuna AB, Malik M, Laguna P. 2004. Characterization of QT interval adaptation to RR interval changes and its use as a Risk-Stratifier of arrhythmic mortality in Amiodarone-Treated survivors of acute myocardial infarction. IEEE Transactions on Biomedical Engineering 51:1511–1520. DOI: https://doi.org/ 10.1109/TBME.2004.828050 Tomek et al. eLife 2019;8:e48890. DOI: https://doi.org/10.7554/eLife.48890 24 of 48 Research article Cell Biology Computational and Systems Biology Pueyo E, Husti Z, Hornyik T, Baczko´ I, Laguna P, Varro´ A, Rodrı ´guez B. 2010. Mechanisms of ventricular rate adaptation as a predictor of arrhythmic risk. American Journal of Physiology-Heart and Circulatory Physiology 298:H1577–H1587. DOI: https://doi.org/10.1152/ajpheart.00936.2009,PMID: 20207815 Pueyo E, Corrias A, Vira´g L, Jost N, Sze´ l T, Varro´ A, Szentandra´ssy N, Na´na´si PP, Burrage K, Rodrı ´guez B. 2011. A multiscale investigation of repolarization variability and its role in cardiac arrhythmogenesis. Biophysical Journal 101:2892–2902. DOI: https://doi.org/10.1016/j.bpj.2011.09.060,PMID: 22208187 Robinson P, Griffiths PJ, Watkins H, Redwood CS. 2007. Dilated and hypertrophic cardiomyopathy mutations in troponin and alpha-tropomyosin have opposing effects on the calcium affinity of cardiac thin filaments. Circulation Research 101:1266–1273. DOI: https://doi.org/10.1161/CIRCRESAHA.107.156380,PMID: 17932326 Scriven DR, Asghari P, Schulson MN, Moore ED. 2010. Analysis of Cav1.2 and ryanodine receptor clusters in rat ventricular myocytes. Biophysical Journal 99:3923–3929. DOI: https://doi.org/10.1016/j.bpj.2010.11.008, PMID: 21156134 Shannon TR, Wang F, Puglisi J, Weber C, Bers DM. 2004. A mathematical treatment of integrated ca dynamics within the ventricular myocyte. Biophysical Journal 87:3351–3371. DOI: https://doi.org/10.1529/biophysj.104. 047449,PMID: 15347581 Stern MD. 1992. Theory of excitation-contraction coupling in cardiac muscle. Biophysical Journal 63:497–517. DOI: https://doi.org/10.1016/S0006-3495(92)81615-6,PMID: 1330031 Streeter DD, Spotnitz HM, Patel DP, Ross J, Sonnenblick EH. 1969. Fiber orientation in the canine left ventricle during diastole and Systole. Circulation Research 24:339–347. DOI: https://doi.org/10.1161/01.RES.24.3.339, PMID: 5766515 Taggart P, Sutton PM, Opthof T, Coronel R, Trimlett R, Pugsley W, Kallis P. 2000. Inhomogeneous transmural conduction during early ischaemia in patients with coronary artery disease. Journal of Molecular and Cellular Cardiology 32:621–630. DOI: https://doi.org/10.1006/jmcc.2000.1105,PMID: 10756118 Taggart P. 2001. Transmural repolarisation in the left ventricle in humans during normoxia and ischaemia. Cardiovascular Research 50:454–462. DOI: https://doi.org/10.1016/S0008-6363(01)00223-1 Tomek J, Rodriguez B, Bub G, Heijman J. 2017. b-Adrenergic receptor stimulation inhibits proarrhythmic alternans in postinfarction border zone cardiomyocytes: a computational analysis. American Journal of Physiology-Heart and Circulatory Physiology 313:H338–H353. DOI: https://doi.org/10.1152/ajpheart.00094. 2017 Tomek J, Tomkova´ M, Zhou X, Bub G, Rodriguez B. 2018. Modulation of cardiac alternans by altered sarcoplasmic reticulum calcium release: a simulation study. Frontiers in Physiology 9:1306. DOI: https://doi.org/ 10.3389/fphys.2018.01306,PMID: 30283355 Tomek J. 2019. Repository for the cardiac model ToR-ORd. GitHub. 4ffab13. https://github.com/jtmff/torord Torres A, Torres D, Enriquez S, Leo´ n EPde, Dı ´az E. 2012. Evolutionary Multi-Objective Algorithms. In: Real-World Applications of Genetic Algorithms. In Tech. Tucker CR, Winkle RA, Peters FA, Harrison DC. 1982. Acute hemodynamic effects of intravenous encainide in patients with heart disease. American Heart Journal 104:209–215. DOI: https://doi.org/10.1016/0002-8703(82) 90194-6,PMID: 6808816 van Oosterom A, Hoekema R, Uijen GJ. 2000. Geometrical factors affecting the interindividual variability of the ECG and the VCG. Journal of Electrocardiology 33 Suppl:219–227. DOI: https://doi.org/10.1054/jelc.2000. 20356,PMID: 11265725 Weerapura M, Nattel S, Courtemanche M, Doern D, Ethier N, He´ bert TE. 2004. State-dependent barium block of wild-type and inactivation-deficient HERG channels in xenopus oocytes. The Journal of Physiology 526:265– 278. DOI: https://doi.org/10.1111/j.1469-7793.2000.t01-1-00265.x Weiss JN, Karma A, Shiferaw Y, Chen PS, Garfinkel A, Qu Z. 2006. From pulsus to pulseless: the Saga of cardiac alternans. Circulation Research 98:1244–1253. DOI: https://doi.org/10.1161/01.RES.0000224540.97431.f0, PMID: 16728670 Weiss JN, Garfinkel A, Karagueuzian HS, Chen PS, Qu Z. 2010. Early afterdepolarizations and cardiac arrhythmias. Heart Rhythm 7:1891–1899. DOI: https://doi.org/10.1016/j.hrthm.2010.09.017,PMID: 20868774 Woosley RL, Romero KA. 2015. QTdrugs list. CredibleMeds. https://www.crediblemeds.org Zacchia M, Abategiovanni ML, Stratigis S, Capasso G. 2016. Potassium: from physiology to clinical implications. Kidney Diseases 2:72–79. DOI: https://doi.org/10.1159/000446268 Tomek et al. eLife 2019;8:e48890. DOI: https://doi.org/10.7554/eLife.48890 25 of 48 Research article Cell Biology Computational and Systems Biology ToR-ORd). The first scenario applied an AP clamp based on ORd AP morphologies for control and for 50/50% I Na /I NaL block (as shown in Figure 3B). The second scenario considered AP clamps based on ToR-ORd AP morphologies in control and 50/50% I Na /I NaL block (as shown in Figure 3A). For each of the four models and the two AP clamp scenarios, the ratio of calcium transient amplitudes (sodium-block versus control) was computed (Appendix 1—table 1). This allowed quantifying the importance of differences in AP morphology and updates to the I CaL representation with regard to CaT amplitude changes under sodium block. Appendix 1—table 1. Dissection of improvement in reaction to sodium blockade. Ratios of calcium transient amplitudes for AP clamps corresponding to baseline and sodium-blocked models, arising from ORd and ToR-ORd. Row/column numbering used in the text does not include the header and the first column which describes the models (i.e. the text refers to four rows and two numerical columns). Model ORd AP clamps ToR-ORd AP clamps M1 (ORd) 1.366 1.061 M2 (ORd with ToR-ORd activity coefficients) 1.301 0.999 M3 (ORd with ToR-ORd activation curve) 1.326 1.015 M4 (ToR-ORd) 1.193 0.985 A value below one in Appendix 1—table 1 indicates the expected decrease in calcium transient amplitude with sodium block compared to control. This was achieved only by using the ToR-ORd AP clamp with the ToR-ORd model (row 4, column 2), and also but to a lesser extent with the ToR-ORd AP clamp and the ORd with ToR-ORd activity coefficients (row 2, column 2). Considering the ToR-ORd AP morphology (versus I CaL activation coefficients or activation curves) has the strongest effect in changes in calcium transient amplitude with sodium blockade. This is demonstrated by the universally lower values in the second column of Appendix 1—table 1 compared to the first column. The update of ionic activity coefficients alone leads to a reduced increase in CaT amplitude (in ORd), or a reduction in CaT amplitude (in ToR-ORd) following sodium blockade (Appendix 1—table 1, row 2 versus 1). Updated activation curve leads to a less pronounced increase in calcium transient amplitude in both models (Appendix 1—table 1, row 3 versus 1). Using the fully updated ToR-ORd model for simulations gave the most pronounced reduction in how much the calcium transient is increased upon sodium blockade using ORd AP clamps (Appendix 1— table 1, row 4 versus 1; column 1). It also yielded the most pronounced reduction in calcium transient amplitude following sodium blockade when using ToR-ORd AP clamps (Appendix 1—table 1, row 4 versus 1; column 2). AP morphology and peak AP levels regulate calcium transient amplitude through I CaL , and specifically through its voltage-dependency of activation curve, driving force and inactivation. It is important to note key differences between ToR-ORd and ORd models in this respect. In the ToR-ORd model, the AP morphology and the I CaL activation curve (Figure 2C) mean that I Na blockade leads to a reduction of peak I CaL activation. Conversely, in ORd, the activation curve is flat from 15 mV on (Figure 2C), and activation is not affected when sodium block reduces peak and early-plateau potential from 40 mV to 30–35 mV. At the same time, lowered peak and early-plateau potentials following sodium blockade increase I CaL driving force and weaken the voltage-driven inactivation, enhancing total I CaL . Furthermore, in the ToR-ORd model, AP with and without sodium block reach a near-identical notch and earlyplateau potential soon after the peak (Figure 3A). However, in the ORd, the difference in early-plateau membrane potentials lasts almost 30 ms after the peak (Figure 3B), prolonging the effect of increased driving force and reduced inactivation on I CaL , which is furthermore not compensated by reduction in activation. To further explore the impact of AP morphology on I CaL properties, I CaL activation and driving force was computed with the ToR-ORd model, again under four AP clamps considering AP morphologies shown in Figure 3A,B (obtained with ToR-ORd control and 50/ 50% I Na /I Na L block, and ORd control and 50/50% I Na /I NaL block). The results are presented in Appendix 1—figure 5. The activation with ToR-ORd AP clamp is reduced with sodium block Tomek et al. eLife 2019;8:e48890. DOI: https://doi.org/10.7554/eLife.48890 32 of 48 Research article Cell Biology Computational and Systems Biology versus control AP clamp morphologies (Appendix 1—figure 5A, yellow vs purple), given the difference in peak membrane potentials (Figure 3A) and their position on the activation curve (Figure 2C). However, the activation over time is similar with ORd AP morphology in control versus sodium block (Appendix 1—figure 5A, red vs blue), as the membrane potential is so high that full activation is reached in both cases. Appendix 1—figure 5. Effect of AP morphology on I CaL variables under sodium blockade. (A) I CaL activation. (B) ICaL driving force. In (C) are shown ratios of curves from (B) red pointwisedivided by blue (shown in light blue) and purple pointwise-divided by yellow (shown in green). Furthermore, the I CaL driving force obtained with the four AP clamps is shown in Appendix 1—figure 5B (more negative values correspond to a greater driving force), and Appendix 1—figure 5C displays the ratio of the driving forces between sodium-blocked AP clamp and control AP clamp for ORd and ToR-ORd morphologies. It is clear that the driving force elevation under sodium-block AP morphology versus control is much smaller and shorter lasting when clamping to ToR-ORd than ORd AP morphologies (Appendix 1—figure 5C). Therefore, the combination of ToR-ORd AP morphology and I CaL activation curve (which is consistent with experiments and with lower peak and early plateau voltages compared to ORd) result in reduction in both activation and driving force of I CaL with sodium bock. This explains how I Na block causes a smaller increase in I CaL and thus also in calcium transient amplitude using the ToR-ORd versus ORd models (Figure 3E and F). This smaller calcium increase due to I Na block is further compensated by the reduction in calcium entry and loading caused by I NaL block. The effect of I NaL on the simulated cell’s calcium loading is mediated by the NCX (reduced sodium influx via I NaL reduces calcium influx via NCX) and by I CaL (APD shortening induced by I NaL reduction also shortens I CaL , reducing calcium influx). Given that the effect of fixed-amount I NaL reduction on APD is stronger in ToR-ORd compared to ORd (Figure 3A,B), I NaL loss reduces calcium transient amplitude more in ToRORd. 11. Sodium block in fibre In order to assess how cell coupling affects the sodium block behaviour, we simulated halfblocks of I Na and I NaL in fibre (Appendix 1—figure 6). The effect of the respective blocks is generally consistent with single-cell behaviour (Figure 3) in that I Na block increases calcium transient amplitude, and I NaL block reduces it. The effect of both half-blocks combined shows the same trend as the single-cell but is of greater magnitude: ToR-ORd shows a greater reduction in the calcium transient amplitude, while ORd shows a greater increase. In this section, a version of ORd with updated I Na to allow for good propagation was used (Passini et al., 2016), and the tissue conductivity was set to achieve the conduction velocity of 63 m/s. Tomek et al. eLife 2019;8:e48890. DOI: https://doi.org/10.7554/eLife.48890 33 of 48 Research article Cell Biology Computational and Systems Biology Appendix 1—figure 6. Effect of sodium blockers on ToR-ORd and ORd in fibre. Membrane potential in fibre at four simulated conditions (control, half block of fast sodium, half block of late sodium, and half block of both fast and late sodium currents) in ToR-ORd (A) and ORd (B). In (C, D) is given the corresponding calcium transient. The measurements were taken in the node at the half of the fibre; consequently, action potentials in the conditions including I Na block are delayed to the other traces, given corresponding conduction delay. 12. Calibration: I CaL block and I Kr replacement During the model development (with the original ORd formulation of I Kr being used at the stage), we observed that in models with good AP morphology and response to sodium block, 50% I CaL block resulted in APD prolongation, which is in conflict with experimental evidence. Dissecting the cause of this phenomenon, we noticed considerably diminished I Kr density upon such treatment. Given that I Kr is the dominant repolarising current in human ventricle and human-specific computer models, this led to a net APD prolongation. The diminished I Kr density is due to the lack of activation of I Kr following from the combination of data-like AP morphology and formulation of fast and slow time constant of activation (Figure 3A in O’Hara et al., 2011). These time constants indicate a rapid activation between 25 and 40 mV, but when the membrane potential is reduced below 15 mV, activation time increases steeply (both fast and slow time constant are >1000 ms at 0 mV). Therefore, when the 50% I CaL block reduces plateau potential (a reduction of ca. 10 mV from control data-driven values of ca. 20 mV at 10 ms, 10 mV at 100 ms), I Kr activation is slowed greatly, producing a reduction in current density. Subsequently, we sought to understand why this problem was not present in the original ORd model, which was shown to respond well to I CaL block with regard to APD (O’Hara et al., 2011;Dutta et al., 2017a). This follows from the high-plateau AP configuration of ORd (ca. 38 mV at 10 ms and 20 mV at 100 ms); when such a high plateau is lowered by 10 mV, it still falls within the zone of rapid I Kr activation. Under such conditions, the loss of I CaL shortens the APD much more than the minor loss of I Kr prolongs it, with the net effect being APD shortening. Importantly, the article describing the Lu-Vandenberg model of I Kr (Lu et al., 2001) used in the final version of ToRORd also provides I Kr measurements under AP clamp. The study utilises AP clamps of both normal-plateau ventricular AP and low-plateau Purkinje AP (which can be taken as a proxy for I CaL blocked ventricular cell). The experimentally measured peak I Kr is slightly reduced in low-plateau AP clamp versus the high-plateau one (reduction to Tomek et al. eLife 2019;8:e48890. DOI: https://doi.org/10.7554/eLife.48890 34 of 48 Research article Cell Biology Computational and Systems Biology 70%). This is generally captured by simulations using our version of the Lu-Vandenberg model (reduction to 57%, Appendix 1—figure 7). However, peak I Kr in ORd reduces to mere 15% in the low-plateau AP clamp versus high-plateau one (Appendix 1—figure 7), which is due to the slow time constant of activation as described above. The fact that I Kr is capable of rapid activation (or recovery from inactivation, mimicking activation) is independently corroborated also by the elegant study of sinusoidal voltage pulse protocols by Beattie et al. (2018), where I Kr responds noticeably to membrane potential changes around 0 mV. Appendix 1—figure 7. AP clamps and resulting I Kr . (A,B) Endocardial and Purkinje AP used as an AP clamp (5th beat of the simulation shown). (C) I Kr produced by the endocardial AP. (D) I Kr produced by the Purkinje AP. Shown are data manually extracted from Lu et al. (2001), as well as simulations using ORd and ToR-ORd (L-V, coding for Lu-Vandenberg hERG model). A normalisation was performed in (C), (D) to facilitate the relative loss of I Kr using the lowplateau AP shape. In (C), the three traces were normalised to 0–1 range via division by maxima of the traces: max Data , max ORd , max LV . Subsequently, in (D), the data trace was divided by max Data , ORd simulation by max ORd , and the ToR-ORd model simulation by max LV . Thus, both traces produced by a given model (or measured experimentally) are divided by the same number, allowing the assessment of relative reduction of I Kr at the lowplateau AP. 13. Validation: drug block and APD In addition to four drugs and their effect on APD reported in the main manuscript Dutta et al. (2017b), also studied the effect of BaCl 2 , modelled as a 90% IK 1 block. However, we excluded considering this drug fundamentally due to the uncertainty on its multichannel effects on several potassium currents, in addition to I K1 . Specifically, it was previously shown that BaCl 2 also blocks I Kr and I Ks at the used concentration (Weerapura et al., 2004;Gibor et al., 2004). Furthermore, as shown in Kristo ´f et al. (2012), even a lower dose of BaCl 2 prolongs the AP plateau and early phase of repolarization, thus at membrane potentials at which IK 1 is not activated. This corroborates Tomek et al. eLife 2019;8:e48890. DOI: https://doi.org/10.7554/eLife.48890 35 of 48 Research article Cell Biology Computational and Systems Biology that BaCl 2 is likely to block I Kr and I Ks , in addition to I K1 , but to an uncertain degree. Furthermore, pure block of I K1 in the ORd model results in APD prolongation from around 50 mV only, which indicates a clear discrepancy with the experimental data. To further explore this issue, we simulated the effect of pure 90% I K1 block (which is consistent with the concentration of BaCl 2 used in the experiments) as well as additional I Kr block using the ToR-ORd. The results are shown in Appendix 1—table 2. Appendix 1—table 2. Effect of BaCl 2 on APD90 in experiments and simulations using the ToR-ORd model with different degrees of I K1 and I Kr block. Basic cycle length Experimental data (O’Hara et al., 2011) Simulated 90% I K1 block Simulated 90% I K1 + 15% I Kr block Simulated 90% I K 1 + 30% I Kr block 500 ms 271 (50) ms 250 ms 271 ms 297 ms 1000 ms 306 (50) ms 281 ms 307 ms 341 ms 2000 ms 361 (50) ms 286 ms 311 ms 342 ms For the three basic cycle length tested, the simulations are in overall agreement with the experimental data (O’Hara et al., 2011), with an improved match for 90% I K1 + 15% I Kr block. It is however uncertain (albeit likely) whether the block of other currents, such as for example I Ks contributes to explain the effect of BaCl 2 on the APD. A second matter discussed in this section is whether the assessment of nisoldipine in Figure 5D classifies as a calibration or validation result. With respect to I CaL block, the genetic algorithm optimisation contained a criterion promoting APD shortening in response to a 50% calcium blockade; however, the degree of block was different (90% when assessing nisoldipine), and the optimiser did not aim specifically for the APD reported by Dutta et al. (2017a). Additionally, the I Kr formulation substitution was motivated by the poor performance of the original I Kr model in response to 50% I CaL blockade. After the substitution was performed and the I Kr conductance scaled to achieve similar APD at 1 Hz pacing at control conditions, no further optimization towards the experimental data in Figure 5D was carried out. Specifically, at no time during development was the model adjusted to provide predictions close to Dutta et al. (2017a) in reaction to 90% I CaL . Therefore, APD shortening with I CaL blockade could be considered a partially calibration criterion, but quantitatively a validation one. 14. Validation: propagation at extreme hyperkalemia For extracellular potassium higher than 9 mM, a fibre was too inexcitable to support the propagation of a full AP. However, a partial activation was observed to propagate throughout the fibre, reaching ca. 30 mV and lasting ca. 50 ms (Appendix 1—figure 8). In case of localised hyperkalemia, such as during acute myocardial infarction, it is possible that even such low-peak wave might reactivate nonischemic cells upon exit of the ischemic zone, triggering re-entry and arrhythmia. Tomek et al. eLife 2019;8:e48890. DOI: https://doi.org/10.7554/eLife.48890 36 of 48 Research article Cell Biology Computational and Systems Biology Appendix 1—figure 8. Subthreshold fibre activation under extreme hyperkalemia. This trace was measured at the centre of the simulated fibre. 15.Simulation methods 15.1 Simulation details Cells were stimulated using a conservative potassium stimulus of 53 mV for 1 ms. Fixed concentrations in the model are: K o = 5 mM (extracellular potassium), Na o = 140 mM (extracellular sodium), Ca o = 1.8 mM (extracellular calcium), Cl o = 150 mM (extracellular chloride), Cl i = 24 mM (intracellular chloride). 15.1.1 Single-cell simulation protocols All results are shown at the end of a 1000-beat long pre-pacing train of stimuli. Fixed ionic concentrations were as stated above with two exceptions, where we matched concentrations to the respective studies (this was done both for ToR-ORd and ORd models when generating respective figures). The first exception is the replication of the IV relationship and steadystate inactivation of I CaL , based on Magyar et al. (2000), where Na o = 140.7 mM, K o = 5.4 mM, Cl o = 146.2, and Ca o = 2.5. Particularly, the calcium concentration is of importance, given that it affects the driving force of calcium and thus the IV relationship directly. In addition to the changed values of extracellular concentrations, we fixated the intracellular concentrations of sodium (to 0 mM), potassium (to 153 mM), and chloride (to 130 mM) to match concentrations in the patch-clamp pipette, given that in the long-term, diffusion would distribute the ions in approximately such way. We decided not to fix calcium levels in the model to the pipette concentrations (which was calcium-free), given that calcium entry via I CaL and the resulting release from the SR might last throughout the inactivating pulse, affecting I CaL function in a persistent way. A second exception is the protocol of induction of early afterdepolarisations (EADs) based on Guo et al. (2011). In this case, Na o = 137, Cl o = 148, and Ca o = 2. Third, for the P2/P1 protocol (Appendix 1-7), the concentrations were used as in Fu ¨lo ¨p et al. (2004): Cl i = 130 mM, Cl o = 157 mM, Na o = 144 mM, Ca o = 2.5 mM, K o = 5.6 mM, K i = 153 mM, Na i = 0 mM. Tomek et al. eLife 2019;8:e48890. DOI: https://doi.org/10.7554/eLife.48890 37 of 48 Research article Cell Biology Computational and Systems Biology 15.1.2 Hypertrophic cardiomyopathy and hyperkalemia representation 15.1.2.1 Hypertrophic cardiomyopathy Most of ionic remodelling was based on a comprehensive study by Coppini et al. (2013). All remodelling of HCM cells versus healthy cells was preferentially based on current density, protein expression if current density was not available, or mRNA expression if neither of the previous two measurements was available. The degree of remodelling was based on data means in the Coppini et al. study; where two subunits of a complex underlying a current were altered, an average of the two values corresponding to increase/reduction was used. Based on current density, I NaL was increased by 165%, I CaL by 25%, and I to was reduced by 80%. In addition, consistent with the study, the time constant of fast inactivation (both voltage and calcium-driven) was increased by 35%, and the time constant of slow inactivation was increased by 20%. Based on protein expression, NCX was increased by 50%, J up was reduced by 35%, and J rel by 30%. Based on mRNA, I K1 was reduced by 30%, I Kr by 35%, and I Ks by 55%. In addition, I NaK was reduced by 30% (Passini et al., 2016), and calcium affinity of troponin was increased slightly by reducing the K trp by 7% (Robinson et al., 2007). 15.1.2.2 Hyperkalemia Extracellular potassium level was raised from 5 mM with an increment of 0.5 mM up to 10 mM. Effective refractory period was determined using S1-S2 protocol, recording the highest S2 coupling interval duration which did not lead to propagation. Hyperkalemic conditions were simulated in a fibre, see below. 15.1.3 Fibre simulation The 1D fibre tissue simulations (hyperkalemia and Appendix 1-11) were conducted using the open-source software CHASTE. The homogeneous 1D tissue fibres had a length of 4 cm consisting of 200 nodes, with a conductivity of 1.64 mS/cm to achieve a diffusion coefficient of 1.171 cm 2 /s (Bueno-Orovio et al., 2008). All the fibres were paced for 20 beats before analysis. For the purpose of measurement of APD or effective refractory period, a stimulus was considered ‘propagating’ if the peak of membrane potential at the half-point of the fibre was above 0. Stimuli which did not propagate did not have their APD or refractory period measured. 15.1.4 Populations of models and drug safety assessment To assess the performance of the newly developed model for studies of populations of models, we constructed two populations of models using the ToR-ORd model as baseline, and the experimentally-calibrated population of models methodology (Britton et al., 2013;Muszkiewicz et al., 2016;Passini et al., 2017). A total of nine ionic conductances were randomly varied: fast and late Na + current (G Na and G NaL respectively), transient outward K + current (G to ), rapid and slow delayed rectifier K + current (G Kr and G Ks ), inward rectified K + current (G K1 ), Na + -Ca2 + exchanger (G NCX ), Na + -K + pump (G NaK ) and the L-type Ca2 + current (P Ca ). Two variability ranges were considered: [50-150]% and [0–200]% with respect to the baseline values, consistent with the reference study (Passini et al., 2017). The number of models passing calibration to experimental biomarkers was 2666 out of 3000 models generated for the [50-150]% population, and 817 out of 3000 for the [0–200]% population. To assess the suitability of the newly developed model for prediction of drug-induced pro-arrhythmic risk, we also replicated the results of our previous study on human in silico drug trials in population of human models (Passini et al., 2017), by testing the same 62 reference compounds at multiple concentrations on the [0–200%] population of models described above. Results were compared in terms of occurrence of drug-induced Tomek et al. eLife 2019;8:e48890. DOI: https://doi.org/10.7554/eLife.48890 38 of 48 Research article Cell Biology Computational and Systems Biology repolarisation abnormalities (RA) in the population and TdP score. CredibleMeds (Woosley and Romero, 2015) was used as gold standard for TdP risk, dividing the reference compounds in four categories: 1, high risk; 2, possible risk; 3, conditional risk; 4, safe (when not included in CredibleMeds). All the population of models described above were constructed using the Virtual Assay software (v3.2.1119, 2019 Oxford University Innovation) (Passini et al., 2017). Further analysis of the results was performed in Matlab (Mathworks Inc Natwick, MA). 15.1.5 Transmurality, 3D simulation, and pseudo-ECG measurement The formulation of transmurality (representation of endocardial, midmyocardial, and epicardial cells) was retained as in O’Hara et al. (2011), with two exceptions. First, the increase of I CaL conductance from endocardium to midmyocardium was reduced to 2-fold, which corresponds to the mean of the data (O’Hara et al., 2011, Figure 10). Second, the increase of I to from endocardium to both midmyocardium and epicardium was reduced to 2fold to avoid excessive notch following the peak of the AP. Membrane potential traces of the three cell types at 1 Hz pacing are given in Appendix 1—figure 9. Appendix 1—figure 9. Action potential morphology for three cell subtypes at 1 Hz pacing. The propagation of the electrical activity in the human ventricles was modelled using the heart bidomain equations and solved with the Chaste software (Pitt-Francis et al., 2009). A cardiac magnetic resonance (CMR)-informed torso ventricular model was used. The ventricular element size was set to 0.4 mm to ensure numerical convergence of the finite element software Chaste for electrophysiological simulations (Dutta et al., 2016;PittFrancis et al., 2009). The presented ToR-ORd model was used to represent the membrane kinetics. Transmural and apex to base cell electrophysiological heterogeneities based on experimental and clinical data from Okada et al. (2011);Drouin et al. (1995);Taggart (2001);Boukens et al. (2015) were incorporated in the biventricular heart model. Transmural heterogeneities were modelled using three layers as in Lyon et al. (2018) of endocardial (45% of the transmural width), mid-myocardial (25%) and epicardial cells (30%) with different AP properties as shown above. Apex-to-base heterogeneities were modelled by including a gradual increase of I Ks conductance from base to apex resulting in APD differences of 25 ms. Tomek et al. eLife 2019;8:e48890. DOI: https://doi.org/10.7554/eLife.48890 39 of 48 Research article Cell Biology Computational and Systems Biology Heart fibre directions were generated using the Streeter rule-based method (Streeter et al., 1969), and tissue conductivities were set to generate realistic conduction velocities for the ToR-ORd cell model, numerical scheme, and mesh resolution. The longitudinal intracellular conductivity of 1.64 mS/cm in a 1D fibre model resulted in a physiological conduction velocity of ~65 cm/s (Taggart et al., 2000). As CardoneNoott et al., 2016 used 1.5 mS/cm, intracellular orthotropic conductivities and the axisymmetric extracellular conductivities were obtained by using a scaling factor of 1.0934 to the ones in Cardone-Noott et al., 2016. Sinus rhythm was simulated at 1 Hz using a phenomenological activation model with early endocardial activation sites and a fast endocardial layer representing a tightly-packed endocardial Purkinje network (Cardone-Noott et al., 2016) and adapted to the ventricular geometry (Minchole ´et al., 2019). Simulation of the pseudo-ECG was computed by calculating the extracellular potentials in the virtual standard 12-lead electrode positions (Lyon et al., 2018) from the ventricular geometry following the dipole model (Gima and Rudy, 2002) as: feð Þ¼ZW DrVm r 1 re k k   dr where e¼ex;ey;ez  are the electrode position coordinates, Dis the diffusion tensor, and Vmis the membrane potential. The integral is calculated over the whole myocardium volume, W. 15.1.6 Timely killing of crashing simulations during multiobjective GA We found it crucial to modify the behaviour of @ode15s in our simulations, using Matlab ODE Events to limit 1 s of simulated time to 5 s of runtime (with the runtime being <0.5 s in normal conditions), killing a simulation upon exceeding the limit. When parameters of the models are randomly perturbed during GA fitting, it is possible to achieve a combination which leads to model instability and eventual crash. However, as ode15s attempts to reduce time step in such a case, it takes up to 6 hr to crash, blocking CPU cores and stalling the whole generation. With the timely killing of unstable simulations, we could run 30 generations with population size of 2500 in ca. 30 hr on an Azure virtual machine with 64 cores, using parallel fitness evaluation. 15.2 Evaluation pipeline and HTML reporter In addition to the model code itself, we also provide a model evaluation pipeline. This allows simulation of the single-cell calibration and validation criteria using a standard model function (I.e. all except I CaL properties such as I-V relationship, or disease models, where modified model codes are used.). After the simulation of a comprehensive set of protocols, an HTML report is produced (Appendix 1—figure 10), providing a clickable mapping between criteria and figures showing how the model fulfils them. Icons adjacent to the criteria allow rapid assessment of whether the model behaves well in the given criterion. Any model with appropriate interface of inputs and outputs can be simulated. It is particularly easy to simulate variants of the ToR-ORd model, for example with changed conductances of currents, as the same structure of parameters that is passed to the model simulation wrapper can be passed to this evaluation pipeline. Given that the evaluation code and subsequent report generation are fully automated, it is easy to extend the pipeline to compute and visualise other evaluation criteria. Tomek et al. eLife 2019;8:e48890. DOI: https://doi.org/10.7554/eLife.48890 40 of 48 Research article Cell Biology Computational and Systems Biology Appendix 1—figure 10. Evaluation pipeline. The left menu shows a list of evaluation criteria with the adjacent icon indicating whether they are fulfilled by the model. The right pane contains figures visualising the model performance in the criteria. The figures can be either viewed by scrolling or by clicking the corresponding entry in the left menu. The report header contains a timestamp of the report generation and a slider which controls the size of figures in the right pane. The whole report is written in HTML, as automatically generated from Matlab code, and can be viewed in any internet browser. The whole process of report generation (and storage of included plots) can be run separately for calibration and validation criteria (with separately stored figures), facilitating unbiased model development, where the model can be assessed using the calibration criteria only, without the user observing any validation results. It is nevertheless possible to also generate a report for only validation criteria, or both calibration and validation. Below are given the criteria evaluated, along with how the rating of the model is defined (PASS = good behaviour in a given criterion, FAIL = problematic behaviour, NA = feature not automatically compared to data, usually because of difficulty of accurate rating; this code is also used when a plot is only an illustration, such as examples of action potential shape at different pacing frequencies). Unless specified otherwise, FAIL is given when PASS is not fulfilled. 15.2.1 Calibration criteria 1. Action potential morphology. PASS ~AP from 10 to 500 ms is within the 10–90% quantile range of the Szeged dataset. The first 10 ms are ignored so as not to limit the model in its peak potential, as this is strongly modulated by cell coupling (the available AP data are based on small-tissue samples). 2. Sample traces showing the AP and calcium transient at different pacing frequencies. Always NA, this is just an illustrator figures useful for an eyeballing-type of assessment. Tomek et al. eLife 2019;8:e48890. DOI: https://doi.org/10.7554/eLife.48890 41 of 48 Research article Cell Biology Computational and Systems Biology Jrel ¼1:5378 1rel;CaMK  Jrel;NP þrel;CaMK Jrel;CaMK 15.3.16 J up Jup;NP ¼0:005425 Cai 0:00092 þCai Jup;P¼2:75 0:005425 Cai 0:00092 0:00017 þCai Jleak ¼0:0048825 CaNSR 15 Tomek et al. eLife 2019;8:e48890. DOI: https://doi.org/10.7554/eLife.48890 48 of 48 Research article Cell Biology Computational and Systems Biology