Modelling and Simulation for Preclinical Cardiac Safety Assessment of Drugs with Human iPSC-Derived Cardiomyocytes

Abstract

As a potentially life threatening side effect, pharmaceutical compounds may trigger cardiac arrhythmias by impeding the heart’s electrical and mechanical function. For this reason, any new compound needs to be tested since 2005 for its proarrhythmic risk both during the preclinical and the clinical phase of the drug development process. While intensive monitoring of cardiac activity during clinical tests with human volunteers constitutes a major cost factor, preclinical in vitro tests with non cardiac cells and in vivo tests with animals are currently under serious debate because of their poor extrapolation to drug cardiotoxicity in humans. For about five years now, regulatory agencies, industry and academia are working on an overhaul of the cardiac drug safety paradigm that is built a) on human heart muscle cells, that can be abundantly bioengineered from donor stem cells without ethical concerns (human induced pluripotent stem cell derived cardiomyocytes, hiPSC-CMs), and b) on computational models of human cardiac electrophysiology both at the cellular and the organ level. The combined use of such human in vitro and human in silico models during the preclinical phase is expected to improve proarrhythmia test specificity (i.e. to lower the false-positive rate), to better inform about the need of thorough heart monitoring in the clinic, and to reduce or even replace animal experiments. This review article starts by concisely informing about the electrical activity of the human heart, about its possible impairment due to drug side effects, and about hiPSC-CM assays for cardiac drug safety testing. It then summarizes the mathematical description of human cardiac electrophysiology in terms of mechanistic ODE and PDE models, and illustrates how their numerical analysis may provide insight into the genesis of drug induced arrhythmias. Finally, this paper surveys proarrhythmic risk estimation methods, that involve the simulation of human heart muscle cells, and addresses opportunities and challenges for future interdisciplinary research.

1 Introduction

Heart muscle cells have the ability to generate and conduct electrical signals referred to as action potentials. The concerted production and propagation of these action potentials throughout the heart forms the basis of its rhythmic contraction and relaxation. This rhythmic activity can be adversely affected by pharmaceutical compounds such that the blood pumping function of the heart is greatly diminished and a life-threatening situation occurs. In fact, drug induced cardiac arrhythmias belong to the most serious drug side effects, have led to the post-marketing withdrawal of several medicines and are a major reason for the discontinuation of candidate compounds during the drug development process. Consequently, the assessment of the proarrhythmic risk is an integral and mandatory part of both preclinical and clinical safety testing of new pharmaceutical compounds. In that regard, a major problem to the industry is that preclinical tests based on animal models show a significant false positive rate and hence may lead to the development stop of otherwise promising compounds based on rather non-specific markers of proarrhythmia.

Mathematical modelling of cardiac electrophysiology [235], [64] looks back on more than 50 years of history and has become a prime example of computational systems biology. While nonlinear ODE systems form the established tool for describing action potentials of single heart muscle cells, complex PDE systems of the reaction diffusion type model the propagation of the electrical excitation wave through interconnected heart muscle cells, cardiac tissue or the whole heart. Recently, the field gained additional momentum driven by international and multidisciplinary research initiatives [36], [98] that aim at overcoming the above mentioned non-specificity problem of preclinical drug cardiotoxicity testing. As one pillar of these efforts, in vitro data about drug effects on transmembrane ionic currents shall be integrated into a mechanistic action potential model of a native human heart muscle cell [138]. This model shall then serve as the basis for proarrhythmic risk prediction by means of machine learning techniques. Finally, these predictions shall be scrutinised by in vitro experiments with human stem cell derived heart muscle cells (hiPSC-CMs), that keep the genetic features of their human donor and provide - without ethical constraints - an abundant stock of human material for drug safety testing [34].

The overall goal is to better support the decision whether intensive (and hence expensive) electrocardiography (ECG) monitoring [262] is still necessary during late-stage clinical trials in humans. This decision might be further facilitated by means of in silico drug trials with 3D high-fidelity multiscale and multiphysics models of the human heart [234], [10], [207]. In return, the respective computational techniques also offer the opportunity to design and analyze experiments with 2D monolayers or 3D structures engineered from hiPSC-CMs as, e.g., macroscale human myocardia [240] or human minihearts in a jar [136].

Several reviews of computational methods for cardiac safety testing of drugs have been previously published [154], [270], [48], [75]. One distinctive feature of the article at hand is that it contrasts ODE and PDE modelling of drug induced cardiac arrhythmias with arrhythmic events that can also be experimentally observed in hiPSC-CMs. Furthermore, this paper covers the association of cardiac arrhythmias with bifurcations and pattern formation in dynamical systems, and surveys simulation based approaches for the classification of drugs according to their proarrhythmic liability. Section 2 gives a short introduction to cardiac electrophysiology and presents how its distortion by pharmaceutical compounds may lead to cardiac arrhythmias. Section 3 informs about the Comprehensive In Vitro Proarrhythmia Assay (CiPA) initiative, that is led by the US Food and Drug Administration (FDA) and works on defining a new safety paradigm, and illustrates how hiPSC-CMs are currently used in drug cardiotoxicity testing. Section 4 deals with the mathematical formulation of cardiac action potentials and their spatial propagation through cardiac tissue and organ. Furthermore, it addresses the rich dynamical repertoire of these models that corresponds with diverse types of arrhythmic behaviour both at the cellular and at the tissue level. Section 5 discusses methods that integrate mechanistic modelling of drug effects into statistical learning tools for the assessment of the proarrhythmic risk. Finally, Sect. 6 lists challenges related to the predictive ability of risk classifiers and to the credibility of cardiac electrophysiology models, in particular of hiPSC-CMs, and also touches on some topics that so far are hardly considered in simulation based drug cardiotoxicity testing. For the convenience of the reader, an appendix gives basic examples of cardiac AP models in full detail and features two prominent logistic regression models for proarrhythmic risk classification.

2 Drug Induced Cardiac Arrhythmias

2.1 Electrical Activity of the Heart

The human heart consists of billions of muscle cells (cardiomyocytes, CMs), that rhythmically contract in order to provide the mechanical force for pumping blood throughout the body [161]. The rhythmic contraction of the heart is regulated by a wave of electrical activation, which is based on the excitability of cardiac cells and their ability to generate and propagate action potentials. An action potential (AP) is a characteristic transient change in the transmembrane voltage that starts with a depolarization phase and ends with a repolarization phase, see Fig. 1 (left). Pacemaker cells located in the sinoatrial (SA) node regularily depolarize by themselves and generate spontaneous APs, a property that is referred to as automaticity. In contrast, atrial and ventricular CMs are quiescent, have a resting state and can generate an AP only if the transmembrane voltage is raised above a certain threshold value by an external stimulus. Still, the CMs are electrically connected via gap junctions, which are ion conducting channels formed by proteins. Hence, the currents created during an AP of a cell may flow to and stimulate an AP in the neighbouring cells. Once a cell has produced an AP, there is a refractory period during which the cell can generate no further AP. Together, excitability and refractoriness of cardiac tissue ensure that the wavefront of depolarization, initiated by the pacemaker cells, is propagated in a unidirectional manner from the atria through the atrioventricular (AV) node to the ventricles.

Fig. 1

left: Typical AP profiles of CMs from different regions of the human heart and their link to the components of a normal ECG signal. Figure credit [161] (reproduced with permission). right: ECG traces corresponding with a normal sinus rhythm (top) as well as life threatening TdP (middle) and VF (bottom) arrhythmias. Figure credit United Medical Education (reproduced with permission)

The propagation of APs throughout the heart causes changes in the potential on the body surface, which can be recorded in an electrocardiogram (ECG). Figure 1 (right, top) displays an ECG that represents a normal sinus rhythm of about 70 beats per minute for a resting person. The small P wave, see Fig. 1 (left), reflects the depolarization of the atria, and the larger QRS complex corresponds with the depolarization of the ventricles. The end of the T wave indicates their full repolarization, such that the QT interval is a measure of the ventricular AP duration (APD), see Fig. 2 (left).

Fig. 2

left: Correspondence of the QT interval of the ECG with the ADP of a ventricular CM. Figure credit [61] (reproduced with permission). right: AP regulation by voltage gated ion channels and their corresponding ionic currents into and out of the cell. Figure credit [196] (reproduced with permission)

The AP of a CM is orchestrated by voltage gated ion channels, which are membrane bound proteins that open and close in dependence on the transmembrane voltage [84]. These ion channels ensure a selective and chronological transmembrane flow of ions such as calcium, potassium and sodium into and out of the cell, see Fig. 2 (right). Of particular importance is the voltage gated $\text{K}_{\text{V11.1}}$ potassium channel, that carries the rapidly activating delayed rectifier potassium current $\text{I}_{\text{Kr}}$ and is encoded by the human ether a go-go related gene (hERG) [221]. Mutations of hERG may impair ventricular repolarisation and prolong the QT interval, a condition referred to as congenital long QT syndrome (LQTS) [44].

When the calcium channels are opened during an AP, calcium flows into the cell and triggers the quick release of calcium from the sarcoplasmic reticulum (SR) into the cytosol [15]. The rapid increase of cytosolic calcium leads to the contraction of the myofilaments, and the subsequent relaxation of the cell is enabled by a pump that drives calcium back into the SR. Oscillatory calcium dynamics not only are involved in the excitation-contraction coupling, but also are implicated in the automaticity of the pacemaker cells [277]. Accordingly, the pacemaking function is maintained by two coupled oscillatory mechanisms termed membrane clock, formed by the ionic currents across the membrane, and calcium clock. As part of the latter, the SR spontaneously releases calcium at local sites within the cell, which is in contrast to the global calcium release triggered by the AP.

2.2 Cardiac Arrhythmias

