The hERG channel base-model used in this study is the model by Li et al. (2016), as shown in Figure 1A in black. It is a six-state Markov model with two inactivated closed (IC) states, two closed (C) states, an inactivated open (IO) state, and an open (O) state. Scaling the ionic conductance of the hERG channel base-model, which is equivalent to inhibiting the hERG current, gives the conductance scaling drug block (CS) model. Extension with drug bound states of open bound (O*) state, inactivated open bound (IO*) state, and closed bound (C*) state, to the hERG channel base-model gives the SD model, as shown in Figure 1A (both black and blue components).

The SD model is based on the model by Li et al. (2017). The transition rate parameters for the hERG channel base-model are taken from Li et al. (2016), Li et al. (2017). The pharmacodynamic component is described by the drug bound states of open bound (O*) state, inactivated open bound (IO*) state and closed bound (C*) state, and their respective transition rates as shown in Figure 1A.

In this model, compounds can bind to the O or IO states of the hERG channel and transition to the O* or IO* states, respectively, with a binding rate (K _b ) given by K b = K u × E max D , (1) = K u × K max × D n D n + E C 50 n , (2)  Where K _u is the unbinding rate and E _max is a sigmoid function describing a drug’s response at concentration D. The function is scaled by K _max, the maximum response of a drug when it is at its saturating concentration. Finally, n is the Hill coefficient and EC₅₀ is the drug concentration when the binding rate is at 50% of its maximum. The compound can be “trapped” in the C* state with a trapping rate (K _t ) that is fixed at K _t = 3.5 × 10⁻⁵ ms⁻¹ (Li et al., 2017). The rate at which the drugs get “untrapped” (K _n ) is defined as K n = K t × X V , (3) = K t × 1 1 + e − V − V half - trap / 6.789 , (4)  where V _half−trap is the membrane voltage when K _n is half the trapping rate K _t . According to Li et al. (2017), the function X(V) is adapted from the hERG channel’s steady state activation of the O’Hara et al. (2011) model. A trapped drug can be modelled with either a high trapping rate (K _t ) or a low untrapping rate (K _n ), and in this case V _half−trap is the only parameter that controls the ratio of trapping to untrapping rate. Finally, the transition rate from the IO* state to the IO state, K _r , is defined to maintain the microscopic reversibility of the system.

To reveal how the SD model alters I _Kr and the AP, we compared the SD model with a simple instantaneous full-state drug block model, the CS model. To ensure that the SD model and the CS model are comparable, we use simulated data from the SD model to calculate the appropriate scaling of the CS model to ensure that the peak I _Kr of both models is equal, as detailed below.

The CS model is equivalent to modelling the current inhibition by scaling the ionic conductance thereby reducing the overall flow of ions across the channel. This model assumes that the half maximal inhibitory concentration (IC₅₀) of a drug on the hERG channel does not change with the experimental protocols. We note that the IC₅₀ of a drug can be protocol-dependent as demonstrated in the later sections, and we acknowledge the differences in the results. However, investigation into such a dependence is not the aim of this study, therefore we use Milnes’ protocol (Milnes et al., 2010, see Section 2.6) to simulate the models as in Li et al. (2017).

The ionic conductance of the hERG channel was scaled to capture the drug response, such that the peak current of the CS model under Milnes’ protocol matched that of the SD model. A schematic of the fitting procedure is shown in Figure 2. We assumed that the SD model captures the experimental data presented in Li et al. (2017), which can be used as a reference point. For each “synthetic drug” (see Section 2.5), the SD model I _Kr was simulated at a range of drug concentrations using Milnes’ protocol. The peaks of these I _Kr were normalised and then fitted to the Hill curve. The Hill curve parameters—the Hill coefficient and IC₅₀—were determined by minimising the mean squared error between the peak currents and the Hill curve. The ionic conductance of the CS model was then scaled with the drug response output obtained from the Hill curve function. The peak hERG current of the CS model will then have the same peak current as the SD model for a given drug and concentration, as required.

The optimisations of the Hill curves to the peak of hERG currents were performed with the covariance matrix adaptation-evolution strategy (CMA-ES) algorithm (Hansen et al., 2003) in the PINTS software (Clerx et al., 2019).

