1. Cell Biology
  2. Computational and Systems Biology
Download icon

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

  1. Jakub Tomek  Is a corresponding author
  2. Alfonso Bueno-Orovio
  3. Elisa Passini
  4. Xin Zhou
  5. Ana Minchole
  6. Oliver Britton
  7. Chiara Bartolucci
  8. Stefano Severi
  9. Alvin Shrier
  10. Laszlo Virag
  11. Andras Varro
  12. Blanca Rodriguez  Is a corresponding author
  1. University of Oxford, United Kingdom
  2. University of Bologna, Italy
  3. McGill University, Canada
  4. University of Szeged, Hungary
Research Article
Cite this article as: eLife 2019;8:e48890 doi: 10.7554/eLife.48890
19 figures and 3 tables


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 (IKr). The transition rates are given in Appendix 1-15.3.5.

Action potential, calcium transient, and ICaL 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 ToR-ORd 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). ICaL,tot is the sum of currents corresponding to all ions (Ca2+, 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 ICaL,tot is shown to demonstrate that the current reversal is not due to other ions. We note that the difference in total amplitude of ICaLin Figure 2F follows predominantly from different action potential shape in ToR-ORd vs ORd, consistent with the I-V relationship.

Sodium blockade in ToR-ORd and ORd.

Simulated effect of sodium current block on the action potential (A, B) and calcium transient (C,D) using the ToR-ORd (left) and ORd (right) models. Control simulations are shown as blue traces, whereas results for 50% block of INa and INaL are shown as red traces. Panels (E,F) show values of changes in calcium transient amplitude with respect to control (‘rel. CaTA’) for varying degrees of INa/INaL block for ToR-ORd and ORd respectively (one is no change, values > 1 correspond to an increase the calcium transient amplitude).

EADs and alternans.

(A) Experimentally observed EADs produced in nondiseased human myocytes exposed to 0.1 µM dofetilide (∼85% IKr block) paced at 0.25 Hz (Guo et al., 2011) and the corresponding simulation using the ORd model (O'Hara et al., 2011), figure reproduced as allowed by the CC-BY license). See Appendix 1-15.1.1 for ionic concentrations used. (B) An EAD produced by the ToR-ORd model at the same conditions. (C,E) APD and calcium transient amplitude versus pacing base cycle length (bcl); the bifurcation indicates alternans. (D,F) Example of APD and calcium transient alternans at the bcl of 260 ms.

Drug block and APD.

All four panels contain mean and standard deviation of experimental data (black) as measured in O'Hara et al. (2011) for three basic cycle lengths (bcl), predictions of the Dutta optimised dynamic-IKr ORd model (blue, Dutta et al., 2017a), and the predictions of the ToR-ORd model (red). (A) 1 µM E-4031 (70% IKr block). (B) 1 µM HMR-1556 (90% IKs block). (C) 10 µM mexiletine (54% INaL, 9% IKr, 20% ICaL block). (D) 1 µM nisoldipine (90% ICaL block). Drug concentrations and their effects on channel blocks are based on Dutta et al. (2017b).

APD accommodation and S1-S2 protocol.

(A) APD accommodation measured experimentally (Franz et al., 1988) and in simulation of the ORd model (reproduced as allowed by the CC-BY licence from O'Hara et al., 2011). (B) APD accommodation in ToR-ORd. (C) S1-S2 restitution curve (S1 = 1000 ms) in ToR-ORd fibre, ORd fibre (including INa modification to facilitate propagation Passini et al., 2016) and human tissue samples data (O'Hara et al., 2011).

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 [50-150]% (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 drug-induced 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).

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% IKr block induces EADs in HCM cell, but not in a healthy one.

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.

Appendix 1—figure 1
Diagram of how activation curves are extracted from the I-V relationship.

The activation curves (top right) are obtained by dividing values in the I-V curve (top left) at each pulse potential by the driving force considering either V-ECa-based driving force (bottom right) or the Goldman-Hodgkin-Katz (GHK) flux equation (bottom left).

Appendix 1—figure 2
Effect of the proposed changes in the activation curve and activation coefficient on the I-V relationship of ICaL in ORd versus ToR-ORd.
Appendix 1—figure 3
Goldman-Hodgkin-Katz equation with colour coding of components.
Appendix 1—figure 4
P2/P1 protocol for two pulse durations.
Appendix 1—figure 5
Effect of AP morphology on ICaL variables under sodium blockade.

(A) ICaL activation. (B) ICaL driving force. In (C) are shown ratios of curves from (B) red pointwise-divided by blue (shown in light blue) and purple pointwise-divided by yellow (shown in green).

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 INa block are delayed to the other traces, given corresponding conduction delay.

Appendix 1—figure 7
AP clamps and resulting IKr.

(A,B) Endocardial and Purkinje AP used as an AP clamp (5th beat of the simulation shown). (C) IKr produced by the endocardial AP. (D) IKr 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 IKr using the low-plateau AP shape. In (C), the three traces were normalised to 0–1 range via division by maxima of the traces: maxData, maxORd, maxLV. Subsequently, in (D), the data trace was divided by maxData, ORd simulation by maxORd, and the ToR-ORd model simulation by maxLV. Thus, both traces produced by a given model (or measured experimentally) are divided by the same number, allowing the assessment of relative reduction of IKr at the low-plateau AP.

Appendix 1—figure 8
Subthreshold fibre activation under extreme hyperkalemia.

This trace was measured at the centre of the simulated fibre.

Appendix 1—figure 9
Action potential morphology for three cell subtypes at 1 Hz pacing.
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.


Table 1
Criteria and human-based studies used in ToR-ORd calibration and validation.
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)
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)
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).

ModelORd AP clampsToR-ORd AP clamps
M1 (ORd)1.3661.061
M2 (ORd with ToR-ORd activity coefficients)1.3010.999
M3 (ORd with ToR-ORd activation curve)1.3261.015
M4 (ToR-ORd)1.1930.985
Appendix 1—table 2
Effect of BaCl2 on APD90 in experiments and simulations using the ToR-ORd model with different degrees of IK1 and IKr block.
Basic cycle lengthExperimental data (O'Hara et al., 2011)Simulated 90% IK1 blockSimulated 90% IK1 + 15% IKr blockSimulated 90% IK1 + 30% IKr block
500 ms271 (±50) ms250 ms271 ms297 ms
1000 ms306 (±50) ms281 ms307 ms341 ms
2000 ms361 (±50) ms286 ms311 ms342 ms

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).

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

Downloads (link to download the article as PDF)

Download citations (links to download the citations from this article in formats compatible with various reference manager tools)

Open citations (links to open the citations from this article in various online reference manager services)