A cardiac arrhythmia is a pathophysiological deviation from the normal sinus rhythm and can be caused by perturbations in the excitation of the tissue or in the conduction of the excitation wave. Figure 1 (right, middle) displays an ECG that corresponds with a Torsade de Pointes (TdP) arrhythmia, which is a particular form of polymorphic ventricular tachycardia (VT) with repeated twisting of the QRS complexes. TdP may degenerate into ventricular fibrillation (VF), see Fig. 1 (right, bottom), which is a life threatening arrhythmia as it can lead to sudden cardiac death within minutes after onset. The basic electrophysiological mechanisms of cardiac arrhythmias [7], [248] can be roughly devided into a) focal activity, that arises in a localized tissue region due to abnormal automaticity or triggered activity of CMs, and b) reentry activity, where an action potential wave fails to terminate itself and reexcites a region of nonrefractory tissue.

An example of a premature trigger is given by early afterdepolarizations (EADs), which are abnormal depolarizations during the repolarization phase of the AP, see Fig. 3 (left). EADs can lead to ectopic beats outside the region of the sinoatrial node and may initiate polymorphic VT and TdP [265], [267]. EADs emerge in a condition referred to as reduced repolarization reserve [215], which is an ionic current imbalance due to reduction of outward and/or an increase of inward currents. Another example for abnormal activity is given by repolarization alternans, which are beat-to-beat oscillations in APD and/or calcium transients. Cellular alternans may manifest as T-wave alternans in the ECG and are involved in arrhythmogenesis including TdP formation [205].

Fig. 3

Electrophysiological mechanisms of cardiac arrhythmias. left: Early afterdepolarizations (EADs) are voltage oscillations during repolarization and may result from an imbalance of ionic currents, e.g., caused by an inhibition of the hERG current $\text{I}_{\text{Kr}}$. Figure credit [60] (reproduced with permission). right: Spiral waves of excitation underlie reentry activity and, depending on the pattern of movement, are associated with different types of cardiac arrhythmias. Figure credit [102] (reproduced under CC-BY license)

Examples of reentry activity are formed by two-dimensional spiral and three-dimensional scroll waves, where action potentials travel in circular motion around an anatomical obstacle and reactivate their site of origin. Such spiral movement is also possible in absence of a structural obstacle, in which case the wave rotates around a functional obstacle formed by its own refractory tail, see Fig. 3 (right). Spiral waves may become self-sustained sources of high-frequency excitation and can drive atrial and ventricular tachycardia and fibrillation. In particular, spiral waves that rotate along a circular trajectory are linked to monomorphic VT, while spiral waves that meander through space are related with polymorphic VT and TdP, see Fig. 3 (right). As a result of conduction block, spiral waves may also break into several spiral waves that continuously terminate and redevelop in a chaotic manner, and such spiral wave breakups are considered as a possible basis for the transition from VT to VF [202]. Recently, it was demonstrated that also a single spiral wave can display chaotic tip trajectories that are consistent with fibrillation [143].

2.3 Drug Cardiotoxicity

Pharmaceutical compounds - not only cardiovascular drugs but also chemotherapeutics, antipsychotics, antibiotics and others - may interfere with the function of the heart and trigger life-threatening cardiac arrhythmias such as TdP [278], see Fig. 1 (right, middle). In the 1990s, a sequence of reports about the occurence of QT interval prolongation and TdP - presumably or demonstrably induced by drugs and in some cases with fatal outcome - led to the withdrawal of 14 drugs from the market [232]. By that time, it was already known that patients with congenital LQTS are predisposed to TdP [263]. As a number of drugs were found to block the hERG channel, see Fig. 4 (left), a mechanistic link [221] between the congenital (or inherited) LQTS and the acquired (or drug induced) LQTS was finally established. This further supported the association between drug induced LQTS and TdP [278].

Fig. 4

left: hERG channel blockade by a compound. Figure credit [257] (reproduced with permission). right: Causal relationship between hERG channel blockade and TdP. Figure credit [216] (reproduced with permission)

According to the trigger-substrate model [105], APD and QT interval prolongation not only facilitate the occurence of EADs but also promote heterogeneity of repolarization across the tissue. A disrupted spatial dispersion of repolarization may form a substrate which maintains reentrant arrhythmias once initiated by EADs or other triggers. Based on this paradigm, illustrated in Fig. 4 (right), the International Conference on Harmonization (ICH) in 2005 released two mandatory guidelines for preclinical (ICH S7B) and clinical (ICH E14) cardiac safety testing of new pharmaceutical compounds, see [2] and [1]. ICH S7B first prescribes an in vitro assay in order to assess the potency of the test compound to block the hERG channel heterologously expressed in host cells. The compound is considered as a torsadogenic hERG blocker if its hERG $\text{IC50}$ value is below the 30-fold of its effective free therapeutic plasma concentration (EFTPC) [211], see also Fig. 25 (right) in the appendix. $\text{IC50}$ denotes the drug concentration at which the current through the channel is reduced by 50$\%$. The second component of ICH S7B is an in vivo QT assay in animal models in order to assess whether the test compound induces QT prolongation. With respect to human trials, ICH E14 specifies a thorough QT study in order to assess a test compound’s impact on the heart rate corrected QT (QTc) interval of healthy volunteers during clinical phase I. This test rates the compound as possibly torsadogenic if the mean prolongation of QTc is above 10 ms. The ICH S7B and ICH E14 guidelines are still in use today and have proven to be successful - no drug that accordingly was approven to market has ever since been withdrawn for cardiac safety concerns.

However, since 2005 more than a billion dollar was spent on hundreds of intensive ECG monitorings in clinical phases IIb and III due to safety concerns raised during ICH S7B and ICH E14 testing [190]. Furthermore, it is meanwhile well established that, on the one hand, QT prolongation not necessarily leads to TdP [125], and on the other hand, hERG blockers not necessarily prolong QT as they may have compensating effects on other ion channels. Hence, hERG blockade and QT interval prolongation are no longer considered to be specific markers to assess a drug’s proarrhythmic liability. As a negative consequence, potentially benefitial compounds may be discontinued at an early stage, mostly based on their impact on hERG alone rather than due to actual proarrhythmic effects.

More complex preclinical assays include ventricular wedge preparations and isolated Langendorff-perfused hearts model from canine or rabbit [125], [105]. While they may grant more informative data about a drug’s proarrhythmic liability, they only offer low throughput recordings and do not bypass the error-prone animal-to-human extrapolation. Although also adult human primary cardiomyocytes, isolated from non-transplantable human donor hearts, form a preclinical model for proarrhythmia risk testing [167], [5], ethical concerns and limited availability do currently not render them suitable for high-throughput studies.

3 Cardiac Safety Testing with Human IPSC-Derived Cardiomyocytes

In contrast, human iPSC-CMs provide easily accessible human test material without ethical concerns and hence have become a promising tool for cardiac arrhythmia research [86]. In the context of preclinical cardiac safety testing of drugs, their potential is currently explored in multidisciplinary research projects run by, e.g., the Japanese IPS Cardiac Safety Assessment (JiCSA) consortium [98] and the Comprehensive in vitro Proarrhythmia Assay (CiPA) initiative [218], [36], [61], [70], see Fig. 5 (left). Within CiPA, hiPSC-CM technology shall be combined with computational modelling of human ventricular CMs [138], [139] see Sects. 4 and 5 for details, in order to better assess a drug’s torsadogenic liability before going into clinical trials. For that purpose, 28 reference compounds were classified based on expert opinion into low, medium and high proarrhythmic risk groups and divided into 12 training and 16 validation compounds, see Table 1. These CiPA compounds include the non-cardiovascular drugs terfenadine and cisapride, which are hERG blockers withdrawn from the market due to cardiac toxicity, but also the non-torsadogenic drug verapamil, whose hERG blocking effects are compensated by a simultaneous block of the calcium channel.

Fig. 5

left: The four pillars of the CiPA initiative, figure credit [262] (reproduced under CC-BY license). right: Ventricular-like APs of a spontaneously beating single hiPSC-CM recorded in manual patch clamp, data by courtesy of Drug Discovery Sciences at Boehringer Ingelheim, Biberach

Table 1 The 28 CiPA training and validation compounds [61]

Human iPSC-CMs possess a heterogeneous AP phenotype characterized by spontaneous activity [99], see Fig. 5 (right) for an example of a ventricular-like AP from a single hiPSC-CM. For studying drug cardiotoxicity, multicellular hiPSC-CM settings are more suitable and range from 2D monolayers over 3D organoids to engineered tissue [135] as well as miniature heart chambers [136]. The associated assays [34] comprise the recording of electrical intra- and extracellular signals, calcium transients, beating movements and contractile forces. For instance, Fig. 6 (left) displays local extracellular APs (LEAPs) with drug induced EADs, while Fig. 6 (right) shows extracellular field potentials (FPs) and cell motion images from 3D human iPS cell-engineered cardiac tissue during drug induced reentrant acticity [102], [227], compare also with Fig. 3. Two recent studies [16] and [17] confirmed the utility of hiPSC-CMs for proarrhythmic risk prediction by showing that high and intermediate CiPA risk compounds consistently evoke EADs-like arrhythmic events across different electrophysiological readout platforms and different commercially available hiPSC-CM lines, see Fig. 7. Note that according to a scoring system suggested in [108], the occurence of EADs in hiPSC-CMs immediately attributes very high proarrhythmic risk to a compound, regardless of all other endpoints measured.

Fig. 6

Examples of experimental modelling of drug induced cardiac arrhythmias with hiPSC-CMs. left: Drug induced EADs in hiPSC-CMs, LEAP data by courtesy of Axion Biosystems. right: Drug induced reentrant activity in 3D human iPS cell-engineered heart tissue. Figure credit [102] (reproduced under CC-BY license)

Fig. 7

Different types of drug induced arrhythmic events in hiPSC-CMs. left: Extracellular FP recordings from microelectrode assays (MEA) and fluorescence signals from voltage sensitive optical platforms (VSO). Figure credit [16] (reproduced with permission). right: Simultaneous recordings of extracellular FPs and impedance as a surrogate for contractility. Figure credit [17] (reproduced with permission)