To study the effect of the hERG channel model at the level of APs, the AP model by Dutta et al. (2017) was used. The model was based on the Li et al. (2017) model, and to allow better quantification of each individual current on the AP, the ionic conductances of five currents were scaled (Dutta et al., 2017), including I _Kr, the slow rectifier potassium current (I _Ks), the inwardly rectifying potassium current (I _K1), the L-type calcium current (I _CaL), and the late sodium current (I _NaL). The AP model, when its hERG channel component is the SD model, is referred to as the AP-SD model; the AP model with its hERG channel component replaced with the CS model is referred to as the AP-CS model.

The binding dynamics on the hERG channel of all of the twelve drugs in Li et al. (2017) were taken as the synthetic drugs in this study; the parameters of the 12 synthetic drugs are provided in the Supplementary Material. Here, dofetilide and verapamil were chosen as examples to compare trapping and non-trapping properties of the hERG channel, respectively. Although verapamil is a multi-ion channel blocker, we considered only its effect on the hERG channel, so any multi-channel effects on the AP were deliberately neglected, reflecting the focus of this study being only on the hERG channel binding effects on the AP. The synthetic drug with parameters describing dofetilide is referred to as example drug T, while the synthetic drug with parameters describing verapamil is referred to as example drug N. Only these two drugs out of the 12 synthetic drugs tested are shown in the main results, and the results of the remaining synthetic drugs are provided in the Supplementary Material.

Four voltage protocols were used in this work (Figure 6A). Firstly, Milnes’ protocol as modified by Li et al. (2017) from that of Milnes et al. (2010). This modified Milnes’ protocol aims to identify trapped and non-trapped drugs. Experimental data obtained from the stimulation with this protocol was used to fit the parameters of the SD model. Therefore, this modified Milnes’ protocol was used as a reference for the model comparison.

The remaining voltage protocols, Pneg80, P0, and P40, were taken from Gomis-Tena et al. (2020). These protocols were designed so that the channel’s state occupancy of a certain state is maximised. When the hERG channel model is stimulated by the Pneg80, P0, and P40 protocols, the channel will most likely be in the closed, open, and inactivated states, respectively, as discussed in Gomis-Tena et al. (2020). Details of all protocols are given in the Supplementary Material.

The behaviour of the drug-related component of the SD model is governed by its transition rates, which in turn are defined by the five drug-dependent parameters. To understand the role of each drug-related parameter in the SD model, we performed a sensitivity analysis on the APD₉₀ difference between the AP-SD and the AP-CS model.

First, we simplified the SD model by normalising the drug concentration D input to the EC₅₀ value. The normalised drug concentration is therefore defined as D ̃ = D E C 50 . (5)  This reduces the number of parameters to four.

Second, we observed that the Hill coefficient n affects only the range of drug concentration where AP prolongation and early after depolarisation (EAD)-like behaviour are observed. To confirm this, we varied the value of n for each synthetic drug and monitored the APD₉₀ differences between the two models. The APD₉₀ differences were quantified by the root mean square difference (see Section 2.8). The root mean square difference (RMSD) of APDs fluctuated within a range of only 25 ms as n was varied (see Supplementary Material). The model was thus further simplified to three parameters, namely, the V _half−trap, K _max, and K _u . We then performed a detailed sensitivity analysis on this 3-dimensional space.

In the sensitivity analysis, we repeated the model comparison procedure as shown in Figure 2 for a range of parameter values in the simplified SD model. The three parameters of the simplified SD model describe a hypothetical drug which is termed a “virtual drug”. Sweeping through the parameter space effectively performed a virtual drug screening for all possible drug binding kinetics of hERG, within the SD model.

The two models, the AP-SD model and the AP-CS model, were considered to be similar if the RMSD of the APD₉₀ values are small. The parameter space around the boundary surface, where the RMSD values between the two models are less than 30 ms, were sampled more densely to obtain higher resolution.

