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

1Human-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. They exploit mature human-based modelling and simulation technology to perform virtual testing of pharmacological therapies or devices.

2Human 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. Its key strengths are the 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. Consequently, ORd was selected by a panel of experts as the model best suited for regulatory purposes.

3Most 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 and in data from additional studies using human cardiomyocytes. 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. 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). 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 in-depth re-evaluation of the L-type calcium current (ICaL) 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 (IKr), the dominant repolarisation current in human ventricle, under conditions that reflect experimental data-driven plateau potentials.

4Using 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.

5Table 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 ToR-ORd 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 IK1 block, based on evidence in a range of studies summarised in. 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 1-1). The validation criteria include response to rate changes, drug action and disease, to explore the predictive power of the model under clinically-relevant conditions.

6We initially performed the evaluation of the ORd model 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 IK1 block. The results are later demonstrated in Figures 2 and 3, and Methods: Calibration of IK1 block and resting membrane potential. Secondly, we attempted parameter optimisation using a multiobjective genetic algorithm. 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.

7We then proceeded to reevaluate the ionic current formulations based on experimental data and biophysical knowledge. Key currents included ICaL and specifically its driving force and activation, as well as the INa, IKr, IK1 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.

8Details 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 was used for comparison with the ToR-ORd model.

9The 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.

10The ICaL 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.

11The ICaL 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-Eion), 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, the driving force of ions via ICaL in cardiac models is modelled based on the Goldman-Hodgkin-Katz (GHK) flux equation. The driving force based on the GHK equation is:where z is the charge of the given ion, V is the membrane potential, F,R,T are conventional thermodynamic constants, and [S]i, [S]o are intracellular and extracellular activities of the given ionic specie., where is the ionic activity coefficient and the concentration (in either the intracellular or extracellular space, yielding or ).

12In order to compute the ionic driving force via the GHK equation, it is necessary to know the ionic activity coefficients of the intracellular (γi) and extracellular (γo) space. The Luo-Rudy model and other Rudy-family models use γo = 0.341 for extracellular space and γi = 1 for the intracellular space. Models based on the Shannon model use 0.341 for both intracellular and extracellular space, but we were unable to find the motivation for this change.

13The Debye-Hückel theory is commonly used to compute the activity coefficients. We used the Davies equation, which extends the basic Debye-Hückel equation to be accurate for ionic concentrations found in living cells:where A is a constant (~0.5 for water at 25°C, ~0.5238 at 37°C), zi is the charge of the respective ion, and I is the ionic strength of the solution. The ionic strength is defined as:where mi is the concentration of the i-th ionic specie present. For concentrations in a study measuring properties of ICaL, I is ca. 0.15-0.17. This warrants the use of Davies equation, which was shown to be accurate for I up to 0.5, unlike the basic Debye-Hückel equation, which is accurate for I up to 0.01 only.

14We 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.

15Throughout 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 γi = 1 in the Luo-Rudy model is unclear, as by the Davies (or by any Debye-Hückel variant) equation, I would have to be zero, which is possible only when there are no ions present.

16An additional improvement in the ICaL 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 ICaL 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.;;;;; ): whereas the Nernstian driving force of (V-ECa) is used to derive the activation curve, the GHK driving force is then used to calculate ICaL. Indeed, experimental studies reporting the activation curve of ICaL generally use the Nernstian driving force of (V-ECa) with ECa being the experimentally measured reversal potential of approximately 60 mV. This is explicitly stated in, and also the activation curve by used in the ORd model is consistent with dividing the IV relationship with (V-60).

17In 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 ICaL was obtained by dividing the I-V curve from 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. The following capped Gompertz function (a flexible sigmoid) was found to be the best fit to the resulting steady-state activation curve:where V is the membrane potential.

1820% of ICaL was placed in the main cytosolic space, consistent with the literature. This increases the plateau-supporting capability of ICaL, given that the myoplasmic ICaL is subject to a weaker calcium-dependent inactivation than ICaL in the junctional subspace. Other minor changes are given in Appendix 1-15.3.3.

19The 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 IKr formulation was used (further details in Appendix 1-12). ICaL block is known to shorten APD experimentally 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 IKr activation, which is however not consistent with relevant experimental data. We therefore considered alternative IKr formulations and specifically the Lu-Vandenberg Markov model (Figure 1B). The Lu-Vandenberg IKr 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 ICaL block in a model which manifests experimental data-like plateau potentials and which contains the ORd IKr 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 IKr formulations with experimental data, demonstrating the good agreement with experimental data of the Lu-Vandenberg formulation but not the ORd.

20Following the inclusion of the Lu-Vandenberg IKr formulation, all models generated during model calibration exhibited APD shortening in response to ICaL block.