4 Modelling and Simulation of Cardiac Electrophysiology and Arrhythmias

4.1 Modelling of Cardiac APs

Computational models of single cell APs typically take the form of nonlinear systems of ordinary differential equations (ODEs). Based on viewing the cell membrane as a capacitor with capacitance $C_{m}$, these equations describe the time course of the transmembrane voltage $V$ according to

$$ C_{m} \frac{dV}{dt} = - \sum \limits _{\text{ion}} I_{\text{ion}} \quad \left ( + I_{\text{stim}} \right ), $$

(1)

along with the dynamics of the voltage dependent ionic transmembrane currents $I_{\text{ion}}$ into and out of the cell, compare with Fig. 2 (right). Furthermore, the equations may include a description of intracellular processes such as calcium handling by the SR. Over the past 50 years, more than 200 cardiac AP models have been developed [141] for animal and human species, including human SA node [55], human Purkinje fibre [231], [247], human atrial [42], [76], [229], [38], human ventricular [118], [251], [77], [177], [12], [243], [78], [244] CMs and human iPSC-CMs [185], [107], [104], [182]. The level of complexity ranges from phenomenological descriptions with a handful of state variables and model parameters to detailed mechanistic models based on experiments with dozens of state variables and hundreds of parameters, see the appendix for simple examples. Typically, the resulting ODE systems are stiff such that explicit methods with small step sizes or implicit methods are needed for their numerical integration.

Usually, the transmembrane currents flowing through the voltage gated ion channels are modelled [103] as

$$ I_{\text{ion}}(V,t) = N_{\text{ion}} \cdot O_{\text{ion}} (V,t) \cdot i_{ \text{ion}}(V) $$

(2)

with the current-voltage relationship $i_{\text{ion}}(V)$ of a single open channel and the total number $N_{\text{ion}}$ of channels of the specific type under consideration. In most cases, the probability $O_{\text{ion}}$ that the channel is open, is described by a monomial of voltage dependent gating variables, which are part of the state vector of the underlying cardiac AP model and take values between $[0,1]$, see the appendix for an example. An alternative is to derive $O_{\text{ion}}$ from the Kolmogorov differential equation associated with a Markov process that describes voltage dependent random fluctuations between different open and closed configurations of the channel [62].

Viewing the channel as an ohmic resistor, the function $i_{\text{ion}}(V)$ is commonly described by the linear relationship

$$ i_{\text{ion}}(V) = g_{\text{ion}} (V - V_{\text{ion}}). $$

(3)

Here, $g_{\text{ion}}$ denotes the single channel conductance, such that the maximum conductance of the macroscopic transmembrane current $I_{\text{ion}}$ is given by

$$ G_{\text{ion}} = N_{\text{ion}} g_{\text{ion}}. $$

(4)

At the Nernst potential

$$ V_{\text{ion}} = \frac{RT}{zF} \text{ln} \left ( \frac{[\text{ion}]_{e}}{[\text{ion}]_{i}} \right ) $$

with universal gas constant $R$, absolute temperature $T$, valence $z$ of the ion and Faraday constant $F$, the transmembrane voltage is balanced by the concentration difference of the ion across the membrane, resulting in a zero net flow of the ion through the channel. An alternative model of the current-voltage relationship is the Goldman-Hodgkin-Katz (GHK) equation

$$ i_{\text{ion}}(V) = p_{\text{ion}} \frac{z^{2} F^{2}}{RT} V \frac{[\text{ion}]_{i} - [\text{ion}]_{e} \text{exp} \left ( \frac{- z F V}{RT} \right ) }{1 - \text{exp} \left ( \frac{- z F V}{RT} \right ) }, $$

(5)

where $p_{\text{ion}}$ is the permeability of the membrane to the ion, normalized by the number of open channels. It can be derived by integrating the Nernst-Planck equation, which describes the ion flow in dependence of the concentration gradient and of the electric field through the membrane, under the assumption of a constant field. As in case of the ohmic model (3), the ionic current (5) is zero at the Nernst potential $V = V_{\text{ion}}$.

Typically, APs of spontaneously active cells are modelled by autonomous ODE systems that feature corresponding stable limit cycles. Models of quiescent CMs rather possess a stable steady state in accordance with the resting membrane potential (RMP) of the cell. In this case, the simulation of an AP requires to provide a suprathreshold current stimulation $I_{\text{stim}}$, see (1), that drives the system into an area of phase space from where it returns to the steady state only after a large detour. Figure 8 (left) shows APs computed with the Paci model [185] for a spontaneously beating hiPSC-CM, and Fig. 8 (right) illustrates the AP response of the O’Hara Rudy (ORd) model [177] of a human ventricular CM to a periodic stimulation with a pacing cycle length (CL) of 1000 ms. This pattern is referred to as a 1:1-rhythm as, after the vanishing of possible transients, each stimulus evokes one and the same AP.

Fig. 8

left: hiPSC-CM models feature limit cycle behaviour that corresponds with AP automaticity, simulation run with the Paci model [185]. right: Models of quiescent CMs require suprathreshold stimulation in order to produce APs. The figure shows a simulation of the ORd model [177] of the human ventricular mid-myocardial cell that corresponds with a 1:1 pattern in response to a pacing cycle length (CL) of 1000 ms

Common readout parameters of cardiac AP simulations include the maximum upstroke velocity $\text{(dV/dt)}_{\text{max}}$, the duration of the AP at $x \%$ of repolarization (APDx), or the AP amplitude (APA), see Fig. 9 (left). Likewise, characteristics of the calcium dynamics may be of interest. For models of quiescent CMs, another important feature is the APD restitution APD = r(DI), which describes the nonlinear dependency of APD (typically APD90) on the previous resting interval referred to as the diastolic interval (DI), see Fig. 9 (right). As discussed in the following, steepness of the restitution curve r is implicated in arrhythmogenesis both at the cellular and the tissue level.

Fig. 9

left: Common readout parameters, also referred to as metrics or features, from an AP simulation include the resting membrane potential (RMP), the maximum upstroke velocity $\text{(dV/dt)}_{\text{max}}$, AP duration at $30\%$ of repolarization (APD30), APD50, APD90, the diastolic interval (DI) and the AP amplitude (APA). right: Example (from [110]) of an APD restitution curve APD = r(DI) with CL = APD+DI, cobweb iterations are generated by $\text{APD}_{n+1} = \text{r}(\text{DI}_{n})$. If CL is chosen such that $\text{r}'(\text{DI})<1$, the intersection of r with the straight line CL-DI defines a stable fixed point that corresponds with a 1:1 pattern. Crossing $\text{r}'(\text{DI})=1$ by lowering CL, the stability of the fixed point is lost and APD alternans of the 2:2 type are generated via a period doubling bifurcation

4.2 Nonlinear Dynamics of AP Models and Cellular Arrhythmias

Cardiac AP models may possess a rich repertoire of nonlinear dynamics [110] that includes proarrhythmic behaviour such as APD alternans and early afterdepolarizations (EADs), see Fig. 10. These abnormal activities are associated with bifurcations [100] in dynamical systems, which are qualitative changes in the dynamics upon variation of model parameters [121]. In case of APD alternans, their onset is linked with a period doubling bifurcation of the iterated map

$$ APD_{n+1} = \phi (APD_{n}), $$

with $\phi (\text{APD}_{n}) = r(\text{CL}-\text{APD}_{n})$ and $\text{CL} = \text{DI}_{n}+\text{APD}_{n}$. For large values of the parameter CL, the fixed point $\text{ADP}_{*}$ of the map $\phi$ lies on the flat part of the restitution curve r and is stable. This situation corresponds with the 1:1 rhythm shown in Fig. 8 (right). A reduction of CL moves the fixed point $\text{APD}_{*}$ towards the steeper part of r until $r'(\text{DI}_{*}) > 1$ (which corresponds with $\phi '(\text{APD}_{*}) < -1$) and the stability of $\text{APD}_{*}$ is lost. Consequently, the period of the map $\phi$ is doubled from 1 to 2 and one obtains $\text{APD}_{n+2} = \text{APD}_{n}$ in correspondence with the 2:2 rhythm shown in Fig. 10 (left). As CL is further reduced, more complex m:n (m stimuli evoke n APs) or irregular patterns of APD alterations may occur [110].

Fig. 10

left: Simulation of APD alternans that correspond with a 2:2 pattern as every two supra-stimuli produce two different AP responses. right: Simulation of APs with EADs that correspond with a $1^{3}$ mixed mode oscillation (MMO) as 1 large amplitude oscillation alternates with 3 small scale oscillations

Cardiac APs with EADs as illustrated in Fig. 10 (right) are an example of mixed mode oscillations (MMOs) [49], which are characterized by an alternation between oscillations of large and small amplitudes. MMOs may result from the multiple time scales at which the ODE state variables operate, and a separation of the variables into slow and fast ones can be used for their analysis. Using a parsimonious AP model with one slow channel gating variable, the occurence of EADs has been attributed to a supercritical Hopf bifurcation [245], [204] in the fast subsystem, where the slow variable then serves as the bifurcation parameter, see Fig. 11 (left). As the slow variable increases and passes through the bifurcation point, the trajectory of the full system may momentary coil around the emerging stable limit cycles, which consequently causes the appearance of EADs. In [114], a delayed subcritical Hopf bifurcation was identified as an alternative EADs generating mechanism, where the twisting of the trajectory only sets in after the slow variable has passed the region of the limit cycles that now are unstable and oriented towards the opposite direction. Using the same AP model but rather working with two slow gating variables than just one, a folded node (FN) singularity was detected in the slow subsystem [116]. Then, MMO theory implies that, in vicinity of the folded node, the trajectories of the full system are organized by twisted slow manifolds. This defines a further dynamical mechanism of EADs [116] and also explains their dependency on the pacing frequency in case of external stimulation [264]. In addition, FNs may jointly occur with saddle-focus (SF) equilibria of the full AP system, see Fig. 11 (right). In the SF mechanisms behind EADs [114], [116], the trajectory first approaches the equilibrium along its stable manifold and then spirals away on the unstable manifold associated with a pair of complex conjugate eigenvectors with positive real part. Viewing the transmembrane voltage as the only fast variable, sufficient parametric conditions for the genesis of EADs have recently been derived in [31], while isolas of periodic orbits have been held responsible for the different EADs patterns in [11].

