The impact of uncertainty in hERG binding mechanism on in silico predictions of drug-induced proarrhythmic risk

Background and Purpose Drug-induced reduction of the rapid delayed rectifier potassium current carried by the human Ether-à-go-go-Related Gene (hERG) channel is associated with increased risk of arrhythmias. Recent updates to drug safety regulatory guidelines attempt to capture each drug’s hERG binding mechanism by combining in vitro assays with in silico simulations. In this study, we investigate the impact on in silico proarrhythmic risk predictions due to uncertainty in the hERG binding mechanism and physiological hERG current model. Experimental Approach Possible pharmacological binding models were designed for the hERG channel to account for known and postulated small molecule binding mechanisms. After selecting a subset of plausible binding models for each compound through calibration to available voltage-clamp electrophysiology data, we assessed their effects, and the effects of different physiological models, on proarrhythmic risk predictions. Key Results For some compounds, multiple binding mechanisms can explain the same data produced under the safety testing guidelines, which results in different inferred binding rates. This can result in substantial uncertainty in the predicted torsade risk, which often spans more than one risk category. By comparison, we found that the effect of a different hERG physiological current model on risk classification was subtle. Conclusion and Implications The approach developed in this study assesses the impact of uncertainty in hERG binding mechanisms on predictions of drug-induced proarrhythmic risk. For some compounds, these results imply the need for additional binding data to decrease uncertainty in safety-critical applications.

I Kr (Sanguinetti et al., 1995).The hERG channel is highly susceptible to blockage or functional inhibition by a variety of pharmaceutical small molecules (Sanguinetti & Tristani-Firouzi, 2006;Vandenberg et al., 2012), we simply refer to these as 'compounds' throughout.Reduction of the rapid delayed rectifier potassium current, I Kr , can lengthen the action potential (AP), which is associated with increased risk of cardiac arrhythmias, including Torsade de Pointes (Curran et al., 1995;Heist & Ruskin, 2010), although this risk is strongly modulated by multichannel block (Mirams et al., 2011).In vitro studies including hERG are part of the regulatory guidelines for proarrhythmic risk assessment of novel compounds (ICH, 2005).Recently, the Comprehensive in-vitro Proarrhythmia Assay (CiPA) initiative and updates of regulatory guidelines consider a more nuanced characterisation of a drug's proarrhythmic potential, by combining in vitro assays for multiple ion channels with in silico simulations.For hERG, this approach examines the details of each compound's binding mechanism, particularly the degree to which each compound is trapped (cannot unbind) when the channel closes (ICH, 2022;Li et al., 2019;Sager et al., 2014).
The inner cavity of the hERG channel is its principal drug compound binding site (Mitcheson, 2008).The molecular structure of the channel suggests that compounds only bind when channels are not closed (Butler et al., 2020), consistent with the observation that most of the compounds do not bind at negative (resting/repolarised) membrane potentials when hERG channels are in a nonconducting closed state (Li et al., 2017).After binding to the channel, compounds may also unbind when channels are in different states and/or at different voltages (Mitcheson et al., 2000;Thouta et al., 2018).However, the unbinding process for some compounds can be impeded when the channels close, and the compounds remain bound, or 'trapped' within the central cavity (Mitcheson et al., 2000;Stork et al., 2007;Thouta et al., 2018;Windisch et al., 2011), such as bepridil (Pareja et al., 2013) and dofetilide (Milnes et al., 2010), as opposed to, for example, cisapride (Milnes et al., 2010) and verapamil (S.Zhang et al., 1999), which unbind when the chancels close.According to the modulated receptor hypothesis, the difference in affinity determines the preferential binding of a compound to one of the states (Carmeliet & Mubagwa, 1998;Hille, 1977;Hondeghem, 1987;Hondeghem & Katzung, 1977), and its special case, the guarded receptor hypothesis, suggests the possibility compounds bind to a particular state only (Starmer & Courtney, 1986;Starmer et al., 1990Starmer et al., , 1991)).Both of these have been applied to computational (in silico) models of hERG binding (Gomis-Tena et al., 2020;Lee et al., 2017;Thurner et al., 2014;Veroli et al., 2013).
To capture all the possible consequences of state-dependent binding for a new compound within in silico models, all the above possibilities for drug binding to hERG should be considered.However, this raises the question of whether one size fits all-does the hERG binding model used in CiPA (Li et al., 2017) account for all the binding properties of interest, and do the model parameters reflect the underlying binding mechanisms?Furthermore, how sensitive is the proarrhythmic risk classification metric to any uncertainty in the binding mechanisms and also the underlying physiological hERG model (i.e., the description of channel gating in the absence of compounds)?
In this study, we propose a method for considering a wide range of drug/ion channel binding models for hERG and rule out those that do not plausibly fit the available experimental data, resulting in an ensemble of models for each compound.The ensemble represents uncertainty in each compound's hERG binding mechanism.We then apply this method to the CiPA v1.0 validation study (Li et al., 2019) and examine how robust risk classifications are to this uncertainty in hERG binding mechanism.The same approach can be used whenever fitting drug binding data to ensure that we do not ignore this important source of uncertainty.

| METHODS
In this study, we designed a set of possible pharmacological binding models for the hERG channel that account for the known or postulated binding mechanisms (which we will describe in Sections 2.1 and 2.2).These models were calibrated to voltage-clamp electrophysiology data under the Milnes et al. (2010) protocol (Section 2.3.1)from Li et al. (2017Li et al. ( , 2019)), and plausible binding models that can explain the data well were selected (Section 2.3.2).We then used these binding models with different hERG physiological models to calculate the 'CiPA v1.0' risk metrics to assess the possible impact of different hERG binding mechanisms on such predictions (Section 2.4).