21The INa current formulation was replaced by an alternative human-based formulation, given established limitations of the original model with regards to conduction velocity and excitability, comment on article from 05 Oct 2012). The Grandi INa model was updated to account for CaMKII phosphorylation (Appendix 1-15.3.1).

22Also from the Grandi model, we added the calcium-sensitive chloride current I(Ca)Cl and background chloride current IClb formulation. 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. The motivation to add these currents was to facilitate the shaping of post-peak AP morphology (via I(Ca)Cl), with IClb 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 IK1 block as explained in the next subsection.

23The IK1 model was replaced with the human-based formulation by, as it was shown to be key for simulations of hyperkalemic conditions. The IK1 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.

24When evaluating the baseline ORd model against the selected criteria, we observed that a reduction in IK1 results in hyperpolarisation of the cell (from −88 to −88.16 mV at 1 Hz pacing). However, it is established that IK1 reduction depolarises cells experimentally. 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 IK1 block, consistent with experimental data.

25We applied a multiobjective genetic algorithm (MGA, @gamultiobj function in Matlab, ) 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 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.

26To 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 (; 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/.

27We 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 and OpenCOR). Furthermore, the Myokit library enables conversion of the CellML file to other languages (such as C or Python).

28The 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 (, Figure 6) and isolated human cardiomyocytes. We note that compared to the Szeged-ORd dataset, our model manifests a slightly increased peak membrane potential in the single-cell form, similar to single-cell experimental data. 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.

29Both 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, 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).

30As described in Materials and methods, the ToR-ORd ICaL activation curve was extracted from experimental data, using the Goldman-Hodgkin-Katz formulation of ionic driving force, ensuring theoretical consistency, unlike the ORd ICaL 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 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 ICaL steady-state inactivation curve provides a good fit to the same experimental data. This difference follows from the formulation of calcium-dependent inactivation of ICaL (see Appendix 1-5 for details).

31We observed that in cases of elevated ICaL (e.g. in midmyocardial cells), ORd reverses current direction towards positive values, which is an unexpected behaviour given its reversal potential of 60 mV. Conversely, the ToR-ORd model yields negative ICaL values in such conditions, consistent with it being an inward current (Figure 2F). This is a direct consequence of the updates to the extracellular/intracellular calcium activity coefficients (as explained in Appendix 1-6), which supports their credibility and it is important for cases of elevated ICaL, such as under ß-adrenergic stimulation.

32We have also simulated a P2/P1 protocol as measured experimentally by, where two rectangular pulses are applied with varying interval between them. Both ORd and ToR-ORd qualitatively agree with the experimental data (Appendix 1-7).

33Figure 3A-D illustrates AP and calcium transient changes caused by block of sodium currents in ToR-ORd (left) and ORd (right). As sodium blockers act on channel Nav1.5 mediating both the fast (INa) and late (INaL) sodium current, we simulate the effect of combined partial INa and INaL block. The ToR-ORd model manifests a small reduction in calcium transient amplitude (Figure 3C), unlike ORd, which gives a sizeable increase (Figure 3D); ToR-ORd is thus consistent with the observed negative inotropy of sodium blockers. This is a major improvement in the ToR-ORd model, as sodium current reduction is involved in a range of disease conditions in addition to pharmacological block.

34Experimental evidence shows that the ratio of INa and INaL block is drug and dose-dependent, with INaL usually being blocked more than INa (Appendix 1-9). Figure 3E,F illustrates the change in calcium transient amplitude obtained with the ToR-ORd and ORd models, respectively, for several combinations of INa and INaL availability. Both models show a similar general trend where reduced INa availability increases calcium transient amplitude and reduced INaL availability diminishes it; however, the models differ strongly in relative contributions of these components. The ToR-ORd model shows negative inotropy for almost all combinations of blocks. A mild increase in inotropy may be achieved only under near-exclusive INa block. Conversely, ORd shows a general tendency for increased calcium transient amplitude; a reduction occurs only when the sodium current block targets near-exclusively INaL. ToR-ORd presents a greater calcium transient amplitude reduction than ORd in response to INaL block, as the current has a greater role in indirect modulation the cell’s calcium loading via APD change. At the same time, ToR-ORd shows a much smaller calcium transient amplitude increase in response to INa block than ORd because of the updated ICaL activation curve (Figure 2C), as well as closer-to experimental data AP morphology (Figure 2A) and its effect on ICaL. A detailed explanation is given in Appendix 1-10.