Fig. 11

left: Bifurcation diagram of the fast AP subsystem with the slow gating variable $x$ as continuation parameter. Stable limit cycles with growing diameter arise at a supercritical Hopf bifurcation and terminate at a saddle-homoclinic bifurcation. The solid and dashed black curves denote stable and unstable fixed points of the fast subsystem. The green and red curves show projections of two AP trajectories of the full model onto its phase space, the difference is due to respective parameter settings. The green curve passes by the basins of attraction of the fast subsystem and corresponds with a normal AP without EADs. The red curve is attracted towards the stable fixed points and limit cycles, and temporarily coils around them such that EADs are generated. Figure credit [114] (adapted under CC-BY license). right: Projection of an AP trajectory with EADs onto the model phase space. First, the EADs are organized by the twisted slow manifold due to a folded node (FN), then EADs are maintained by the spiraling along the unstable manifold of the saddle-focus equilibrium (SF) with tangential space $M_{u}$. $F^{+}$ denotes the fold line associated with FN, $e_{s}$ denotes the eigenvector corresponding with the negative eigenvalue of SF. Figure credit [116] (adapted under CC-BY license)

In an alternative approach to EADs analysis that does not exploit multiple time scales, the dynamical mechanisms of EADs generation and termination in periodically paced models of the human ventricular AP were shown to depend on the bifurcation structure of the non paced model counterparts [120], [119]. Of note, AP models may also feature EADs alternans with an alternating number of small scale oscillations [223], [264] or alternating EADs amplitudes [183], as also observable in experiments, see Fig. 6 (left). Furthermore, AP models may possess multiple stable EADs dynamics [275], [249], [111] that differ in the number of small scale oscillations, where only the initial conditions used for the simulation decide which type of dynamics manifests itself. Likewise, multistable AP dynamics with different m:n patterns are possible [237], [236]. In particular, multistability serves as an explanation for transient arrhythmogenesis [274] and intermittent occurence of cardiac arrhythmias [275], [111].

While APD alternans and EADs are considered to be proarrhythmic and sometimes are referred to as cellular arrhythmias, the simulations shown in Fig. 10 correspond with periodic, hence rhythmic AP model behaviour. However, deterministic AP models may also feature chaotic dynamics, i.e., aperiodic long-term behaviour that is highly sensitive to the initial conditions. For instance, Fig. 12 (left) displays an AP trace with an irregular sequence of APDs, obtained with the Beeler Reuter AP model under periodic pacing at high rates [201]. Chaotic AP dynamics is also possible in the presence of EADs, as illustrated in Fig. 12 (right) for the periodically paced Luo-Rudy AP model [245]. In [245], [223] the appearance of chaotic EADs is linked to the homoclinic bifurcation in the fast subsystem, at which the stable limit cycles of the fast variables terminate, see Fig. 11 (left). In [115], it was demonstrated that this homoclinic bifurcation is not a necessary condition for chaotic EADs dynamics and that the latter may also result from a cascade of period doubling bifurcations of limit cycles in the full system. This period-doubling route to chaos was detected both in periodically paced models and in an unpaced model [115] of a spontaneously active cell. Furthermore, bifurcation analysis showed that chaotic EADs dynamics may coexist [115] in a stable manner with periodic AP behaviour without EADs, where only the initial conditions decide which type of dynamics occurs. Finally, irregular EADs dynamics are also associated with random ion channel fluctuations [223], that are modelled by Markov processes or stochastic ODEs [57]. Random dynamical systems [8] may show behaviour that is not present in corresponding deterministic systems, see [230] for a recent analysis of noise-induced EADs.

Fig. 12

Examples of chaotic AP dynamics under high frequency pacing. left: Simulation of the Beeler Reuter AP model with a chaotic alteration of the AP duration, CL = 100 ms. right: Simulation of the Luo Rudy AP model with chaotic occurence of EADs, CL = 700 ms. Figure credit [115] (adapted under CC-BY license). Here, the chaotic state coexists with regular AP dynamics (not shown), obtained with identical model parameters and CL but different initial conditions

4.3 Modelling of AP Wave Propagation

The propagation of a wave of electrical excitation through some cardiac tissue preparation or the heart organ $\Omega$ (or some part of it such as the ventricular chambers) is commonly described by the bidomain model [250], [235]

$$\begin{array}{rlrlrlrl}\frac{\partial s}{\partial t}& & & =& & f(s,V,t)& & x\in \mathrm{\Omega },\\ \mathrm{\nabla }\cdot ({\sigma }_i\mathrm{\nabla }V)+\mathrm{\nabla }\cdot ({\sigma }_i\mathrm{\nabla }{u}_e)& & & =& & C_m\frac{\partial V}{\partial t}+\sum I_{ion}(s,V)& & x\in \mathrm{\Omega },\\ \mathrm{\nabla }\cdot ({\sigma }_i\mathrm{\nabla }V)+\mathrm{\nabla }\cdot (({\sigma }_i+{\sigma }_e)\mathrm{\nabla }{u}_e)& & & =& & 0& & x\in \mathrm{\Omega },\end{array}$$

(6)

with the transmembrane potential $V(x,t)$, the extracellular potential $u_{e}(x,t)$ and a vector $s(x,t)$ of cellular state variables. The nonlinear ODE system in (6) describes gating processes and electrochemical reactions in the cardiac cells, while the two partial differential equations (PDEs) describe the diffusion of the electrical signal. $\sigma _{i}$ and $\sigma _{e}$ denote the intra- and extracellular conductivities that are orthotropic due to the organization of the tissue into fibres and laminar sheets. Furthermore, cardiac tissue is heterogeneous such that model properties and components such as excitability, conductivity, ionic currents, channel conductances may be nonuniform and varying across $\Omega$. Different action potential models (1) may be incorporated into (6) in order to model different cell types that occupy spatially distinct but coupled regions as, e.g., in case of the transmural heterogeneity of the ventricular wall. The bidomain model yields a macroscopic description of cardiac electrical activity, in which both the intracellular and the extracellular domain are continuous and occupy all of $\Omega$, so both $V$ and $u_{e}$ are defined on all of $\Omega$.

The numerical simulation of the bidomain model [235], [64], in particular on realistic three-dimensional heart geometries with detailed cellular models, requires advanced discretization and solution techniques in order to properly capture the spatial and temporal variations of the wavefront. Under the assumption that the intra- and extracellular conductivity tensors are aligned, i.e., $\sigma _{e} = \lambda \sigma _{i}$, the bidomain model reduces to the monodomain model

$$\frac{\partial s}{\partial t}=f(s,V,t)x\in \mathrm{\Omega },$$

(7)

$$\frac{\lambda }{1+\lambda }\mathrm{\nabla }\cdot ({\sigma }_i\mathrm{\nabla }V)=C_m\frac{\partial V}{\partial t}+\sum I_{ion}(s,V)x\in \mathrm{\Omega },$$

(8)

which can be solved at less costs but still requires high performance computing in realistic applications, see for instance [25] for a coupled monodomain solver with optimal memory usage. Figure 13 shows the simulation of a normal AP wave propagation across the left and right ventricular chambers of a human heart model, using the monodomain equations with anisotropic conductivity and different AP models for human ventricular endocardial, mid-myocardial and epicardial cells as well as Purkinje fibre cells in order to account for regional specificity [40]. Numerical comparisons [65], [198], [19] between the mono- and the bidomain model indicate differences in wave propagation, while a theoretical result of [41] states that the error between the mono- and the bidomain solutions is bounded by the error between the respective conductivity operators.

Fig. 13

Normal propagation of the excitation wave across the ventricular chambers of a human heart model [40], snapshots are taken between the beginning of the QRS complex and the end of the T wave, compare with Fig. 2 (left). The domain is discretized with approximately 7 million linear hexagonal finite elements of edge length 0.3 mm, which gives rise to 270 million internal variables. The Purkinje fibre network is modelled with approximately 40000 linear cable elements and 800000 internal variables. The time series plot shows the pseudo ECG corresponding to 5 heart bearts and calculated by means of (9) with an electrode located 2 cm away from the left ventricular wall. Figure credit [40] (adapted under CC-BY license)

A common postprocessing step both for the monodomain and the bidomain model is to calculate pseudo-ECG signals [69], [9], [91], [152] by a projection of the transmembrane voltage gradient according to

$$ ECG(t) = - \int _{\Omega }\nabla V(x,t) \cdot \nabla \left ( \frac{1}{\|x_{e} - x \|} \right ) \, dx, $$

(9)

where $x_{e}$ is the location of an electrode placed in some distance from $\Omega$. More realistic simulations of ECG signals recorded at the body surface can be obtained by coupling the bidomain equation with a model of the human torso [235], [152] that views the latter as a passive volume conductor. Alternatively, ECG signals can be derived via the boundary element method [199], that involves integrals of $u_{e}$ over the torso surface, or via the lead field method [197], that exploits the Green’s function of the bidomain diffusion operator.

4.4 Modelling of Cardiac Arrhythmias at the PDE Level

The mono- and bidomain equations are reaction-diffusion models of excitable cardiac tissue and can show a variety of spatiotemporal patterns, ranging from normally propagating waves over spiral waves to spatiotemporal chaos. Patterns associated with cardiac arrhythmias result from a collision of travelling waves with mutual annihilation due to refractoriness, and they even may arise in isotropic and homogeneous tissue, where all model properties are uniform and do not change across $\Omega$.

APD Alternans and EADs in Cardiac Cables