The RMSD and the mean difference (MD) between the APD₉₀ of the AP-SD and the AP-CS models were used to quantify the model comparison. The RMSD measures the magnitude of the APD₉₀ difference between the models, and is given by R M S D = ∑ i = 1 N SD i − C S i 2 N , (6)  where N is the total number of data points, excluding data points where EAD-like behaviours are observed, and SD and CS are the APD₉₀ values of the AP-SD model and the AP-CS model, respectively.

The MD measures the actual difference between the two models, M D = ∑ i = 1 N SD i − C S i N . (7)  A positive value of MD implies that the sum of the APD₉₀ values of the AP-SD model are higher than the AP-CS model, while a negative value implies the opposite.

The signed RMSD was used to indicate the magnitude and direction of the APD₉₀ difference between the AP-SD model and the AP-CS model: R M S D signed = R M S D × M D | M D | . (8)

All simulations, including voltage-clamp and current-clamp experiments, were run with Python 3.8 using Myokit 1.33 (Clerx et al., 2016) with the CVODE solver (Hindmarsh et al., 2005). The absolute tolerance and relative tolerance were set to 10⁻⁷ and 10⁻⁸, respectively.

Figures 1B, C show that the SD model can capture the trapping properties of a drug, observed in the lack of recovery from block and the state occupancy of drug bound states. However, it is not clear whether a simpler AP-CS model could replicate the AP prolongation by the AP-SD model that has state dependency and a trapping mechanism included. We systematically compared the APs of the two models for the example drug T and the example drug N using the methods described in Section 2.3; Figure 2.

Figure 3A shows the I _Kr of both the SD model and the CS model with an example drug T present there throughout under Milnes’ protocol. By design, the peaks of the I _Kr of the two models are the same at various drug concentrations of the example drug T, which give the same dose-response effect. The APs of the two models (the AP-SD model and the AP-CS model) with the example drug T are given in Figure 3B, showing two pulses of the AP at steady state, together with their corresponding I _Kr. As the drug concentration increases, the amplitude of the I _Kr decreases, and the APD increases. The two models show similar AP behaviours, although the AP-SD model simulated I _Kr at 10 nM of example drug T has a bigger amplitude and a shorter period when the I _Kr is positive. The differences give rise to a slightly shorter APD because I _Kr aids in the repolarisation of APs. Moreover, both models show an EAD-like behaviour at the same drug concentration (300 nM) of the example drug T.

In Figure 3C, the APD₉₀ values of APs simulated from the AP-SD model and the AP-CS model for an example drug T were calculated. APs that show EAD-like behaviours are indicated with an asterisk. The APD₉₀ values of both models are similar, consistent with the APs shown in Figure 3B. Additionally, the qNet values, a proarrhythmic risk marker, of example drug T shown in Figure 3D are similar for both models. Therefore the CS model can replicate the SD model for an example drug T at steady state.

We further compared the APs of the models for an example drug N, a “non-trapped” drug; the non-trapping phenotype of the drug is captured by the small value of V _half−trap. Figure 4A shows the I _Kr of the two models under the effect of the example drug N. The I _Kr simulated by the SD model show dips after the initial increase at drug concentrations higher than or equal to ∼300 nM. Since the peak current is used to define the inhibition level of a drug, which is commonly done for IC₅₀ calculations, the total amount of I _Kr simulated by the CS model is higher than the SD model.

Figures 4B, C compare the APs and their corresponding I _Kr of the AP-SD model and the AP-CS model for the example drug N. The I _Kr and the APs are significantly different between the two models. At high drug concentrations the AP-SD model shows a longer APD and a lower I _Kr amplitude than the AP-CS model. The AP of the AP-SD model is prolonged more than the AP-CS model and displays EAD-like behaviours at lower drug concentration. The example drug N, for which the untrapping rate is higher and the binding rate is lower than the example drug T, generated lower APD₉₀s with the AP-CS model. Similarly, the qNet values of example drug N (Figure 4D) are different between the two AP models. Example drug T and example drug N are taken as examples to show the difference in APD₉₀ and qNet for drugs with varying trapping tendency. The same analysis was repeated for other drug compounds (see Supplementary Material).