35Fibre simulations carried out to assess the effect of cell coupling on the effect of sodium block are consistent with the single-cell simulations (Appendix 1-11). The difference in response to half-block of INa and INaL between ToR-ORd and ORd is even larger, as ToR-ORd in fibre predicts a greater reduction in calcium transient amplitude than in single cell (−14% vs −6% respectively), while ORd in fibre predicts a slightly greater increase in calcium transient amplitude than in single cell (+25% vs +24%).

36With the ToR-ORd model, the 14% reduction in CaT amplitude in the electrically coupled fibre with 50% block of both INa and INaL is generally consistent with clinical data on sodium blockers: Encainide reduced stroke work index by 15% and cardiac index by 8%. In another study using encainide, the cardiac index was reduced by 18% and the stroke volume index by 28%. Flecainide reduced left ventricular stroke index by 12% and the left ventricular ejection fraction by 9%. Simulations with the ToR-ORd model show overall agreement with the clinical data. However, a direct quantitative comparison is challenging given the different indices of contractility measured (CaT amplitude versus clinical indices) and that it is not possible to estimate the exact ratios of INa and INaL block in clinical data (Appendix 1-9).

37EADs are an important precursor to arrhythmia, manifesting as a membrane potential depolarisation during late plateau and/or early repolarisation. They are thought to arise mainly from ICaL current reactivation. The ToR-ORd model manifests EADs at conditions used experimentally in nondiseased human endocardium (; Figure 4A). The amplitude of simulated EADs is 14 mV (Figure 4B), which matches the maximum EAD amplitude shown by. We also note that the experimental data by Guo et al. manifest early plateau potential of ca. 23 mV (which is matched by ToR-ORd), in line with other studies we referred to previously regarding this matter.

38Repolarisation alternans is another established precursor to arrhythmia, facilitating the formation of conduction block. It is induced by rapid pacing and it is mostly thought to arise from calcium transient amplitude oscillations being translated to APD oscillations, although purely voltage-driven mechanism was also proposed. Alternans in the ToR-ORd model is calcium-driven and appears via the same mechanism as in the ORd: sarcoplasmic reticulum calcium cycling refractoriness. It occurs at rapid pacing, in both calcium and APD (Figure 4C–F). The peak APD alternans amplitude (difference in APD between consecutive beats) is 12 ms, which is matches the value 11 ± 2 ms reported in human hearts without a structural disease. Direct quantitative comparison is however slightly limited by the fact that the data were recorded in RV septum, which may or may not differ from endocardial cells in alternans amplitude.

39Figure 5 illustrates simulations of drug action using the ToR-ORd model (red traces), compared to experimental data (black traces) and to simulations with the ORd model reparametrised by (blue dashed lines). APD is shown in the presence of IKr block (E-4031, Figure 5A), IKs block (HMR-1556 Figure 5B), multichannel block of INaL, ICaL, IKr (mexiletine, Figure 5C), and a ICaL block (nisoldipine, Figure 5D), at base cycle lengths of 500, 1000, and 2000 ms. We note that while the Dutta et al. model was specifically optimised for response of APD to these drug blocks, no such treatment was applied to the ToR-ORd model, making the results presented here an independent validation. Appendix 1-13 contains further details on the choice and use of the drug data.

40The predictions produced by the ToR-ORd model are in good agreement with experimental data, particularly given the lack of optimisation towards this result. Simulating E-4031, ToR-ORd provides a prediction similar to the experimental data mean and the Dutta model (Figure 5A). This is crucial, given the key role of IKr in the repolarisation reserve of human cardiomyocytes. The response to IKs blockade via HMR-1556 is even better in ToR-ORd than in the Dutta model, which is also within standard deviation of the data, but carries a clear trend towards AP prolongation (Figure 5B). When simulating the multichannel blocker mexiletine, ToR-ORd prediction is within standard deviation of the experimental data, with the Dutta model giving similar or closer-to-mean predictions at 0.5 and 1 Hz (Figure 5C). The predicted effect of the calcium blocker nisoldipine in the ToR-ORd model matches well the experimental data mean (Figure 5D), even better than the Dutta model (also within standard deviation). We note that the good performance of the simulated nisoldipine effect critically relies on the IKr replacement (Materials and methods and Appendix 1-12).

41Experimental measurements in human cardiomyocytes show how the APD shortens upon increase in pacing frequency, and then prolongs again, as the pacing frequency returns to control (Figure 6A, top). APD adaptation dynamics with changes in heart rate are regulated by changes in sodium homeostasis, and their manifestation in QT adaptation have been shown to be useful for arrhythmia risk prediction. While simulations with the ORd model capture the general trend of APD accommodation, there are differences compared to the experimental data (Figure 6A). First, changes in pacing rate are followed by slow-dynamics (~30 s) APD prolongation not present in the experimental recordings. Second, the time constant of accommodation is generally slow. Conversely, the ToR-ORd model reproduces the pattern of accommodation well, where the change in APD soon after change in frequency is relatively fast, and then gradually slows down (Figure 6B). This suggests that the ionic balance in ToR-ORd is likely to have been improved compared to ORd.