A cardiac cable occupying $\Omega = [0,L]$ is the simplest example of a spatially extended system, in which APD oscillations may vary throughout the domain. APD alternans are called concordant if the whole domain shows the same long-short pattern in APD. In contrast, discordant APD alternans are characterized by a spatial gradient of APD from short to long that, on the next beat, is inverted and gets long to short. Depending on the length $L$ of the cable and the pacing period CL [58], discordant APD alternans may lead to a conduction block, where only every second beat propagates along the entire cable, see Fig. 14 (left).

Fig. 14

left: A homogeneous cardiac cable of length L = 20 cm is paced at $x=0$ with CL=230 ms. Discordant APD alternans form and lead to a conduction block, as only every second pulse is propagated to the other end of the cable. Example taken from [58]. right: A heterogeneous cable of length L = 9 cm where a domain $[0,0.5]$ of cells with repolarization abnormalities is coupled to a domain $[0.5,9]$ of spontaneously active cells without abnormalities. EADs form at the left end and trigger premature APs that propagate towards the other end. Example previously unpublished

If EADs are large enough, they may form premature triggers that propagate as premature ventricular complexes (PVC) into neighbouring areas. Figure 14 (right) illustrates PVC propagation through a heterogeneous cardiac cable, in which cells with repolarization abnormalities on the left are coupled to spontaneously active cells with normal repolarization on the right. For an example with a homogeneous cable, where multiple PVCs form due to partial regional synchronization of chaotic EADs and maintain uncoherent electrical activity even after termination of external stimulation, see [224].

Spiral Wave Break up in Cardiac Tissue

Figure 15 displays the break up of a spiral wave and the onset of chaotic spatiotemporal reentry activity in a 2D simulation of homogeneous cardiac tissue with $\Omega = [0,L] \times [0,L]$. In this example [228], the spiral is first produced by partially blocking a plane wave that is traveling from left to right. The broken wave gains velocities in both longitudinal and transverse directions such that it curls and builds a spiral with multiple turns. As the spiral core meanders through space, a part of it dissociates and collides with the remaining part of the spiral. This leads to spontaneous break up and initiation of multiple smaller spirals. These spirals drift through space and interact in an irregular and self-sustained manner, that causes chaotic high frequency excitation of the tissue, until they finally disappear and the tissue comes to a rest.

Fig. 15

left: Formation of a spiral wave in response to a conduction block, and subsequent break up. Multiple smaller waves arise and interact in a self-sustained and chaotic manner for a long period of time before they finally die out. Example taken from [228]. The time series plot shows the corresponding pseudo ECG calculated by means of (9) with an electrode located above the center of the domain. Figure previously unpublished

Different mechanisms of spiral wave break up have been identified in simulation experiments [58], [228] with the monodomain model, that may be further promoted by tissue heterogeneities, anisotropy or three spatial dimensions. They include steepness of APD restitution, and dynamically caused dispersions in APD, refractoriness and conduction velocities, but also localized heterogeneities with impaired conduction or spatial gradients in material properties. The liability of APD alternans in 2D for spiral wave break up in both homogeneous and heterogeneous tissue has been addressed in several computational studies [203], see [251], [266] for examples with human anisotropic but homogeneous tissue and a focus on sodium channel function. The role of the latter in APD alternans-mediated chaotic reentry activities is also highlighted in [79] using an isotropic and homogeneous 2D tissue with a parsimonious rabbit AP model.

The particular role of EADs in the generation of cardiac arrhythmias has been investigated in a series of computational studies using the monodomain model in a variety of settings that include animal or human cell models, 2D or 3D geometries, isotropic or anisotropic conductivities $\sigma _{i}$, and homogeneous or heterogeneous tissue, see [260] for a recent review. Typically, a part or all of $\Omega$ is made prone to EADs by reducing the repolarization reserve of the underlying cell model, before the response of the affected tissue to an excitation wave is studied. The observed arrhythmic behaviour includes meandering of reentrant waves due to EADs formation at the spiral core [88], focal activity due to local synchronization of chaotic EADs [224], and spiral fibrillation with multiple small rotations due to wave breaks caused by EADs [259], [287], [260]. The occurence of EADs induced PVCs depends on a variety of factors that include the size of the EADs cell cluster, the coupling strength between the EADs cells but also the pacing cycle length CL [286]. In heterogeneous and anisotropic models of the human ventricular chambers, EADs may promote reentrant activity [52] and induce TdP arrhythmias via focal and/or reentrant activity [258] in dependence on the difference in repolarisation reserve between the heterogeneity and the surrounding tissue.

The majority of the computational electrophysiology studies of cardiac arrhythmias focuses on ventricular cells and chambers, as TdP and VF can lead to sudden cardiac death, see Fig. 1 (right). However, given the health burden of atrial fibrillation (AF) as the most common sustained type of arrhythmia [32], also the interest in modelling of atrial cells [83] and atrial chambers [37], [256] is continuously growing.

5 Simulation Based Cardiac Drug Safety Testing

Computational models of cardiac electrophysiology find applications [172] in cardiovascular disease research or medical device and treatment design, but they are also extensively used in the context of cardiac drug safety testing. The most common way of describing the impact of a drug on cardiac ion channels [21], [158] is the pore-block model. There, the maximum ion conductance (4) (or the maximum ion permeability) in the channel current formulation (2) is simply reduced by a factor

$$ \left ( 1+\left [ \frac{[D]}{\text{IC50}_{\text{ion}}} \right ]^{h_{ \text{ion}}} \right )^{-1} $$

(10)

in dependence on the applied drug and its concentration $[D]$. Here, $h_{\text{ion}}$ denotes the Hill coefficient, where $h_{\text{ion}}=1$ assumes that one drug molecule is sufficient to block one ion channel. Typical non-integer values of $h_{\text{ion}}$ lie between 0.4 and 1.6. For instance, the simulation of a compound that blocks both the hERG and the calcium channel requires knowledge of values of $\text{IC50}_{\text{Kr}}$, $h_{\text{Kr}}$ and $\text{IC50}_{\text{CaL}}$ $h_{\text{CaL}}$, to be derived from experiments or read out from data bases, in order to calculate the corresponding scaling factors for $G_{\text{Kr}}$ and $G_{\text{CaL}}$ of that drug. More elaborate descriptions of drug-channel interaction address the kinetics of drug binding to the channel [50], [158], [129], [137], [238] and are based on Markov models of the channel that involve drug bound states.

At the single cell level, one usually studies changes in AP morphology, such as APD90 prolongation or the occurence of EADs and repolarization failure, and in calcium transients, that arise in response to parameter modifications in the AP model (1) motivated by the presence of a drug, see [29], [175], [181] for examples with human ventricular and hiPSC cardiomyocyte models. In silico drug trials are also run on cardiac cell models adapted to mimic disease states such as LQTS and ischemia [149], [30] or on populations of cell models [22], [183], [160], [182] see Fig. 16 (left), that are obtained by random sampling of the parameter space, in order to take intra- and intersubject variability into account. Based on the fraction of models for which abnormalities occur under the influence of reference compounds with known TdP risk, risk scores can be derived from multiple simulations that possess higher accuracy than animal models in predicting proarrhythmicity [191]. Recently, blinded studies on AP model populations also suggested that the minimum set of ion channels necessary for reliable TdP risk predictions is given by Nav1.5(peak), Cav1.2 and hERG [283].

Fig. 16

left: In silico drug trial on a population of human ventricular AP models. First, the population is generated by random sampling of ionic conductances in a $[0-200] \%$ range of the baseline ORd model (solid line) values and by only keeping those AP profiles that are in correspondance with a normal repolarization. Next, the drug (dofetilide) is applied in terms of the pore block model and the number of resulting AP phenotypes with repolarization abnormalities (RA) is counted. The fraction of models with RA over the total number of models in the population defines a TdP risk score between 0 and 1. Figure credit [191] (reproduced under CC-BY license). right: The CiPA TdP risk metric qNet as a function of drug concentration for the 12 CiPA training compounds of Table 1. qNet is defined as the area under the curve of the net current during a simulated AP of the ORd human ventricular model and consistently classifies each of the compounds into the correct TdP risk category across a wide range of concentrations. Figure credit [262] (reproduced under CC-BY license)

Single cell AP models (1) are also integrated into statistical learning algorithms [82] such as linear determinant analysis (LDA), logistic regression, support vector machine (SVM) or principal component analysis (PCA) in order to perform risk classification of drugs. One observation made is that risk classifiers based on indirect drug features derived from model simulations such as APD90 [153], APD50 and diastolic $\text{Ca}^{2+}$ [122], [112] or EADs tendency [149], [189], [160], [94] may outperform those solely based on direct features [109], [200], [127], [16] obtained from experiments, see also the appendix. Within the CiPA initiative, several derived features have been tested by means of an optimized version [51] of the ORd human ventricular cell model, now incorporating $I_{Kr}$-drug binding dynamics [137], and the 12 training compounds listed in Table 1. The risk metric found to best stratisfy the drugs into the three TdP risk levels across various conditions is qNet, see Fig. 16 (right), which is defined as the charge passed by the net ionic current from the beginning to the end of the AP beat. Furthermore, qNet was shown to correlate with the system’s robustness against EADs, in the sense that low values of qNet make the cell more prone to EADs induction by an $I_{\text{Kr}}$ reduction in addition to the drug block effect.