What is already known
• Compounds interact with the hERG channel with different state dependence and trapping properties.
• Proarrhythmic risk depends on the degree of trapping in the CiPA in silico modelling predictions.

What does this study add
• For many compounds, multiple binding mechanisms are consistent with data gathered under the CiPA guidance.
• Risk predictions are sensitive to these hERG binding mechanisms but more robust to gating model.

What is the clinical significance
• Clinical risk predictions for some compounds vary between plausible hERG binding mechanisms.
• Narrowing down plausible binding mechanisms will be needed to reduce this source of uncertainty.

| hERG physiological models
The basic physiological models of hERG describe the drug-free control of voltage-dependent gating behaviour for the rapid delayed rectifier potassium current (I Kr ) at physiological temperature and conditions.
Two physiological models are used in this study, with the transition rates following p i expðp j VÞ, where p i and p j are physiological model parameters taken from the literature.The first one, physiological model A, used in CiPA studies is shown in Figure 1a (Dutta et al., 2017;Li et al., 2017).The second model, physiological model B, is a symmetric four-state model (Figure 1b) with transition rates k 1 to k 4 , which is equivalent to a Hodgkin & Huxley-style model with one activation gate and one inactivation gate (the 37 C average model in Lei, Clerx, Beattie, et al., 2019;Lei, Clerx, Gavaghan, et al., 2019).

| Pharmacological binding models for hERG
To account for various known or proposed mechanisms for compounds binding to the hERG channel, a set of pharmacological binding models was designed as shown in Figure 2. The hERG physiological model is indicated in black, with common states I and O shown (on the right-hand side of each model from Figure 1), with dots representing the rest of either model.
In Figure 2, Model 1 represents a compound that can bind to both open (O) and inactivated (I) states, the drug is not trapped in the binding pocket, and the open and inactivated states share the same binding and unbinding rates (shown in the same colour).The binding rate is assumed to be proportional to the drug concentration ½D raised to the power of a Hill coefficient, n, and the association rate constant is k on ; the unbinding rate is assumed to be a constant k off .Model 2 allows binding only to the open state; the inactivated version of Model 2 (Model 2i) allows binding only in the inactivated state, using the guarded receptor hypothesis (Escobar et al., 2022;Lee et al., 2019).We have not observed any evidence that existing drugs bind at the resting potential ( ≈ 80 mV) when the channel is closed, therefore we have excluded binding models that bind only to the closed state.Model 12 is taken from Li et al. (2017).Here, rather than the drug binding rate being linearly proportional to ½D n , instead it saturates, and this is represented using a Hill equation: where EC 50 is a half maximal effective concentration, that is, when  2); they also assume a fixed trapping rate k trap ¼ 3:5 Â 10 À5 ms À1 , with V 1=2; trap altering the degree of trapping (Li et al., 2017).
We also included two additional models: Models 0α and 0β, as basic and standard 'conductance block' models for drug effects.We consider an 'all-state-blocker' model that has binding and unbinding rates, k on ½D n and k off , for all channel states.The degree of block is then independent of state occupancy and is equivalent to having a modulating fraction of unavailable channels or 'b gate' such that we scale (multiply) the drug-free current by ð1 À bÞ.b itself follows the equation with an initial condition of b ¼ 0 at the time a compound is first introduced, we call this Model 0β.Model 0α is a simpler conductance The ratio k off =k on in Model 0β is known as the dissociation constant and is equivalent to IC n 50 in the Hill equation in Model 0α.All the binding model parameters that need to be calibrated on a compoundspecific basis are listed in Table 1.