42A second indicator of how a model responds to a change in pacing frequency is the S1-S2 restitution protocol. The S1-S2 restitution curve obtained with the ToR-ORd model is given in Figure (Figure 6C), showing a good agreement with the experimental data.

43Drug safety testing is one of the key applications of computer modelling which has yielded highly promising results. To assess the suitability of ToR-ORd for drug safety testing, we replicated the study by, which was carried out using populations of models based on the ORd model. Two populations were created based on ToR-ORd similarly to the original study, altering conductances of important currents within the ranges of 50–150% and 0–200%. Models in both populations are stable under significant perturbation of ionic conductances, which supports the robustness of the model (Figure 7A).

44Prediction of the risk of drug-induced Torsades de Pointes based on simulated drug-induced repolarisation abnormalities using ToR-ORd population yielded similar results to the original study, with predicted risk being correct for 54 out of 62 compounds (87% accuracy). Compared to, the assessment of Mexiletine (a predominantly sodium blocker that is safe) was improved from false positive to true negative. High-dose Mexiletine led to formation of many EADs in ORd, but not in ToR-ORd (Figure 7B), highlighting the importance of the advances on sodium blockers presented in this work. At the same time, Procainamide and Metrodinazole were misclassified as false negatives compared to. However, these drugs are controversial, as Metrodinazole is considered non-torsadogenic by, and this study predicted both the drugs to be non-risky. Torsadogenic risk for all evaluated compounds and the confusion matrix of the classification are given in Figure 7C.

45Hyperkalemia, the elevation of extracellular potassium, is a hallmark of acute myocardial ischemia caused by the occlusion of coronary artery. It was shown that hyperkalemia can significantly inhibit sodium channel excitability following repolarisation, leading to the prolongation of postrepolarisation refractoriness. The dispersion of effective refractory periods (ERPs) between normal and ischemic zones forms a substrate for the initiation of re-entrant arrhythmia. In this new model, we tested the effect of hyperkalemia on tissue excitability using 1D fibres. As shown in Figure 8A, the elevation of extracellular potassium level led to an increase of the resting membrane potential (RMP) and the decrease of AP upstroke amplitude. As a result of weaker upstroke and more depolarised RMP, the APD shortened under hyperkalemia; however, the ERPs were prolonged due to the stronger sodium channel inactivation caused by the elevation of RMP (Figure 8B). Therefore, this new model successfully reproduced the longer post-repolarization refractoriness under hyperkalemia observed in experiments, and it can be used in the simulations of re-entrant arrhythmia under acute ischemia. In this regard, it presents an improvement over the original ORd model, which did not manifest postrepolarisation refractoriness without further modifications.

46Hypertrophic cardiomyopathy (HCM) is among the most common cardiomyopathies, manifesting as abnormal thickening of the cardiac muscle without an obvious cause. Beyond 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 ) to our baseline model using an approach similar to, 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. Ultimately, the HCM variant of our model is more prone to the formation of EADs (Figure 8E), as was shown experimentally. This difference is in line with postulated key role of ICaL and NCX in EAD formation, 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.

47We 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.

48Figure 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. 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. The APD map shows longer APD in the endocardium and the base, and shorter APDs in the epicardium and the apex, respectively (Figure 9D).

49In 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 IKr. 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.

50The AP morphology of ToR-ORd is in agreement with the Szeged endocardial myocyte dataset used to construct the state-of-the art ORd model. 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. The validation of ToR-ORd shows that the model responds well to drug block with regards to APD. Good APD accommodation (reaction to abrupt, but persisting changes in pacing frequency) indicates a good balance between ionic currents. Replication of arrhythmia precursors such as early afterdepolarisations and alternans 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 and hypertrophic cardiomyopathy. The model is also shown to be promising in drug safety testing, and whole-heart simulations demonstrate physiological conduction velocity and produce a plausible ECG signal. Among the improved behaviours compared to the state-of-the-art ORd model, 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, 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 and late sodium current are altered.

51An 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. 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.

52The 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. 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, 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.

53A second major contribution of this work reaching beyond the model itself is the set of observations on modelling of IKr, the dominant repolarising current in human ventricle. We noticed limitations of the ORd IKr 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 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.

54We 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 ICaL calcium influx and calcium release via the utilization of the junctional calcium subspace, a more direct representation of local control, realistic spatially distributed calcium handling, 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.