Likewise, PDE models of human cardiac electrophysiology find application in in-silico drug trials on 1D cardiac cables [246], 2D ventricular wedge preparations [113] and 3D models of the ventricular chambers [159], [280], [40], often also embedded into torso models [269], [281], [282], [178], [92] for computing impacts on ECG signals. Both the mono- [159], [269], [280], [40] and the bidomain [178], [282], [113] model are in use and may quickly pose high computational demands if anatomical geometries, realistic conductivity tensors, and detailed models of drug action and of heterogeneous cellular electrophysiology are applied. Typically, for known safe (e.g. verapamil) and known torsadogenic (e.g. dofetilide) drugs the underlying ion channel models first are manipulated based on corresponding pharmacological profiles. Then, excitation wave propagations and ECG signals are simulated under various drug doses, before concentration thresholds, at which arrhythmic behaviour possibly occurs, are reported, see Fig. 17 for an example. Furthermore, these simulation studies point to mechanisms that possibly are responsible for the different proarrhythmicity of the drugs, despite them having similar channel blocking potencies or impact on QT interval prolongation. They include accumulation and recovery times during frequency dependent sodium channel block [159], flattening of APD restitution [269], [280] or EADs translating into pronounced U waves in the ECG [40]. Given the computational burden of drug parameter space screening with 3D heart models, Gaussian process machine learning is suggested in [39] for building a surrogate model of the QT interval and in [219] for more efficiently constructing the classification boundary between safe and proarrhythmic domains.

Fig. 17

In silico drug trial on a model of the human ventricular chambers. The plot shows the evolution of the transmembrane voltage in response to dofetilide at concentration 30x in terms of the pore block model. Following a significantly prolonged QT interval, several spiral waves form and excite the heart in a chaotic manner, compare with Fig. 13 and Fig. 15. Figure credit [40] (adapted under CC-BY license)

Of particular interest for preclinical cardiac drug safety testing is modelling and simulation of monolayers of hiPSC-CMs, as used in high throughput in vitro recording tools such as MEA or VSD, compare with Fig. 7. These monolayers form a mixture of hiPSC-CMs with atrial-like, ventricular-like and nodal like phenotypes that may not be partitioned into clearly distinct regions and rather are spatially disorganized. In order to take a spectrum of AP morphologies throughout a well-mixed tissue into account, homogenisation is proposed in [20] to correspondingly extend the mono- and bidomain equations, while a probabilistic description of cell heterogeneity using correlation matrices is proposed in [241], [210]. There, the bidomain model furthermore is coupled to an electrode model in order to simulate field potential recordings from MEAs. These FP simulations are used to enrich experimental MEA measurements in order to improve the performance of SVM and LDA risk classifiers of drugs based on a composition of FP features such as FP duration or repolarization amplitude. An alternative and computationally much cheaper tool for the simulation of FPs from MEAs is discussed in [239], where the transmembrane voltage of the single cell ODE model (1) is postprocessed by means of SPICE based software through an electronic circuit that resembles the junction between the cell membrane and the MEA electrode, and that acts as a high pass filter with the FP signal as the output.

6 Challenges and Outlook

Computational cardiology has continuously grown over the past 50 years and has turned into a prominent showcase of systems biology [43]. One major area of application is cardiac drug safety testing [217], [68], [253], which recently gained novel momentum in course of the emergence of hiPSC-CM technology and the initiatives built upon it for an overhaul of the current safety paradigm [226], [61], [98]. While proarrhythmic risk classification of drugs based on mechanistic models of human cardiac electrophysiology [191], [262], [242], [219] gains more and more recognition also within pharmaceutical industries and regulatory agencies, still there are a couple of hurdles on the way to its full integration into critical decision making. Some of these challenges are listed in the following.

Validation, Verification and Uncertainty Quantification

One important aspect in the use of computational models for safety-critical applications is the establishment of their trustworthiness through validation, verification and uncertainty quantification [3]. Here, validation refers to the assessment of the degree to which the mathematical model properly reflects the reality of interest. For instance, a comparison of EADs [87] occuring in AP models with those observed in experiments revealed rather large discrepancies in inter-EAD intervals and attributed them to shortcomings of calcium current formulations. As another example, inconsistencies between AP models and experiments were found in [255] with respect to discordant alternans and linked with limitations in the description of voltage-calcium coupling. Often, cardiac AP models (1) are overparametrized and degenerate [53] such that completely different combinations of model parameters may lead to the same or to a nearly identical AP morphology [222], [96]. Furthermore, most AP models have been built by manually combining submodels based on experimental data, that were obtained from a range of cells, often even across species, and under a variety of environmental conditions. These problems can be addressed by collecting more informative and specific data sets, e.g., from single cardiomyocytes in a series of designed experiments [80] or from clinical trials [146], [112] in order to better constrain the model. This requires to differ between structurally and practically identifiable model parameters [46], [268] as well as to quantify their uncertainties [194], [132]. One approach is to use Bayesian methods in order to infer probability distributions of the AP model parameters instead of individual numbers [47], [95]. Motivated by biological variability [272], [169], an alternative is to build populations of models by statistical sampling of model parameters, where the population is either calibrated to ranges of data values [22] or to data distributions [126]. Uncertainties in both model inputs and outputs affect cardiac electrophysiology modelling across all scales [195], [155]. In particular at the ion channel level [35], computational strategies are under development that define experimental protocols for novel high-throughput recording systems [13], [130], in order to improve model credibility. For a mathematical model, that not only accounts for variability in ion current data but also in the underlying voltage-clamp experiment, see [131].

Verification deals with the correct implementation of the mathematical model and its accurate numerical solution [3]. In the context of computational cardiology, this in particular affects reaction diffusion models based on the mono- and the bidomain equations, see (8) and (6). Their spatial discretization by means of finite elements leads to a nonlinear system of algebraic-differential equations that needs to be integrated with respect to time either in a coupled or a decoupled manner [235], [64]. Decoupled methods, such as operator splitting [235] or semi-implicit time discretization [54], separate the diffusion term from the reaction term or the PDEs from the ODEs, treat them by different numerical schemes in order to increase numerical efficiency, and only have to deal with linear systems at each time step. Coupled methods yield more accurate results than decoupled ones and their time stepping is not subject to stability constraints [64]. However, coupled methods are much more expensive as they require to solve a nonlinear system at each time step. Consequently, they have been applied to PDE models only in 1D [279] or only in combination with phenomenological AP models [163], but so far are prohibitive in higher spatial dimension for physiological AP models with dozens of state variables. Recently, a coupled monodomain solver with optimal memory usage for 2D and 3D applications with mechanistic AP models has been introduced in [25] that exploits the sparsity of local matrices in the assembly of global FEM matrices. Currently, this methodology is extended to the solution of the bidomain equations with space-time adaptive discretization [26], see Fig. 18. While the progression of high performance computing (HPC) paths the way towards realistic simulations of human cardiac electrophysiology in real time at least with operator splitting methods for the monodomain model [171], [213], reduced order modelling of the mono- and bidomain equations [18], [186] and the reduction of physiological AP models [142] form further methodological strategies for bringing down complexity and computational costs. Another approach for the acceleration of mono- and bidomain solvers is the use of graphical processor units (GPU) [165], [273], [150], [97], which for benchmarks with 3D geometries enable speed ups by a factor of 10 to 200 in comparison with CPU implementations.

Fig. 18

Parallel space-time adaptive solution of the bidomain equations in 3D coupled with the ten Tusscher human AP model. top row: Space-time evolution of the transmembrane potential with initiation of cardiac reentry. mid row: Corresponding evolution of the adaptive FEM mesh. The maximum number of finite elements is about 7 million, giving rise to about 27 million internal variables. bottom row: Corresponding distribution of the mesh on 48 cores via dynamic load balancing. Figure credit [26] (previously unpublished)

Simulation Based Proarrhythmic Risk Assessment

As outlined above, several methods for the computational assessment of a drug’s proarrhythmic risk have been suggested that involve mechanistic ODE or PDE models of cardiac electrophysiology. Given the computational burden and the higher degree of uncertainty of PDE models [242], [219], ODE models currently dominate the literature on simulation based cardiac drug safety testing [153], [122], [191]. In a proof of concept study [4], a single cell model was even shown to be equally accurate as a 3D ventricular chamber model in predicting abnormal electrical activity. In general, the performances of the proposed risk metrics and classification schemes are difficult to compare [271], also because nonuniform sets of drugs are considered and because individual compounds are, sometimes even inconsistently, placed into 2 up to 5 different risk categories (from low to high TdP risk) for the training of the respective classifiers. In the CiPA study [51], ternary classification based on proportional odds logistic regression and leave-one-out validation was chosen for the training compounds of Table 1. This determined qNet, see Fig. 16 (right), as the best performing metric out of 13 candidates derived from an optimized version of the ORd model. However, when exposed to uncertainties in the model parameters that describe drug binding and ionic current block, qNet may loose its ability to correctly separate drugs by their TdP risk [27], see Fig. 19 (left). In a study with a phenomenological AP model [116], qNet was shown to be a monotonically decreasing function of the normal distance from the region of EADs occurence in parameter space, see Fig. 19 (right), such that low risk drugs paradoxically are located closer to the critical domain than their high risk counterparts. Also, a global sensitivity analysis [188] of the CiPA ORd model revealed significant differences between the parameter sets that affect qNet and EADs generation, which further challenges the correlation of qNet with the robustness against EADs as postulated in [51]. Recently, a torsade mectric score [140], defined as an averaged qNet value, showed high accuracy in TdP risk prediction but still misclassified some of the validation compounds of Table 1. While further ODE simulation based risk markers, such as drug induced shortening of the electromechanical window [192] or concentration dependent repolarization reserve current [67], are currently under investigation, also a debate has been kicked off [189], [156], [157], [81] as to whether proarrhythmic risk classifiers built on features derived from cellular AP models actually provide increased prediction accuracies, in comparison with statistical learning approaches based on direct features such as ion channel blocking potencies [109], [156] or hiPSC-CM responses to drugs [16], [108], [193].

Fig. 19

left: (Output) uncertainty in the CiPA TdP risk metric for the 12 CiPA training compounds of Table 1, compare with Fig. 16 (right). The shaded areas denote the $95\%$ confidence intervals of qNet and indicate that the ability to correctly separate the drugs may get lost in the presence of (input) uncertainties in model parameters of drug-channel interactions. Figure credit [27] (reproduced under CC-BY license). right: Normal distance to EADs as a proarrhythmic risk metric for a simplified AP model. The red shaded area denotes the parameter region in which EADs occur, its boundaries are defined by bifurcations. Distance from the EADs area may serve as a safety parameter. However, approaching the area along a normal direction, the value of qNet increases. This is in contrast with Fig. 16 (right), where high values of qNet are associated with low proarrhytmic risk. Figure credit [116] (adapted under CC-BY license)