| Electrophysiology data
The voltage-clamp electrophysiology data were taken from the openly available 28 compound CiPA training and validation studies (Li et al., 2017(Li et al., , 2019) ) where manual patch-clamp experiments were performed on HEK293 cells (CLS Cat# 300192, RRID: CVCL_0045) stably expressing hERG1a subunit at 37 C-see Data Availability section.
Data were originally collected using a modified Milnes' protocol (Milnes et al., 2010); the protocol had 10 sweeps measured with 10-ms time-point interval, and each repeat consisted of 1 s (with leak step) at the resting potential (À80 mV), followed by a long voltage step to 0 mV for 10 s, before returning to the resting potential for 14 s.For the details of the experiments, please refer to Li et al. (2017).

| Calibration of pharmacological binding models
We calibrated/fitted each of the pharmacological binding models in Figure 2 with each of the hERG physiological models in Figure 1 independently to the voltage-clamp electrophysiology data for each compound.The available data have been normalised to the control current, giving the fraction of unblocked current, which reveals the change in current due to the drug binding.Furthermore, only the data at the voltage step 0 mV (time interval 1.1-11 s) for the 10 sweeps were used; the remaining data at the resting potential (À80 mV) will have almost no current, giving little information about the drug binding.All models, except Model 0α, were calibrated by minimising the root-mean-square difference (RMSD) of the percentage current between the model output and the mean experimental data from 10 sweeps for four different concentrations of a compound, giving 4 Â 10 Â 990 ¼ 39600 data points for each compound.The optimisation was performed with logarithmic-transformed model parameters, except V 1=2; trap in Models 12 and 13.This allows a better search across a wide parameter range, especially for the rate parameters.The optimisation was performed using the covariance matrix adaptationevolution strategy (CMA-ES) algorithm (Hansen, 2006) in PINTS (Clerx et al., 2019) and was repeated 10 times from different initial guesses sampled from wide boundaries (k on , k off ½10 À7 , 1 ms À1 , k trap , k untrap ½10 À9 , 10 3 ms À1 , kon ½10 À6 , 10 9 ms À1 Á nM Àn , EC n 50 ½1, 10 9 nM n , V 1=2;trap ½0, 200 mV, and n ½0:2, 2) to ensure we arrived at a global minimum.The calibration was repeated for all 28 compounds listed in Li et al. (2019).
Model 0α is a special case of Model 0β where drug binding is approximated as instantaneously reaching steady state.
Therefore Model 0α should not be calibrated to transient experimental data.Hence, the parameters of Model 0α were directly taken as the steady state of Model 0β-IC n 50 of Model 0α in Table 1 was calculated as the dissociation constant k off =k on of Model 0β, as Equation (3) suggests.
However, since each pharmacological binding model in Figure 2 was designed to model a specific mechanism, not all the models are expected to be able to describe all the compounds that were tested.
Therefore each calibrated binding model was assessed by comparing its fitted RMSD to the RMSD of bootstrapped samples of the data which were computed as follows: since each of 10 repeated T A B L E 1 The model parameters for all the binding models in Figure 2. The asterisk indicates the rate parameter was determined by microscopic reversibility.

Model
Binding parameters Trapping parameters experiments at each concentration was independent, for one compound across all four concentrations, there are 10 4 permutations.We randomly selected (with replacement) one trace out of the ten available for each concentration to obtain a set of four traces for all concentrations; the RMSD of this set to the mean of the experimental data was computed, giving the RMSD of a 'bootstrapped' data trace.
This process was repeated 1000 times to get a range of RMSDs based on these bootstrapped samples of the data.This essentially compares how well the models fit to the mean data relative to an individual data trace's fit to the mean data, providing a compound-dataset-specific measure of goodness-of-fit.
Although if the only source of error were normally distributed noise, we could compare the model fits against the standard error of mean, this would probably rule out all the models (including the CiPA v1.0 model).Therefore, a more relaxed threshold was chosen to allow for some degree of imperfection in our models of both where β a 0 and β b 0 are the intercepts, and β 1 the linear coefficient (of the torsade metric), of the linear equations that map the torsade metric through a logistic function to the cumulative probabilities of the risk categories within the logistic regression model.Because the action potential shape changes slightly with a different physiological hERG model, then q net at control and derived thresholds also change.The metric values for each physiological model were normal- where x can be either q net or the torsade metric score, and m denotes the physiological model (either model A or B), such that the boundary between high-risk and intermediate/low is at 0 and between low-risk and high/intermediate is at 1 to facilitate comparison of classification between physiological models.

| Nomenclature of targets and ligands
Key protein targets and ligands in this article are hyperlinked to corresponding entries in https://www.guidetopharmacology.org, and are permanently archived in the Concise Guide to PHARMACOLOGY 2021/22 (Alexander et al., 2021).

| Multiple binding mechanisms can explain the same drug effects
All binding models (Figure 2) were calibrated to the experimental data for each compound independently and compared against the maximum RMSD for bootstrapped samples of the data.Figure 3 shows the calibrated binding models (left) and their RMSDs to the mean experimental data (right) across all binding models, for three The other two example compounds in Figure 3, terfenadine and verapamil, had a weaker to no trapping tendency compared to dofetilide; that is, current recovers significantly between 0 m V pulses, suggesting unbinding rather than trapping at À80 mV.Interestingly, more binding models failed to fit to the experimental data of these two drugs.For terfenadine, a weakly trapped drug (Kamiya et al., 2008;Stork et al., 2007;T. Yang et al., 1995) with slow binding rates (Li et al., 2017), only Models 7 and 8 (variants of the independent open and inactivated binding model) and Models 6 and 10-13 (variants of flexible-trapping models) were able to explain its drug effect as well as the bootstrap samples of the data, although a few more nontrapping models were merely marginally ruled out.For verapamil, a nontrapped drug (S.Zhang et al., 1999), our approach successfully disqualified all simple-trapping models (Models 4, 5, 5i, and 9), as well as the conductance blocking and the all-state binding models (Models 0α and 0β).However, it left the non-trapping models (Models 1-3, 7, and 8) and the flexible-trapping models (Models 6 and 10-13) as plausible candidates for the binding mechanism.We observed no obvious pattern between the proarrhythmic risk of drugs and the type of binding models ruled implausible.A few more nontrapped drugs were identified simply based on the ruling out of the simple-trapping models (Models 4, 5, 5i, and 9) as in the case of verapamil (Figure 3), such as cisapride (Milnes et al., 2010) and droperidol (Stork et al., 2007;Windisch et al., 2011).For most of the drugs, multiple binding models were able to explain the observed drug effects using the experimental data collected through CiPA-Milnes protocol.

| Inferred binding rates can be strongly dependent on binding mechanism
Studying drug binding kinetics, such as the binding and unbinding rates, is important for understanding the drug effects on the channels and predicting behaviour in new situations.Taking verapamil as an example, the inferred binding rate parameter k on was similar for most of the successfully calibrated binding models (Figure 3).The two binding rate parameters k on;O (empty marker) and k on,I (filled marker) for the plausible models-Models 7, 8, and 10-were similar too, and their average roughly equalled k on of the other models.The similarity suggested that the independence assumption may be superfluous in this case, because if the two inferred rates were the same, the models would be equivalent to those without the extra degree(s) of freedom (Models 1, 3, and 6).
Similarly, the unbinding rates k off were similar for all binding models, except the k off;O (empty marker) and k off;I (filled marker) in Model 6; all of the inferred Hill coefficients n were within the range of 1.3-1.8.It is worth noting that the k off of Model 12 and the CiPA v1.0 model (the same binding model structure but under different calibration schemes) shared similar inferred values for verapamil but their kon were not shown due to different units.
For metoprolol, the inferred k on were heavily binding-modeldependent, with a coefficient of variation (defined as the ratio of the standard deviation to the mean, for the plausible models) of 607%, as compared to 141% for verapamil.The other two parameters k off and n were similar across the calibrated binding models.Unlike verapamil, the inferred k off of Model 12 and the CiPA v1.0 model were inconsistent, using the same set of experimental data but with different calibration schemes; see also mexiletine in the Supporting Information.
The differences in the inferred parameters may raise questions about the calibration data and the complexity of the model structure used, leading to potential parameter unidentifiability issues (Whittaker et al., 2020; see also Section 4).

| q net with different binding mechanisms can diverge
Thus far, we compared the binding models of hERG when calibrated to voltage-clamp experimental data for various compounds.Here, we show the results of our investigation on the impacts of these binding models on APs and risk simulations (Section 2.4), whilst accounting for multichannel effects, following Li et al. (2019).Figure 6 shows the metric q net at various C max levels for all binding models, with implausi- Various degrees of q net spread were observed across the binding models for different compounds.
For example, diltiazem is a strong I CaL blocker relative to its I Kr effects, therefore the drug effect on APs measured through q net reflects mostly the drug block of I CaL and the effects of different I Kr binding models are insignificant.Nifedipine and verapamil are also multi-channel blockers of I CaL and I Kr , with similar block levels for each current.However, the uncertainty in their resulting q net predictions were drastically different-nifedipine showed tight q net predictions across the binding models, whilst verapamil produced a wide spread.Furthermore, in general, we also observed an increasing spread of the q net predictions from plausible binding models at higher C max levels; verapamil is one of the most obvious examples.This phenomenon showed the differences between the binding mechanisms were amplified with higher concentrations of the drug, and at the top concentrations verapamil spanned all risk categories.However, the torsade metric uses only 1-4Â C max of q net -the metric that was used to train the classification model (ordinal logistic regression model) for producing the decision boundaries of the risk categories, and we examine risk predictions in this range in the next section.AP for high risk in red.The plausible binding models can result in metric predictions with high uncertainty, even spanning multiple risk categories; for example, domperidone predictions span all three risk categories.S1.The shift in the boundaries was consistent with the change in the net charge carried by physiological model B (and alteration to AP shape, which also affects other currents), which resulted in a new control value of q net .The two new decision boundaries are used for physiological model B in Figure 8b,c.

| The effect of hERG physiological model on risk classification is subtle
In Figure 8b,c, to compare the differences in the metric predictions between the two physiological models, the metrics predicted by each model were normalised to its decision boundaries using Equation ( 6 F I G U R E 6 q net , the net charge carried by currents active in plateau and repolarisation over one beat of AP (I Kr , I NaL , I CaL , I to , I Ks , and I K1 ) at various multiples of C max for all binding models with six example drugs.Here, diltiazem, nifedipine, verapamil, dofetilide, cisapride, and bepridil are shown, revealing a spectrum of behaviours from the binding models.Dashed horizontal lines indicate the decision boundaries for the low-risk (green), intermediate-risk (blue), and high-risk (red) categories for the torsade metric score (average q net at 1-4Â C max ) in Li et al. (2019).The implausible models are shown as transparent dotted lines, and the CiPA v1.0 model is shown as black dashed lines.

F I G U R E 7
The torsade metric predictions, the mean q net of 1-4Â C max , of all binding models with physiological model A for all compounds, and the corresponding AP predictions for seven example drugs.Models 0α, 7 and 11 were used in Windley et al. (2016) to study the effects of compounds such as cisapride.
In Figures 3 and 4, we demonstrated that our approach can be used to distinguish some of the simpler binding mechanisms based on the CiPA-Milnes protocol data for some compounds, such as terfenadine and verapamil; resulting in plausible binding mechanisms consistent with the literature.It was also able to highlight certain similarities between some compounds such as bepridil (Kamiya et al., 2006;Pareja et al., 2013) and terfenadine (Stork et al., 2007) where only flexible-trapping models were deemed plausible to explain the observed data; indeed, studies with specifically designed voltage protocols consider these two compounds to be trapped slow-binders (Kamiya et al., 2008).Also, our approach was able to highlight compounds, such as tamoxifen, loratadine, and nitrendipine, where none of the models (not even the CiPA 1.0 model) were able to fit the data satisfactorily.On closer inspection (Supporting Information) the data of loratadine, and nitrendipine showed a slight increase of the (percentage) current over time during the 0 m V pulses, raising potential data quality issues (Lei, Clerx, et al., 2020;Montnach et al., 2021;Raba, et al., 2013) or the need for methods to account for inadequacy of the models (Lei, Ghosh, et al., 2020;Lei & Mirams, 2021) and/or new (un)binding mechanisms to explain this observation; whilst for tamoxifen, there was a more obvious data quality issue for one of the concentrations.In Figure S26, we also included the results of fitting all the binding models whilst assuming the Hill coefficient (number of binding sites) to be n ¼ 1.However, in this case, most of the binding models failed to fit to many of the compounds-being classified as implausible models-suggesting the importance of the extra degree of freedom provided by the Hill coefficient in explaining these observations, although we do not eliminate the possibility of this having been caused by some experimental artefact resulting in a drift/rundown that the analysis approach mistakes for non-saturating or more quickly/slowly saturating hERG block.
It may be better to think of our 'plausible' binding models as 'not implausible', there is no evidence that any individual model is correct, indeed by their nature some contradict others, it is just that they cannot be ruled out yet based on the available data.It is desirable for the calibration data to contain enough information to rule out as many binding models as possible, in an ideal world leaving just a single plausible model.But with the resulting ensemble at present, we can objectively and quantitatively describe a compound's binding mechanism instead of qualitatively classifying the compound as 'trapped', 'binding to open state', and so forth.However, we observed that the data elicited under the CiPA-Milnes protocol, which was originally designed to differentiate between trapped and nontrapped compounds (Milnes et al., 2010), were not able to distinguish between all the possible binding mechanisms, or indeed even whether trapping occurs, for all compounds (Figure 3).This suggests the need for designing better, richer experimental protocols for model calibration and selection (Lei et al., 2023), for instance gathering data on block onset at different voltages (Gomis-Tena et al., 2020;Lee et al., 2019), and to overcome some of the difficulties in measuring fast timecourses of block (Windley et al., 2017).Indeed a range of protocols may needed for assessing fast-and slow-binding compounds that bind via multiple mechanisms.
The model structure is not only useful for identifying the binding mechanisms; the inferred model parameters, such as binding/unbinding rates, for the plausible pharmacological binding model(s) can be used as a proxy to quantitatively assess the binding behaviour and dynamics.
However, we showed that, with the limitations of the calibration data, not all plausible models recover the same binding and unbinding rates, although their ratios, the dissociation constants k off =k on of the models, were more consistent (Figure 5 and Supporting Information).This leads to a more complicated interpretation of the binding properties.
Without being able to identify the correct binding mechanism(s), we would not be able to study the true binding/unbinding rates of the compound.Moreover, for more complex pharmacological models, the inferred parameters may even be subject to the calibration scheme and procedure: For example, the inferred unbinding rate k off of Model 12 and CiPA v1.0 were inconsistent (Figure 5), although our Model 12 gives lower RMSDs (Supporting Information).The different optimal parameters in CiPA v1.0 versus Model 12 fits could be due to the use of a different objective function-there a tailored weighted sum of squares of residuals, and based on fewer optimisation runs which could perhaps make local rather than global optima more likely (Li et al., 2017).We note that this is not the same as the unidentifiability issue between, in our notation, kon and EC 50 discussed in Li et al. (2019), suggesting this is a difficult optimisation problem.Our findings emphasise the need for careful design of experiments when parameterising complex models (Whittaker et al., 2020).
In Figure 4, we observed that constant instant conductance block (Model 0α) was not classified as a plausible model for most of the compounds, as we would expect, the goodness of fits were inadequate and poor (Figure 3).Yet, the torsade metric predictions by Model 0α were not dissimilar to the other models (Figure 7; see also Han et al., 2019;Mistry, 2019), which was likely due to the difference between using transient data (Milnes' protocol) and predicting steady-state (q net /torsade metric) conditions (Farm et al., 2023).However, we believe the torsade metric predictions of Model 0α were acceptable only because of the steady-state nature of the metric.If Model 0α were to be used to predict transient APs under certain changes of conditions or pacing rates, then kinetic effects would be neglected and bad predictions would likely result.
In Figure 7, we also noticed that the degree of variation tended to be larger for drugs in higher proarrhythmic risk categories.This phenomenon was thought to be due to the multichannel effects, which either compensate for the effects of different hERG binding mechanisms or make them insignificant.However, if our multiplebinding-model approach was applied to other types of current, such as I CaL , it may result in a similar level of uncertainty in their binding mechanisms, leading (correctly) to a higher level of uncertainty in the risk predictions than we observed (Figure 7), particularly for compounds that mainly block I CaL (e.g., diltiazem and nifedipine).Nonetheless, such an observation also implies that it is likely to be most important to determine the correct hERG binding mechanism when predicting drugs with higher proarrhythmic risk.Also, given some of the risk prediction spanned multiple risk categories, we advocate efforts to reduce the uncertainty in binding model.Overall, the results highlight the importance of the details of hERG binding mechanisms.
Furthermore, we have studied the effects of two different hERG physiological channel gating models (Figure 1) on proarrhythmic risk predictions.The CiPA v1.0 model uses hERG physiological model A (Li et al., 2017), which was designed to closely match the I Kr in O'Hara et al. ( 2011) which itself was fitted to the native I Kr of a cardiac myocyte-composed of a mixture of hERG1a/1b channels or heterotetramers (Jones et al., 2004;Sale et al., 2008) and any associated (possibly still unknown) subunits and regulatory components, rather than just hERG1a.In contrast, physiological model B was trained using data from a cell line over-expressing hERG1a (Lei, Clerx, Beattie, et al., 2019), which explains why the two models exhibit marked differences in their predicted kinetics (Figure 8A).Since the drug-binding calibration data were also gathered using hERG1a cell lines, we might We further note the potential of using computational methods for studying molecular structural dynamics and relationships in drug binding, which may also provide insight into the binding mechanisms (DeMarco et al., 2021;Munawar et al., 2019;P.-C. Yang et al., 2020).
In future, structural approaches could enhance our approach, in terms of ruling out further binding mechanisms, particularly when all voltagedependent state structures for wild-type hERG1a and 1b are available (Abi-Gerges et al., 2011;Maly et al., 2022;Y. Zhang et al., 2020).
Finally, to our knowledge, this is the first study that attempts to address the question of whether we need a complex pharmacological model that attempts to nest most/all of these binding mechanisms or a set of multiple possible simpler pharmacological models to better represent uncertainty in proarrhythmic risk predictions.In theory, if parameter unidentifiability were not an issue, the two wings should arrive at the same conclusion-for example, when modelling an untrapped compound, the transition rates to the trapping component of the complex model would approach zero, and the non-trapping (simpler) model would be selected as the only plausible model.
However, it may be better to use simpler (fewer parameter) models that can reproduce the underlying binding mechanism given the difficulty of eliciting information-rich data for calibrating one-size-fits-all complex binding models leading to potential parameter unidentifiability issues (Whittaker et al., 2020), given the inevitable presence of some model discrepancy (Lei, Ghosh, et al., 2020) and residual experimental artefacts (Lei, Clerx, et al., 2020).Pragmatically, we would suggest to select and use all the plausible simple models for prediction as demonstrated here-a type of ensemble model prediction which provides an estimate of uncertainty due to model discrepancy (Murphy et al., 2007;Parker, 2013;Tebaldi & Knutti, 2007).
In conclusion, this study has developed an approach to analyse a set of possible pharmacological small molecule binding models of hERG that is effective in assessing their impacts, as well as the impact of different physiological I Kr models, on the proarrhythmic risk predictions.Determining the details of binding mechanisms, perhaps through the design of an improved calibration protocol, is crucial for mitigating the induced, substantial uncertainty in risk predictions for some compounds.
and is a six-state Markov model with two inactivated closed (IC) states, two closed (C) states, an inactivated (I) state, an open (O) state, and transition rates a 1 to a 7 and b 1 to b 7 Model 3 is a variant of Model 1 where transition between the compound-bound states are also allowed, and happen at the same rate as the (unbound) O ⇌ I transitions.Models 4 to 5i are the trapped equivalents of Models 1 to 2i, where the trapping component is indicated in grey, a 'mirror image' of the hERG physiological model with the same transition rates to admit the possibility of channels closing and preventing unbinding from CD or ICD states (not shown in the schematics) corresponding to C and IC in the drug bound channels.Model 6 relaxes the mirror trapping component of Model 4 by allowing an extra degree of freedom with a trapping rate factor k trap multiplying the original transition rate.Models 7 to 10 allow an extra degree of freedom compared to Models 1, 3, 4, and 6 by assuming independent binding and unbinding rates for open and inactivated states-the modulated receptor hypothesis-whilst enforcing microscopic reversibility by specifying rates indicated by an asterisk as a function of the other rates in the closed loop (Colquhoun et al., 2004).Model 11 introduces independent trapped states for the open and inactivated compound-bound states with transition rates k trap and k untrap , which is inspired by a combination of the 'intermediate encounter complex' model in Windley et al. (2016) and the Li et al. (2017) model.
binding mechanisms and noise processes/artefacts in the data, whilst still distinguishing between plausible and implausible binding mechanisms.A pharmacological binding model was classified as a plausible model if its RMSD was smaller than the maximum RMSD of the bootstrap samples of the data or as plausible as the CiPA v1.0 model if its RMSD was within 1:2Â the RMSD of the reference model.2.4 | AP models, q net , and torsade metric scoresTo assess the impact of the choice of hERG binding model on predictions of drug-induced proarrhythmic risk, we adapted the approach used inLi et al. (2019).We used the optimised CiPA v1.0 AP model byDutta et al. (2017) for predicting the effects on AP due to different hERG binding mechanisms.The 'dynamic hERGbinding model' in the AP model was replaced with one of the successfully calibrated pharmacological binding models in Figure 2. To account for multiple ion channel block effects, the Hill equation characterised by the half-inhibition concentration (IC 50 ) and the Hill coefficient (n)-Model 0α-was used to model the drug effects on three other currents, the fast sodium current (I Na ), the late sodium current (I NaL ) and the L-type calcium current (I CaL ), using their reported median values.The hERG physiological model can be either of the models shown in Figure 1.However, to use the new Figure 1b hERG physiological model in the CiPA v1.0 AP model, we recalibrated its I Kr conductance, in control (drug-free) conditions, by matching the AP duration at 90% repolarisation (APD 90 ) at 0.5-Hz pacing (stimulus amplitude À80 A/F and 0.5-ms duration) at (quasi-)steady state after 1000 paces.We obtained a new I Kr conductance of 0.0912 pA/pF for the hERG physiological model B in the AP model that resulted in the same APD 90 as the CiPA v1.0 model.The hERG-binding-model-replaced-AP models were then used to calculate the q net metric and the proarrhythmic risk prediction.For each model and compound, pacing at 0.5 Hz was initialised from the steady state under control-conditions and continued for 1000 paces after compound addition, at multiples of each compound's maximum therapeutic concentrations (C max ).The q net metric was defined as the net charge over one beat carried by I Kr , I NaL , I CaL , the transient outward potassium current (I to ), the slow rectifier potassium current (I Ks ), and the inwardly rectifying potassium current (I K1 ); this was computed by integrating the sum of the six currents between two consecutive stimuli (with time step 0.01 ms) using the trapezium rule.The proarrhythmic risk prediction was made using the torsade metric score, defined as the mean q net value averaged at 1 Â ,2 Â ,3 Â ,4Â C max .Finally, an ordinal logistic regression model (all-threshold variant,Rennie & Srebro, 2005) was used to estimate the new drug-induced proarrhythmic risk thresholds for low-, intermediate-, and high-risk categories, using the torsade metric as the feature.The classifier was trained with L2 regularisation using only the training compounds and was solved with the L-BFGS-B algorithm in SciPy(Virtanen et al., 2020).The two thresholds for separating (a) the low-risk category from intermediate/high and (b) high from low/intermediate risks were calculated using Threshold example compounds: dofetilide (top), terfenadine (middle), and verapamil (bottom).The percentage current plots in Figure 3 (left) show only the effect of the drug over time for the 10 pulses of holding potential (0 mV) in the CiPA-Milnes protocol.Regardless of the physiological models (A, squares and B, circles), the results of the calibration of the binding models were similar.The results of all the remaining compounds are shown in the Supporting Information.For dofetilide, a trapped drug(Milnes et al., 2010), all binding models except the steady-state conductance scaling model (Model 0α) were able to fit the calibration data, and so we cannot unpick the binding mechanisms of this drug with these data.The four grey horizontal lines in Figure3(left) are predictions of Model 0α, the only implausible model for dofetilide.Model 0α was the only model that showed no dynamics in the percentage current plots as, by definition, it scales only the conductance of the model so its effect on I Kr must be constant over time.Another model that lacked in the drug dynamics was Model 0β, the all-state binding model, which marginally passed the RMSD check.Model 0β modelled the drug effect with the b gate in Equation (2) that is independent of channel state/voltage, producing a single exponential decay over time-even during the long resting potential (À80 mV) of CiPA-F I G U R E 3 Calibration of the binding models and their RMSDs for three example compounds: dofetilide (top row), terfenadine (middle row), and verapamil (bottom row).The left column shows the percentage current of the data (semitransparent line for the mean and highly-transparent area showing AE one standard deviation) and the calibrated binding models (solid, smooth lines) for all four concentrations used during calibration.The right column shows the RMSD of all models compared to the RMSD of the bootstrap samples of the data (box-plot) and the reference model, CiPA v1.0 (red star).Both physiological models A (green squares) and B (orange circles) are shown for comparison.The horizontal red dashed line indicates the maximum RMSD of bootstrapped samples of the data.Grey lines/markers on each panel are the binding models that are ruled out through the RMSD comparison (above the red dashed line)-implausible models.Milnes protocol (although not shown) when the channels were in the closed state(s), resulting in the drop in current between pulses, when the data show (if anything) a slight restoration of current between pulses.

Figure 4
Figure 4 presents a summary of the binding models with physiological model A, which gave plausible fits to the experimental data and/or performed similarly to the CiPA v1.0 model (shown in light/ dark greens and both are considered as plausible models in subsequent analyses); a summary for those with physiological model B are provided in the Supporting Information; the results were very similar with only minor differences for cisapride, ibutilide, and loratadine.
Figure 5 shows inferred binding rate parameters k on , unbinding rates k off , and the Hill coefficients n of the calibrated binding models.All binding models have k on , k off , and n, apart from the conductance scaling model (Model 0α) which has only two parameters IC 50 and n (Table 1).Models 7-10 have independent binding rates for the open and inactivated states, giving two k on shown as empty (for open) and filled (for inactivated) markers; only Model 7 has two (free) k off , shown in the same way, as Models 8-10 have closed-loop states which reduce one degree of freedom due to microscopic reversibility (Figure 2).The CiPA v1.0 model is shown as red stars; the two physiological models, A and B, are shown together as squares and circles, respectively, and the observed results remained similar regardless of the physiological model.The results for all the remaining compounds are shown in the Supporting Information.
ble models shown with dotted lines.The torsade metric decision boundaries for the low-risk (green), intermediate-risk (blue) and highrisk (red) categories in Li et al. (2019) (0.0689C/F and 0.0579C/F, respectively) are indicated as dashed horizontal lines for reference.

F
I G U R E 4 Summary of the selected binding models with physiological model A for all compounds.A binding model (column) is considered to be appropriate for a compound (row)-a plausible model-if coloured in green, where the RMSD of the model to the averaged data is either smaller than the RMSD of the bootstrap samples of data or similar to the CiPA v1.0 model to the averaged data; it also shows the models that are as plausible as the CiPA v1.0 model (darker green) but worse than the bootstrap samples of data.Model 12 is identical to the pharmacological binding component in the CiPA v1.0 model.Compounds are sorted according to the training and validation lists, and their proarrhythmic risks.

Figure 6 (
Figure 6 (bottom row) also shows three more examples: dofetilide, cisapride, and bepridil.These compounds are almost pure I Kr channel blockers at these concentrations.Again, they showed inconsistent spread of q net predictions from the plausible binding models, demonstrating the importance of identifying the correct binding mechanisms or at least narrowing down the possibilities.

Figure 7 (
Figure 7 (left) shows the torsade metric predictions of all binding models with physiological model A for all compounds.The drugs are sorted according to their proarrhythmic risk categories and are split into training and validation lists (Li et al., 2019); the same decision boundaries as Figure 6 are shown as dashed vertical lines, for the

Finally, we compared
the effects of the choice of the I Kr physiological model (Figure1) in predicting the torsade risk classes.The 'O'Hara-Rudy CiPA v1.0' model(Dutta et al., 2017) had hERG physiological model A replaced with physiological model B (solid lines;Lei, Clerx, Beattie, et al., 2019; Lei, Clerx, Gavaghan, et al., 2019), the I Kr maximum conductance was re-calibrated to match the APD 90 of the CiPA v1.0 AP model (dashed lines), as shown in Figure8a(blue lines).The I Kr within the two AP models, shown as orange lines, reveals differences in the dynamics of the two hERG physiological models.Although the total charge carried by the two I Kr models is similar (0.190C/F and 0.144C/F for models A and B, respectively), physiological model B has a smaller current during the early phase of the AP-depolarisation to plateau-and plays a more important role during repolarisation.The two I Kr physiological models show a similar overall transition between the states (Figure 8a): starting from mainly the closed state(s) to the inactivated state(s) during depolarisation/plateau before occupying the open state and back to the closed state(s).For physiological model B, training the ordinal logistic classification model to the torsade metric of the training drugs for all the plausible binding models adjusted the risk category thresholds to be 0.0584C/F and 0.0483C/F for separating the low-risk category from intermediate/high and the high-risk category from intermediate/ low, respectively; the fitted ordinal logistic classifier parameters are provided in Table ), such that the normalised decision boundaries are 0 and 1 for both models; the results for only physiological model B (without normalisation) are shown in Figure S25.Replacing the physiological model of I Kr still produced a similar trend for the drug risk categories, although some of the drugs, such as bepridil and tamoxifen, were clustered closer to the new decision boundaries.The observed differences could result from the new I Kr model, the choice of calibration protocol, and the choice of the metric, all of which were designed for CiPA v1.0.
(Left) Drugs are sorted according to their proarrhythmic risk categories.Dashed vertical lines indicate the CiPA v1.0 decision boundaries for the low-risk (green), intermediate-risk (blue), and high-risk (red) categories.The CiPA v1.0 model is shown as circles, and the conductance scaling model, Model 0α, is shown as squares.(Right) The AP predictions of dofetilide, terfenadine, diltiazem, vandetanib, domperidone, pimozide, and nifedipine at 4Â C max are shown, revealing a range of different behaviours from the binding models, indicated with the same colour code as shown on the left.The plausible models are shown as coloured solid lines, the implausible models as transparent dotted lines, and the CiPA v1.0 model as black dashed lines.The drug-free (control conditions) model is shown as grey dotted lines.In this study, we have designed a set of pharmacological binding models for the hERG channel.After selecting a subset of plausible binding models through calibration to the voltage-clamp electrophysiology data under CiPA-Milnes protocol, we compared their effects, as well as the effects of having different physiological hERG models, on proarrhythmic risk predictions.Our pharmacological hERG binding models accounted for most of the plausible mechanisms by which a compound might bind to the F I G U R E 8 Comparisons of hERG physiological model B in predicting the q net metric and the torsade metric with hERG physiological model A. (a) Comparison of AP, I Kr , and state occupancy using the two hERG physiological models A (dashed lines) and B (solid lines).(b) Normalised q net metric at different C max levels for all binding models with physiological models A (green) and B (orange) for one example of each risk category prediction: dofetilide, terfenadine, and diltiazem.More compounds are shown in the Supporting Information.(c) Normalised torsade metric predictions of all binding models with physiological models A (green) and B (orange) for all compounds.Dashed vertical/horizontal lines indicate the normalised decision boundaries for the low-risk (green), intermediate-risk (blue), and high-risk (red) categories.hERG channel.The choice of these models has covered many of the literature binding models of hERG.For example, our simple binding models (Models 2, 2i, and 7) and trapping models (Models 5, 5i, and 9) are similar to models referred to as 'unstuck' and 'stuck',   respectively, in Gomis-Tena et al. (2020);Escobar et al. (2022).
expect that physiological model B could capture the drug effects better (as the open and inactive state occupancies available for binding should match a hERG1a model's predictions better than a hERG1a/1b model).Yet, since the CiPA v1.0 model and the torsade metric were designed and optimised to work with physiological model A, we might expect it to predict risk better.One could conceivably fit the hERG1a binding data with physiological model B to obtain the binding parameters, then apply these binding parameters to physiological model A to predict changes to native I Kr and therefore clinical risk.However, fewer assumptions would be needed if one directly performed voltage-clamp experiments on hERG1a/1b cell lines (Ríos-Pérez et al., 2021) when calibrating pharmacological binding models, as compounds with differing affinities for hERG1a and 1b have been observed (Abi-Gerges et al., 2011).
. Lei, D. G. Whittaker, and G.R. Mirams conceptualised and designed the study.C. L. Lei and D. G. Whittaker wrote the code.C. L. Lei performed the simulations and data analysis and designed the visualisation and generated the results figures.C. L. Lei and G.R. Mirams wrote the manuscript.All authors revised and approved the final version of the manuscript.
of pharmacological models representing different mechanisms of drug binding, where black (states I and O, and dots) is the physiological model of the hERG channel in Figure 1.Dashed double arrows indicate the rate marked with an asterisk is set by microscopic reversibility.Model 12 is identical to the pharmacological binding component in the CiPA v1.0 model.scaling model which is derived by assuming Model 0β is an instantaneous process, that is, the degree of block b is set immediately to the steady state of Equation (2):