The main difference between the two drug block models is the inclusion of the trapping mechanism in the hERG channel model. The feature of this mechanism is the accumulation of the drug compound in the channel due to “trapping”. Figures 5A–D show the progression of the APs and their I _Kr for the example drug T and the example drug N with the two AP models. In Figure 5A, the I _Kr decreases gradually with time for the AP-SD model with the example drug T, and the AP prolongation increases. In contrast, the I _Kr and the APs for the example drug N in the AP-SD model are relatively stable, as shown in Figure 5B. The AP-SD model displays a shorter transient phase with the example drug N than with the example drug T, which is due to the differences in the trapping phenotype of these synthetic drugs. Figures 5C, D compare the AP-CS model with the example drug T and the example drug N, respectively, and show that the APs and the I _Kr do not change much between pulses.

The APD₉₀ values from Figures 5A–D are given in Figures 5E, F for the example drug T and the example drug N, respectively. For the example drug T (Figure 5E), the AP-SD model showed a longer transient phase, then stabilised to generate the same APD₉₀ as the AP-CS model. By contrast, the two models with the example drug N (Figure 5F) had a short transient phase with different APD₉₀ values at steady state.

In all model comparisons shown in Figures 3, 4, we have assumed that the IC₅₀ of a drug is independent of the protocol used to measure its response. Thus, the IC₅₀ can be used to replicate drugs effect in any drug block models (Mirams et al., 2011). However, the assumption does not always hold (Kirsch et al., 2004; Yao et al., 2005). Here, we demonstrate the change in the APD₉₀ when different protocols are used to generate the I _Kr for drug characterisation. Figure 6A shows Milnes’, Pneg80, P0, and P40 protocols that are used to stimulate the SD model. The normalised peak of the I _Kr simulated with each of the given protocols for both the example drug T and the example drug N are shown in Figures 6B, C, respectively. The example drug T (Figure 6B) yields similar dose-response curves for all of the protocols, except for the P40 protocol, while the example drug N (Figure 6C) displayed different dose-response curves for all protocols. Using the Hill curves in Figure 6C to characterise the drug effect for the CS model (see Section 2.3; Figure 2), the APs of the AP-CS model were simulated and quantified. Figure 6D shows the APD₉₀s of the simulated APs with an example drug N. While using the Pneg80 protocol gave similar APD₉₀s to using Milnes’ protocol, the use of P0 and P40 protocols generate higher APD₉₀ values, causing EAD-like behaviours to appear at lower concentrations of the example drug N. Example drug T and example drug N are taken as examples. The Hill curves for all the protocols with the other drug compounds are given in the Supplementary Material.

A sensitivity analysis was performed on the simplified SD model as described in Section 2.7, where the SD model was simplified to just three drug-related parameters—V _half−trap, K _max, and K _u —without loss of generality. Given a virtual drug from the parameter space, the procedure shown in Figure 2 was repeated (Section 2.3). The APD₉₀ differences between the AP-SD model and the AP-CS model were measured by the signed RMSD (Equation 8).

Figure 7A shows the signed RMSD for combinations of the three parameters. Each point represents a virtual drug, while its color indicates the APD₉₀ difference between the AP-SD model and the AP-CS model when the virtual drug is added. The majority of the virtual drugs resulted in a positive signed RMSD value, indicating that the AP-SD model generated APs that had longer durations than the AP-CS model. It is also shown that the signed RMSD shifts from positive to negative as K _max decreases. The virtual drugs are generally predicted to cause shorter AP prolongation when using the AP-CS model.

The parameter space where the two models give a similar APD₉₀ value (RMSD < 30 ms) is limited, as shown in Figure 7B. For these parameters, the AP-SD model and the AP-CS model are interchangeable when we measure only the APD. The triangles show the positions of the synthetic drugs in the parameter space. Synthetic drugs that lie within the boundary surface are listed in the Supplementary Material.

The results show that V _half−trap plays a small role in affecting the differences between the AP-SD model and the AP-CS model. K _max and K _u are the major driving forces for the changes in the APD₉₀ difference between the two models. Furthermore, at low values of K _max, K _u does not affect the APD₉₀ difference as much; for large K _max, the two models differ the most around K _u ≈ 5 × 10⁻⁴. A different viewing angle of Figure 7B is available in the Supplementary Material.