Modelling of hiPSC-CMs

A vital component of the CiPA initiative [262], [138] is to integrate drug effects on multiple ion channels into (an optimized version of) the ORd model of the human adult ventricular cardiomyocyte and to contrast model predictions with drug effects on human iPSC-CMs, see Fig. 5 (left). However, hiPSC-CMs differ from human adult CMs in morphology, electrophysiology and calcium handling, a situation often referred to as hiPSC-CM immaturity [214], [101], see Fig. 20. Furthermore, hiPSC-CMs possess a quite diverse spectrum of phenotypes, typically characterized by spontaneous beating, and only roughly categorizable into nodal-, atrial- and ventricular-like subtypes [99]. Hence, distinct computational models of hiPSC-CMs have been suggested as a tool to gain mechanistic insight into their immature phenotype [107], to interpret experimental findings in specific hiPSC-CM cell lines [133] and to translate from hiPSC-CMs to human adult CMs [73], [252], [93]. The series of papers [180], [181], [184], [185], [182] features corrective adaptations of a single cell ODE model (1) of a ventricular-like hiPSC-CM with membrane clock driven automaticity. In [107], the model of [181] is complemented by a PDE diffusion equation for the calcium release from the sacroplasmatic reticulum, which accounts for the calcium clock as the second mechanism behind the spontaneous activity of hiPSC-CMs. The dependency of the latter on the potassium current $I_{K1}$ formulation is studied in [56], while [104] suggests a novel hiPSC-CM model that is based on a simplified description of voltage-dependent gating variable rate constants and incorporates experimental hiPSC-CM data from multiple laboratories. Still, all of the available hiPSC-CM models, in parts and due to lack of data, contain ionic current and calcium handling formulations from previously published models of human adult or even animal CM models. Furthermore, they so far only enable oversimplified mono- or bidomain simulations of experiments with hiPSC-CM monolayers [209], [20] or 3D engineered tissue, compare with Figs. 6 and 7, as statistical data about the phenotypic mix in such multi-cellular systems but also corresponding cell models, e.g. for nodal-like hiPSC-CMs, are missing. As a consequence, simulated drug effects on hiPSC-CMs may deviate from in vitro observations, see [117] for a study on the relation between beating frequency and duration of repolarization, such that more effort is needed to increase data availability and to refine hiPSC-CM models. With respect to an improved characterization of individual ionic currents, the dynamic clamp technique [14], [261] allows to inject computationally simulated currents into isolated hiPSC-CMs in real time and to test model-based hypotheses on ion channels in a closed loop setting. Further applications of dynamic clamp include the improvement of hiPSC-CM immaturity [74] and the investigation of drug induced EADs [179].

Fig. 20

Comparison of hiPSC-CMs and human adult CMs. hiPSC-CMs are sometimes considered as moving targets as their characteristics strongly depend on time in culture. Late hiPSC-CMs differ from human adult CMs, e.g., in shape, size and ion channel expressions. Figure credit [214] (reproduced with permission)

Cardiac Electro-Mechanic Coupling

In this paper, we surveyed computational methods in cardiac electrophysiology and their relevance for preclinical drug cardiotoxicity testing with hiPSC-CMs. In doing so, we followed the paradigm according to which cardiac arrhythmias result from disturbances of ion channel function and excitation wave propagation that cause desynchronized contractions of the heart. However, the mechanical forces produced by excitation can themselves impair ion channel function and calcium handling [106]. Mathematical models of excitation-contraction coupling (ECC) and mechano-electric feedback (MEF) exist both at the cellular and the tissue/organ level. In the first case, electrophysiological AP models are coupled to models of the myofilament protein dynamics [89], [170], [212] that describe tension and changes of sarcomere length (SL) in dependence on the intracellular calcium concentration. In turn, the calcium dynamics are affected, e.g., by calcium-troponin binding or stretch-activated ion channels. For human species, such cellular models are used to study the modulation of AP alternans [285] and EADs [284] formation by contractile mechanisms, but they also are integrated into electro-mechanic models of the human heart [123], [10], [124], see Fig. 21. In general, cardiac mechanics has a passive component, corresponding with the tissue response to external forces, and an active component, referring to the force generation for the muscle contraction. While the active stress approach [164], [173] considers the first-Piola-Kirchhoff tensor as the sum of the passive and the active component, the active strain approach [28], [6] is based on a factorization of the deformation gradient tensor into the product of a passive and an active part. The resulting equations of continuum mechanics [33], with the muscle deformation and the stress tensor as unknowns, are coupled to the mono- (8) or bidomain (6) equations and are used, e.g., to study the inotropic effects of ion channel mutations [168] or the impact of contraction on spiral wave activity in cardiac tissue [164], [28]. For a recent simulation study of simultaneous drug effects on contractility and electrophysiology of the 3D human left ventricular chamber, see [147]. As also multicellular hiPSC-CM modelling progresses, these computational tools will also be used in ECC and MEF studies of force-generating human myocardium models [240], [135], [136] that are bioengineered with stem cell derived CMs, see Fig. 22 for examples. In particular, such models of the human myocardium allow to screen compounds not only for side effects on electrophysiology but also on contractility [127], [220], [233], e.g., drug induced sarcomere shortening, and on structure [276], e.g., drug induced changes in morphology.

Fig. 21

Simulation of the excitation-contraction coupling using a weakly coupled electromechanic model of the human heart, where the solution of the monodomain model is used as an input to the equations of continuum mechanics with active stress for describing muscle contraction. Model by courtesy of the Living Heart Project https://www.3ds.com/products-services/simulia/solutions/life-sciences/the-living- heart-project/

Fig. 22

left: A ring shaped human myocardium, bioengineered with about $1.5{\cdot }10^{6}$ stem cell derived CMs, and attached to two force sensing posts. Figure credit [240] (reproduced with permission). right: An electro-mechanically functional miniature human ventricle-like heart chamber, bioengineered with about $10^{7}$ hiPSC-CMs. Figure credit [136] (reproduced with permission)

Further topics, that go beyond the CiPA paradigm, as illustrated in Fig. 5 (left), include drug-drug interactions [151], drug induced disruption of hERG trafficking [174], or hidden drug cardiotoxicity [59], where drugs only negatively affect the diseased heart. Moreover, the prediction of a drug’s proarrhythmic risk might be improved by taking clinical exposures and metabolites under stressed scenarios into account [134]. In tackling these problems, computational models of drug adverse effects on cardiac function can be of great value as they allow to run through thousands of complex scenarios, which is not possible in in vitro or in vivo experiments. Underlying mechanistic models of the human heart are continuously refined and extended, e.g., by consideration of blood flow [234], [10], [207], [206] or of physiology based pharmacokinetics [254]. However, with big data also interfusing the field of drug toxicity testing, these models face the competition with machine learning algorithms [148], [145], [72], [23] and related bioinformatics approaches [208], [144], [162], [128] such as quantitative structure activity relationship (QSAR) [66], [71]. Ultimately, all these in-silico disciplines need to build even stronger bonds with human-relevant testing platforms [187] in order to actually achieve replacement, reduction and refinement (3Rs) of animal testing [176].

Appendix

Foundations of Cardiac AP Modelling

A prominent example of a phenomenological description of excitability is the FitzHugh-Nagumo model [63]. If adapted and reparametrized for cardiac AP simulation, the model reads as

$$\begin{aligned} \frac{dV}{dt} =& \frac{c_{1}}{v_{\text{amp}}^{2}} (V-v_{\text{rest}})(V-v_{ \text{th}})(v_{\text{peak}}-V) \\ & - \frac{c_{2}}{v_{\text{amp}}} (V-v_{\text{rest}}) W + I_{ \text{stim}} \\ \frac{dW}{dt} =& b (V-v_{\text{rest}} - c_{3} W), \end{aligned}$$

where $V$ represents the transmembrane voltage and $W$ is referred to as the recovery variable. Furthermore, $v_{\text{rest}}$, $v_{\text{th}}$ and $v_{\text{peak}}$ denote the resting potential, the threshold potential and the peak potential with $v_{\text{amp}} = v_{\text{peak}} - v_{\text{th}}$, see Table 2 for parameter values.

Table 2 Default parameter values of the reparametrized FitzHugh-Nagumo model [235]

Mechanistic AP models are based on viewing the membrane as an electric circuit with capacitive and conductive branches in parallel, see Fig. 23 for a schematic representation of the Noble model [45], which marks the beginning of cardiac AP modelling. The model describes the Purkinje fibre action and pacemaker potentials and reads as

$$\begin{aligned} C_{m} \frac{dV}{dt} = & -(I_{Na} + I_{K} + I_{L}) \\ \frac{dm}{dt} =& \alpha _{m}(V) (1-m) - \beta _{m}(V) m \\ \frac{dh}{dt} =& \alpha _{h}(V) (1-h) - \beta _{m}(V) h \\ \frac{dn}{dt} =& \alpha _{n}(V) (1-n) - \beta _{n}(V) n \end{aligned}$$

with the ohmic current formulations

$$\begin{aligned} I_{\text{Na}} = & (g_{Na} m^{3}h + 0.14)(V-V_{Na}) \\ I_{\text{K}} = & \left (1.2 e^{\frac{-V-90}{50}} + 0.015 e^{ \frac{V+90}{60}} + g_{K_{2}} n^{4} \right ) (V-V_{K}) \\ I_{\text{leak}} = & g_{L}(V-V_{L}) \end{aligned}$$

and $V_{Na}=40$ mV, $V_{K} = -100$ mV, $V_{L} = -60$ mV, $C_{m} = 12$ $\mu \text{F}/\text{cm}^{2}$, $g_{Na} = 400$ $\text{mS}/\text{cm}^{2}$, $g_{K_{2}} = 1.2$ $\text{mS}/\text{cm}^{2}$, $g_{L} = 0.075$ $\text{mS}/\text{cm}^{2}$. Here, $m$ is the sodium activation variable, $h$ is the sodium inactivation variable and $n$ is the potassium activation variable, and the voltage dependent rate functions are given by

$$\begin{array}{rlrl}{\alpha }_n(V)& =\frac{0.0001(-V-50)}{e^{\frac{-V-50}{10}}-1},& {\beta }_n(V)& =0.002e^{\frac{-V-90}{80}},\\ {\alpha }_m(V)& =\frac{0.1(-V-48)}{e^{\frac{-V-48}{15}}-1},& {\beta }_m(V)& =\frac{0.12(V+8)}{e^{\frac{V+8}{5}}-1},\\ {\alpha }_h(V)& =0.17e^{\frac{-V-90}{20}},& {\beta }_h(V)& =\frac{1}{e^{\frac{-V-42}{10}}+1}\end{array}$$

in units of $1/\text{ms}$. One interpretation [103] of these equations is that, e.g., the sodium channel consists of three gates of type $m$ and one gate of type $h$, each of which can be either closed or open. Under the assumption that the gates operate independently, the fraction of open sodium channels is given by $m^{3}h$. Then, $\alpha _{m}(V)$ and $\alpha _{h}(V)$ denote the transition rates from the closed to the open state of the respective gates, while $\beta _{m}(V)$ and $\beta _{h}(V)$ are associated with the transitions from open to closed. The equations for the gating variables $m$, $h$, $n$, which can be interpreted as open propabilities of subunits that constitute the channel, are often rewritten in the form

$$ \frac{dx}{dt} = \frac{x_{\infty }(V)-x}{\tau _{x}(V)} $$

with the steady state function $x_{\infty }(V)$ and the time constant function $\tau _{x}(V)$ given by

$$ x_{\infty }(V) = \frac{\alpha _{x}(V)}{\alpha _{x}(V)+\beta _{x}(V)}, \quad \tau _{x}(V) = \frac{1}{\alpha _{x}(V)+\beta _{x}(V)} $$

for $x \in \left \{ n,m,h \right \}$.

Fig. 23

The cellular membrane as a capacitor. left: Schematic representation of the Noble model [45] for the Purkinje fibre action and pacemaker potentials with potassium, sodium and leakage currents. Figure taken from the model repository [141] of the physiome project [90]. right: Equivalent electric circuit with capacitive and conductive branches. The capacitive branch corresponds to the insulating lipid bilayer, and the conductive branches stand for the ion channels with the Nernst potentials $V_{Na}$, $V_{K}$ and $V_{L}$ as electromotive forces. Figure generated with Matlab Simulink

The Noble model is a modification of the Hodgkin-Huxley (HH) equations [85] for the description of excitation in neural cells. Today, the HH formalism is present in nearly all of the modern mechanistic AP models and used for the description of most of the cardiac ionic currents, either in combination with the linear current-voltage relation (3) or the GHK equation (5). The HH description can be generalized to Markov models [62] that represent the channel by a discrete and finite set of open and closed states, which allows to also take relationships between gating processes into account. A Markov model is not a unique description of the ion channel function and different model structures that fit the same experimental data can have different interpretations with respect to underlying mechanisms [24]. Furthermore, Markov states are abstract objects that do not necessarily represent structural conformations of the ion channel protein [166]. Still, such phenomenological Markov models are the current method of choice for describing kinetic drug-channel interaction [50], [158], [129], [137] after a model extension by drug bound states.

In addition to detailed descriptions of a number of transmembrane currents, state of the art cardiac AP models (1) also feature ion pumps and exchangers. Furthermore, calcium fluxes into and out of the sacroplasmatic reticulum (SR) as well as calcium buffering are taken into account, which makes, e.g., the SR calcium concentration $[Ca]_{SR}$ and the myoplasm calcium concentration $[Ca]_{i}$ part of the ODE state vector in (1). Figure 24 illustrates the ORd model of the human ventricular CM and the Paci model of ventricular-like hiPSC-CMs and indicates the complexity of cardiac AP modelling.

Fig. 24

Schematic views of AP models with transmembrane currents, ion pumps and exchangers as well as intracellular calcium handling. left: ORd model of the undiseased human cardiac ventricular AP. Figure credit [177] (reproduced under CC-BY license) right: Paci model for the AP of a ventricular-like hiPSC-CM. Figure credit [180] (reproduced with permission)

Examples of Risk Classification Based on Experimental Features

The ICH S7B guidelines [2] for preclinical cardiac safety testing consider a compound to be torsadogenic if its hERG IC50 value, see Fig. 25 (right), and its effective free therapeutic plasma concentration (EFTPC) satisfy

$$ \frac{\text{IC50}_{\text{hERG}}}{\text{EFTPC}} < 30. $$

Figure 25 (right) illustrates that this criterion is too conservative and not specific, as not only the CiPA high TdP risk training compounds would accordingly be rated as torsadogenic, but also most of those with only medium or even low TdP risk. For instance, verapamil is a potent hERG blocker but still considered safe as it is in therapeutic use without any signs of torsadogenicity. One explanation for the high false-positive rate of the hERG assay is that negative hERG blocking effects on cardiac repolarization may be compensated by additional blockage of other channels such as the sodium NaV1.5 or the calcium CaV1.2 channel. In [109], several logistic regression models with $\text{IC50}_{\text{hERG}}$, $\text{EFTPC}$, $\text{IC50}_{\text{CaV1.2}}$ and $\text{IC50}_{\text{NaV1.5}}$ as direct features were trained based on data for 32 torsadogenic (TdP = 1) and 23 non-torsadogenic (TdP = 0) drugs. The multichannel effect (MICE) model, that was found to most significantly reduce false positives and false negatives in comparison with the hERG test, is given by

$$\text{logit}(\text{TdP}) = \beta _{0} + \beta _{1} \cdot \text{log} \left ( \frac{\text{IC50}_{\text{hERG}}}{\text{IC50}_{\text{CaV1.2}}} \right ), $$

Fig. 25

left: Dose-dependent inhibition of the hERG current $\text{I}_{\text{Kr}}$ by the drug dofetilide. At drug concentration $\text{IC50}$, the current is reduced by $50\%$. The solid line represents the sigmoidal function (10) with $\text{IC50}_{\text{hERG}} = 19.56$ nM and $h_{\text{hERG}} = 1.19$, obtained by curve fitting. Data by courtesy of Drug Discovery Sciences at Boehringer Ingelheim, Biberach. right: The ICH S7B criterion misclassifies most of the 12 CiPA training compounds from Table 1 as torsadogenic, because their IC50 value against the hERG channel is below the 30-fold of their EFTPC. Data from [211] and [109]. Figure previously unpublished

with $\beta _{0} = -0.75$ and $\beta _{1} = 3.04$. Note that it only involves the ratio of $\text{IC50}_{\text{hERG}}$ and $\text{IC50}_{\text{CaV1.2}}$. However, Fig. 26 (left) indicates that the model still lacks specificity in particular if applied to compounds not used for model development. In [225], several derivatives of the torsadogenic hERG blocker dofetilide were tested on hiPSC-CMs. While Dofe54 induced EADs, Dofe42 and Dofe44 had hardly any impact on AP morphology even at high dosage. Still, the MICE model would rate all of them as torsadogenic.

Fig. 26

left: MICE model for TdP risk prediction based on the ratio of $\text{IC50}_{\text{hERG}}$ and $\text{IC50}_{\text{CaV1.2}}$. The model assigns high TdP risk to the compounds Dofe42 and Dofe44, though they do not induce EADs (or other types of arrhythmia) in hiPSC-CMs even at high dosage. Data (black dots) for building the logistic regression model (black solid line) are taken from [109], IC50 data for Dofe42, Dofe44 and Dofe54 are taken from [225]. Figure previously unpublished. right: CiPA model prediction of the probability for a drug to have low versus high or intermediate TdP risk, calculated for all of the 28 CiPA drugs of Table 1. Five drugs with high or intermediate TdP risk are falsely classified to be of low TdP risk as they fall below the threshold (red dashed line). Figure credit [16] (reproduced under CC-BY license)

Within the CiPA initiative, logistic regression models for proarrhythmic risk prediction have been tested [16], with up to seven direct features describing electrophysiological responses of spontaneously beating or optogenetically paced [193] hiPSC-CMs to the 28 drugs of Table 1. One model finally chosen reads as

$$\text{logit}(\text{P}) = \beta _{0} + \beta _{1} \cdot \text{F}_{1} + \beta _{2} \cdot \text{F}_{2} + \beta _{3} \cdot \text{F}_{3} $$

with $\beta _{0} = -0.1311$, $\beta _{1} = 0.6583$ (in case of EADs of type A, see Fig. 7 (left)), $\beta _{2} = 0.00687$ and $\beta _{3} = 0.0232$, where $P$ is the probability of a drug to be of high or intermediate TdP risk. Feature $\text{F}_{1}$ is categorical and indicates whether the compound induced arrhythmia-like events at any concentration across all replicates, feature $\text{F}_{2}$ is the maximum compound induced repolarization change observed at any concentration, and feature $\text{F}_{3}$ denotes the drug induced repolarization change at the clinical maximum concentration, with both $\text{F}_{2}$ and $\text{F}_{3}$ averaged across all replicates. A threshold of low TdP risk versus high or intermediate TdP risk was set at 0.8, see Fig. 26 (right), such that all of the drugs in the low-risk category of Table 1 correctly fell below it. However, also one drug from the high TdP risk category and four drugs from the intermediate TdP risk category fell below the threshold.

Examples of misclassification, as those shown in Fig. 25 and in Fig. 26, motivate the use of mechanistic models of cardiac electrophysiology in order to build more accurate risk classifiers that (also) involve features derived from model simulations, see Sect. 5.