Molecular determinants of pro-arrhythmia proclivity of d- and l-sotalol via a multi-scale modeling pipeline.

Drug isomers may differ in their proarrhythmia risk. An interesting example is the drug sotalol, an antiarrhythmic drug comprising d- and l- enantiomers that both block the hERG cardiac potassium channel and confer differing degrees of proarrhythmic risk. We developed a multi-scale in silico pipeline focusing on hERG channel - drug interactions and used it to probe and predict the mechanisms of pro-arrhythmia risks of the two enantiomers of sotalol. Molecular dynamics (MD) simulations predicted comparable hERG channel binding affinities for d- and l-sotalol, which were validated with electrophysiology experiments. MD derived thermodynamic and kinetic parameters were used to build multi-scale functional computational models of cardiac electrophysiology at the cell and tissue scales. Functional models were used to predict inactivated state binding affinities to recapitulate electrocardiogram (ECG) QT interval prolongation observed in clinical data. Our study demonstrates how modeling and simulation can be applied to predict drug effects from the atom to the rhythm for dl-sotalol and also increased proarrhythmia proclivity of d- vs. l-sotalol when accounting for stereospecific beta-adrenergic receptor blocking.

Introduction

Drug stereoisomers are increasingly viewed as potentially useful means to develop and test pharmacological agents [1,2] due to potentially more potent and/or selective action on their intended protein targets. Stereoisomerism also plays a crucial role in drug toxicity profiles with enantiomers possessing different, on- and off-target effects and capacities for metabolic degradation [3]. This phenomenon is well-described in the drug-induced cardiac toxicity [4] for commonly used drugs including an opioid methadone [5,6], an anti-malarial drug quinine and its diastereomer quinidine [7,8], the local anesthetics bupivacaine and ropivacaine [9-13] as well as for the calcium channel blocker verapamil used in the context of cancer chemotherapy [14,15].

A fascinating and puzzling example of isomerism to impact cardiac safety is the antiarrhythmic sulfonamide drug sotalol. The d-isomer was infamously shown in the SWORD (Survival With ORal D-sotalol) clinical trial [16] to increase mortality and risk of sudden cardiac death in patients, leading to its withdrawal [17]. The racemic mixture comprising d- and l-sotalol, however, is widely used as an effective antiarrhythmic – although not entirely without risk [18-29]. For instance, a prevalence of torsades de pointes (TdP) arrhythmia was estimated to be 0.1–7% depending on therapeutic doses, sex (significantly larger for women) and co-morbidities [28,30-35]. The apparent disparity in cardiac risks of the racemic mixture compared to a pure d-enantiomer could be attributed to differences in interactions between sotalol stereoisomers and cardiac K+ channel proteins, or from the higher efficacy of l-sotalol to act as a beta-blocker [36-40]. Here, we developed a multiscale approach combining cellular and tissue functional kinetic modeling with atomistic simulations and experiments to reveal the specific ion channel interactions with d- and l-sotalol and unveil elusive molecular mechanisms underlying stereospecificity of arrhythmia provoking drug-induced cardiotoxicity.

Cardiotoxicity in the form of abnormal heart rhythm is one of the most critical regulatory concerns for drugs and has resulted in withdrawal of a number of therapeutic agents. Drug-induced arrhythmia has been estimated to affect pharmaceuticals from multiple drug classes [41], with an estimated 3% of all prescription drugs worldwide harboring pro-arrhythmic side effects [42]. Indeed, cardiotoxicity is one of the leading causes of drug attrition [43], and accounts for 22–28% of US post-marketing drug withdrawal [44]. The problem of cardiotoxicity has even plagued drugs intended to treat arrhythmias like sotalol, encainide or flecainide [17,45]. Such deleterious drug side effects have been linked to blockade of the potassium current through the cardiac ion channel protein KV11.1, encoded by the human Ether-à-go-go-Related Gene (hERG) [46-48]. Block of cardiac hERG current (IKr), causes prolongation of the QT interval on the surface electrocardiogram (ECG), sometimes causing potentially deadly arrhythmias [48-50]. While not all hERG block and QT prolongation leads to increased risk of arrhythmia, there is no reliable method to distinguish unsafe hERG blockers from safe drugs [51,52], which hinders the pharmacological assessment of cardiotoxicity and may lead to attrition of safe and effective pharmaceuticals.

A recently proposed Comprehensive in Vitro Proarrhythmia Assay (CiPA) initiative is intended to address this need by improving predictions of drug pro-arrhythmia proclivities through the combination of multi-scale modeling and in vitro experimental assays [53,54]. However, the CiPA platform does not yet provide a ready-to-go recipe on how to predict drug arrhythmogenicity. The example of sotalol discussed in this study, where two stereoisomers have different pro-arrhythmia proclivities, is especially challenging. Atomistic molecular dynamics (MD) simulations of hERG – drug interactions, which we employed here, can be used to help identify stereoselective protein – drug interactions [55,56] and thus molecular determinants of drug-induced arrhythmogenesis. However, since cardiac arrhythmia is an emergent phenomenon, manifesting at the cardiac tissue and whole heart scales [57-60], we set out to create a link between atomistic-scale protein and drug structure-based MD simulations and functional kinetic simulations of cardiac cells and tissue. Recently, we developed such an integrative computational modeling approach that spans scales from the atom to the cardiac rhythm and were able to accurately predict pro-arrhythmia proclivities for dofetilide and moxifloxacin, hERG blockers with different cardiac safety profiles, and directly link them to clinical data [61].

Here, we set out to investigate and reveal the stereospecific molecular mechanisms of d- and l-sotalol causing ventricular arrhythmogenesis by applying a multi-scale modeling pipeline from the atomistic to the tissue scale. We have previously developed and validated atomistic force field parameters for cationic and neutral states of sotalol [62], and validated open conducting atomistic models of the hERG channel [61,63,64] based on a high-resolution cryo-EM structure (PDB ID: 5VA2) [65]. Sotalol is known to access the hERG channel pore via an open state, with subsequent hERG inactivation stabilizing the drug – channel interaction [66]. In this work, we used multi-microsecond long unbiased MD simulations with saturating sotalol concentrations to reveal key binding modalities for cationic and neutral d- and l-sotalol. Then, we used enhanced sampling MD simulations [61,63] to compute state-dependent affinities and “on”/”off” rates for open hERG channel – sotalol interactions, which we used as parameters for our functional kinetic model. We optimized the model for inactivated state block to reconcile and predict emergent pro-arrhythmia risks at cardiac cell and tissue scales [61]. To test model predictions, we obtained d- and l-sotalol through enantiomeric separation and conducted electrophysiological recordings of their hERG channel block. We found the affinities of both enantiomers to hERG channel to be comparable. Our MD simulation informed functional cardiac cell and tissue models predicted sotalol induced QT interval prolongation due to a combined effect of hERG channel block and stereoselective beta-adrenergic receptor (βAR) inhibition to be in good agreement with clinical data.

Methods

A brief overview is provided below. See Appendix A Supplementary Materials for a more detailed description.

Atomistic simulations

We used our atomistic structural model of open hERG channel [61], based on the cryo-EM structure (PDB ID: 5VA2) [65] and generated using Rosetta molecular modeling software [67-69] as well as all-atom force field parameters of cationic and neutral sotalol [62]. Standard CHARMM36 all-atom force fields for protein, lipid, ions [70-72] and TIP3P water model [73] were used as well. Our hERG channel model was embedded in a 1-palmitoyl-2-oleoylphosphatidylcholine (POPC) bilayer and solvated by 0.15 M aqueous KCl solution and sotalol molecules using CHARMM-GUI [74] and equilibrated using a staged protocol by gradually reducing restraints [63] using NAMD program [75]. 8 μs long unbiased MD simulations with 23 sotalol molecules in aqueous solution corresponding to an initial 0.05 M concentration were run on Anton 2 [76]. Steered MD simulations with a sotalol molecule pulled from an intracellular aqueous into the hERG channel pore in five restrained 90 ns runs were used to seed umbrella sampling MD (US-MD) simulations [77]. Positional harmonic restraints for pore domain Cα and selectivity filter (SF) backbone non-H atoms to preserve an open hERG channel conformational state were reduced from 1.0 kcal·mol−1·Å−2 in steered MD runs to 0.2 kcal·mol−1·Å−2 during first 5 ns of US-MD runs and kept for the remainder of those and for the entire Hamiltonian replica exchange US-MD (HREUS-MD) [78] simulations. hERG channel residue sidechains and sotalol conformations were not restrained. US-MD simulations ran for 50 ns (including 10 ns equilibration) per each US-MD window, located at −50 ≤ z ≤ −5.5 Å in 0.5 Å intervals using 10 kcal·mol−1·Å−2 harmonic restraints. z = 0 is the center of mass of hERG SF Cα atoms. HREUS-MD simulations followed US-MD runs after 20 ns and continued for 60 ns per each window. Free energy profiles and diffusion coefficients were computed from those simulations and used to calculate dissociation constants (KD) as well as drug “on” (kon) and “off” (koff) rates [61] used as our functional model parameters.

Functional simulations

We used a previously published approach [61] to develop functional kinetic models of hERG channel – sotalol interactions based on the wild-type drug-free hERG Markov model [61]. MD-derived KD, kon and koff values were used to simulate open-state hERG channel block, whereas for the inactivated state block KDi was assumed to be 4.3-fold less in the model 1 based on ref. [79] and 63-fold less than KDo in the model 2 using a top-down optimization approach and ref. [80]. Rates for the inactivated state block were optimized using experimental dose response curves from this work and ref. [80] for models 1 and 2, respectively, and expression k = kon · KD. These hERG – sotalol models were implemented into the O’Hara-Rudy human ventricular cardiac myocyte model [81] to simulate QT interval in pseudo-ECG recordings from a 1-dimensional strand of O’Hara-Rudy cells paced at 1000 ms basic cycle length (BCL) as was done previously [61]. Beta-adrenergic stimulation by isoproterenol (ISO) was implemented through setting protein kinase A (PKA) phosphorylation effect to maximum, and its inhibition by sotalol, i.e. beta-blockade, was modeled as 1/(1 + [sotalol]/Ki) where [sotalol] is drug’s plasma concentration and Ki is its inhibition constant, which is 38,000 ng/ml for d-sotalol and 650 ng/ml for l-sotalol [38]. Sotalol-induced QT prolongation was compared with clinical data at the matching BCL (RR interval) [82].

Experimental methods

Free sotalol base was prepared from sotalol hydrochloride salt through dissolving it in methanol and KOH solution, separation of supernatant, dissolving it in ethyl acetate, filtering and evaporating. High performance liquid chromatography (HPLC) resolution of the two sotalol enantiomers from free base sotalol was performed with an Agilent 1100 HPLC equipped with a Chiralpak IA 4.6 × 250 mm column, heated at 25 °C. An isocratic solvent system consisting of 20% mobile phase A (ethanol with 0.1% diethyl amine) and 80% mobile phase B (hexane with 0.1% diethyl amine) was used for 15 min for each run. Sotalol was detected by its absorption at 220 nm using an Agilent 1100 diode array detector (DAD). Under these conditions (S) sotalol enantiomer has t = 7.99 min, while (R) sotalol enantiomer has t = 11.88 min. Multiple runs using this method furnished about 25 mg of each enantiomer with almost 100% enantiomeric excess. The separated enantiomers were characterized by 1H NMR and proton decoupled 13C NMR (13C-NMR1H) spectra as well as X-ray crystallography using established methodology [83]. Please see details in the Appendix B Data Supplement and linked spectra files therein.

For electrophysiological recordings we used a human embryonic kidney (HEK) 293 cell line stably expressing hERG maintained in minimum essential medium supplemented with 10% fetal bovine serum and 400 μg/ml G418 as previously described [84]. HEK cells were continuously superfused with HEPES-buffered Tyrode solution. Membrane currents were recorded in the whole-cell configuration established using pipette suction [85]. Leak compensation was not used. Experiments were performed at room temperature (22–25 °C). The data were stored on a computer hard disk and analyzed using PatchMaster and Igor Pro 7 (WaveMetrics, Portland, OR). The IC50 for hERG block by sotalol were measured by fitting the data to Hill’s equation as follows: where f(D) is the current ratio in the presence of the drug at a given concentration (D) and h is the Hill coefficient.

The data and datasets generated and analyzed in the current study can be accessed upon request sent to the corresponding author of the manuscript.

Results

hERG channel binding sites and pathways from sotalol flooding simulations

We previously developed atomistic force field parameters for cationic(+) and neutral(0) forms of sotalol and used MD simulations of their lipid membrane partitioning to validate those models showing water – membrane interface accumulation of sotalol(0) as well as significant free energy barrier and membrane perturbation for sotalol(+) [62]. This indicates ionization state-dependent energetics and kinetics of sotalol – lipid membrane interactions and may affect how the drug gains access to the hERG channel pore.

Here, we explicitly explored this using 8 μs long unbiased MD simulations of our open-state hERG channel model embedded in a POPC lipid bilayer and soaked by an aqueous solution with 0.15 M KCl and an initial 0.05 M concentration of sotalol in water, corresponding to 23 molecules in the simulation box. Utilizing these so-called drug flooding MD simulations, we were able to elucidate potential entry pathways and binding sites for d- and l- sotalol(+) and sotalol(0) interacting with the hERG channel. The distributions for cationic and neutral forms of each enantiomer were notably different (Fig. 1 and Appendix A Figs. S6–S7). d- and l-sotalol(0) (Fig. 1B, left panels) were found to embed into the membrane and interact with hERG channel, whereas d- and l-sotalol(+) (Fig. 1B, right panels) remain predominantly in aqueous solution, sampling much more of the solvent space over the duration of the simulation (see drug aqueous concentration time series and atomic density distribution plots, Appendix A Figs. S6 and S7).

Fig. 1.

Distribution of cationic(+) or neutral(0) d- or l-sotalol around the hERG channel from multi-μs long unbiased MD simulations. (A) Chemical structures of neutral(0) and cationic(+) forms of d- and l-sotalol (B) Snapshots of the molecular systems consisting of the hERG channel embedded in the POPC bilayer, solvated with aqueous 150 mM KCl and initial 50 mM sotalol solution, at the end of 8.1 μs MD simulations. For sotalol molecules within 3.5 Å of hERG protein residues non‑hydrogen atoms are shown in the colored space filling representation, non-interacting sotalol molecules are shown as gray sticks. The hERG channel is shown as green ribbons, POPC lipid tails as thin gray sticks, water as aquamarine surface, K+ and Cl− ions are not shown for clarity. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Dominant binding poses of cationic and neutral d- and l-sotalol are shown in Fig. 2, whereas time series of drug binding z positions in the hERG pore are shown in Appendix A Fig. S8. Interestingly, one d-sotalol (0) or l-sotalol(+) molecule was observed to bind deep into the hERG pore, just below the SF region and interact with the canonical drug binding F656 and Y652 hERG residues [86-90] in the pore lining S6 helices (drug molecule M2 in Figs. 2A and D) for most of 8 μs long MD runs, while another drug molecule was observed to transiently bind below, interacting with F656 and/or Y652 from another domain as well as S660 and other residues at the bottom of S6 helices (M1 in Figs. 2A and D). Interestingly, for l-sotalol(0) and d-sotalol(+) systems, we observed only transient binding of one or two drug molecules at the bottom of hERG channel pore interacting with S660 and other S6 helix residues there (Figs. 2B and C and Appendix A Fig. S8B and C) Hence up to two sotalol molecules were able to bind to the hERG pore, in agreement with our electrophysiological data (as described below). Sotalol flooding MD simulations demonstrated that d- and l-sotalol are likely to enter the hERG channel pore through the aqueous intracellular gate and did not reveal any drug entry through lipid facing fenestrations as in the case of NaV channels [91-94] and also suggested for some hERG binding drugs such as ivabradine [95,96]. Thus, sotalol hERG pore binding pathways closely resemble those for dofetilide [61,63,95,96], which also has methanesulfonanilide moieties and fairly similar polarity [97].

Fig. 2.

Binding sites of neutral(0) or cationic(+) d- or l-sotalol around the hERG channel from 8.1 μs long unbiased MD simulations. (A) d-sotalol(0); (B) l-sotalol(0); (C) d-sotalol(+); (D) l-sotalol(+). Left panels: Time-series rendering for binding of one or two sotalol molecules (labeled M1 and M2) within the hERG pore. Sotalol molecules in the frames are shown by colored sticks from the beginning (red) to the end (blue) of each representative binding event. The hERG channel is shown in the initial (transparent green ribbons) and the final (solid green ribbons) conformations. Canonical drug interacting residues Phe656 and Tyr652 as well as selectivity filter (SF) residues are shown as solid or transparent atom-colored ribbons (C – gray, O – red, N – blue). Right panels: Representative binding poses adopted by sotalol molecules (thick atom-colored sticks with C – cyan, S – yellow, others as above) in the hERG channel pore. Interacting hERG residues (within 3.5 Å of any non-H atoms of the drug) are shown as thick atom-colored sticks (C – gray, others are as above). Non-interacting hERG residues Phe656, Tyr652 as well as its SF residues are shown as thin pink, blue and yellow sticks. Hydrogen atoms are not shown for clarity. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

These qualitative observations do not exclude other possibilities since even very long simulation times of 8 μs might not be sufficient to fully sample other potential drug binding poses and pathways, which can be related to distortions of the hERG pore domain (Appendix A Fig. S9) including partial asymmetric pore closure as evidenced by decrease in interchain distances for S6 segment F656 and/or S660 residues (Appendix A Fig. S10) as well as SF distortions revealed by interchain distances of S624 to G628 SF residues (see Appendix A Fig. S11). These structural changes are similar to those observed for a non-conducting state of the hERG channel during a multi-microsecond long MD simulation under an applied transmembrane voltage [63] and may be indicative of a channel transition to a different conformational state, possibly inactivated. Moreover, sotalol contacts with hERG residues during drug flooding MD simulations were asymmetric and highly variable among MD runs (See Appendix A Figs. S12–S13) due to long-lasting interactions (Appendix A Fig. S14).

Energetics of sotalol – hERG channel binding from enhanced sampling simulations

Our multi-microsecond long drug flooding simulations indicated that enhanced sampling MD runs of a drug entry into the channel pore through an intracellular gate are needed to obtain quantitative information about hERG – sotalol interactions. Therefore, we applied umbrella sampling MD (US-MD) [77], and Hamiltonian replica exchange umbrella sampling MD (HREUS-MD) [78] techniques to compute drug binding affinities and rates. Weak harmonic restraints were applied to hERG pore domain backbone atoms to preserve open hERG conformational state.

First, we investigated binding poses of different sotalol forms, corresponding to free energy minima from US-MD runs as shown in Fig. 3. Both d- and l-sotalol occupy similar positions within the hERG pore (Fig. 3A) and interact with the canonical binding residues F656 and Y652 (see Fig. 3B). Sotalol molecules along with interacting protein side chains, which are not restrained during the simulations, adopt different orientations and conformations (see Fig. 3B and Appendix A Figs. S15–S20) to optimize those interactions. d- and l-sotalol(0) were bound slightly deeper into the pore than their cationic counterparts. Interestingly, all sotalol forms were found to bind just below a ring of Y652 residues, interacting with them, as long as other S6 segment residues from multiple chains (see Fig. 3B) in general agreement with previous studies [98,99]. Use of HREUS-MD methodology allowed us to facilitate sampling of multiple drug orientations and conformations in the hERG channel pore (Appendix A Figs. S16-S20).

Fig. 3.

Representative binding poses of neutral(0) and cationic(+) d- and l-sotalol to open hERG channel obtained from US-MD simulations. (A) hERG channel – bound sotalol structures from four US-MD runs corresponding to free energy minima for each simulation are superimposed and represented by different shades. Two opposite chains of the open-state hERG channel structures are shown as green ribbons. Bound sotalol molecules are shown as thick colored sticks: d-sotalol(0) – blue, d-sotalol(+) – purple, l-sotalol(0) – orange, and l-sotalol(+) – red. hERG SF residues are shown as yellow thin sticks, and canonical binding residues F656 and Y652 as thin pink and ice-blue sticks. (B) Close-up views of sotalol hERG binding poses corresponding to a dotted black box location in panel A. Sotalol molecules are shown as thick atom-colored sticks (C – cyan, N – blue, O – red, S – yellow). hERG channel is shown by transparent green ribbons with residues within 3.5 Å of sotalol non‑hydrogen atoms shown by thin atom-colored sticks (C – gray, N – blue, O – red). Non-interacting SF, F656 and Y652 residues are shown by thin colored sticks as in panel A. Hydrogen atoms are not shown for clarity. Box border coloration in panel B corresponds to coloration of each isoform of d- and l-sotalol in panel A. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

US-MD and HREUS-MD simulations allowed us to estimate the free energies of drug binding, ΔGbind, and to compute dissociation constants (KD) [61,63] of cationic and neutral d- or l-sotalol binding to the hERG channel in the open state (see Table 1, Fig. 4 and Appendix A Figs. S21–S22), which shows good agreement between results obtained with US-MD and HREUS-MD methodologies. Our computed dissociation constants, KD, indicate one to two orders of magnitude stronger binding of neutral d- or l-sotalol to open hERG channel compared to their cationic counterparts (Table 1) in agreement with our previous results for dofetilide and moxifloxacin [61]. Combined KD values (taking into account neutral and cationic drug contributions at physiological pH = 7.4) of about 2 mM are very similar for d- vs. l-sotalol, and also for US-MD and HREUS-MD simulations (Table 1). They are also in excellent agreement with racemic sotalol IC50 values for a non-inactivated hERG S620T mutant [79].

Table 1

Open hERG pore sotalol binding data from US-MS and HREUS-MD simulations. Molecular dynamics (MD) computed drug binding free energies (ΔGbind), dissociation constants (KD) and diffusion coefficients at the drug binding site (Dpore) as well as drug “on” (kon) and “off” (koff) rates used for parameterizing functional scale models. Data are means ± standard error of means (SEM) for 10-ns or 20-ns US-MD and HREUS-MD simulation blocks.

		ΔGbind (kcal mol−1)	KD (mM)	Dpore (10−6 cm2 s−1)	kon (μM−1 s−1)	koff (s−1)
d-Sotalol (0)	US-MD	−5.34 ± 0.74	0.17 ± 0.21	0.13 ± 0.02	7.4 · 102	1.3 · 105
	HREUS-MD	−5.03 ± 0.52	0.28 ± 0.24
d-Sotalol (+)	US-MD	−2.32 ± 0.07	23.3 ± 2.6	0.12 ± 0.02	3.6 · 102	8.3 · 106
	HREUS-MD	−2.44 ± 0.33	19 ± 10
d-Sotalol (pH = 7.4)			2.3 (1.5)
l-Sotalol (0)	US-MD	−4.57 ± 0.51	0.60 ± 0.49	0.16 ± 0.02	7.9 · 102	4.7 · 105
	HREUS-MD	−5.08 ± 0.24	0.26 ± 0.10
l-Sotalol (+)	US-MD	−3.60 ± 0.21	2.9 ± 1.0	0.17 ± 0.05	4.4 · 102	1.3 · 106
	HREUS-MD	−2.80 ± 0.08	10.5 ± 1.4
l-Sotalol (pH = 7.4)			2.0 (2.0)

Fig. 4.

Thermodynamics and kinetics of d- and l-sotalol binding to the open hERG pore from US-MD simulations. Free energy, ΔG, profiles (A) and corresponding diffusion coefficient, D, profiles (C) computed from US-MD simulations, and (B) free energy, ΔG, profiles alternatively computed from HREUS-MD simulations for d-sotalol(0) (blue), d-sotalol(+) (purple), l-sotalol(0) orange, and l-sotalol(+) (red). Error bars are computed from block averaging and represent standard errors of mean. (D) Molecular snapshots of d-sotalol(0) simulated systems with two opposite chains of the open hERG channel shown as green ribbons, selectivity filter residues shown as yellow sticks, and sotalol molecule shown in blue space-filling representation at the pore binding site (opaque) or in the bulk aqueous solution (transparent). The reaction coordinate for these simulations is the z-coordinate with respect to the SF backbone Cα center of mass, leading from the intracellular bulk aqueous solution (z = −50 Å) to the bottom of the selectivity filter (z = −5.5 Å), shown as an arrow and a bounding box with a 20 Å width corresponding to a diameter of the flat-bottom restraint in the xy-plane. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

From hERG channel – sotalol US-MD simulations, we also computed position-dependent diffusion coefficient (D) profiles (see Fig. 4C), which are very similar for all the sotalol forms and indicate a steep ~70-fold drop in the drug diffusion as it enters the confinement of the hERG channel intra-cavitary space, as was also shown previously for dofetilide and moxifloxacin [61,63] as well as sotalol – lipid membrane partitioning [62]. Using ΔG and D profiles, we computed ingress or “on” rates, kon, which are very similar for d- and l-sotalol and somewhat smaller for cationic vs. neutral counterparts (see Table 1). Using computed KD and kon estimates we calculated their products corresponding to “off” or egress rates, koff, similar for d- and l- stereoisomers and 1–2 order of magnitude smaller for sotalol(0) compared to sotalol (+) due to stronger binding of the former (Table 1). kon and koff MD estimates will be used as parameters for functional kinetic models described in Section 3.4.

Experimental validations of stereospecificity in sotalol binding to hERG channel

Sotalol is commercially only available as a racemic mixture of its hydrochloride salt. We separated the racemic mixture into the S- and R-isomers known as d- and l-sotalol, respectively, using quantitative chiral HPLC technique to obtain multi-milligram quantities of both enantiomers (see Appendix A Fig. S1), which were tested using whole-cell patch-clamp measurements with hERG channels stably expressed in HEK-293 cells at room temperature with continuous superfusion of drug-containing external solution. We observed comparable hERG current inhibition by 0.3 mM d- or l-sotalol application, both at depolarized potentials and upon hyperpolarization, indicating classical drug-induced channel blockade (Fig. 5, top panels). Fitting of concentration response curves (Fig. 5, bottom panels) revealed IC50 values of 0.286 ± 0.007 mM and of 0.288 ± 0.010 mM for d- and l-sotalol, respectively. The Hill coefficients from these curves were 1.474 ± 0.051 for d-sotalol and 1.370 ± 0.010 for l-sotalol, respectively, suggesting cooperative interaction between more than one sotalol molecule to block hERG channel in agreement with flooding MD simulation data (Section 3.1 and Fig. 2). IC50 values also did not reveal any distinction in the ability of the two sotalol enantiomers to inhibit hERG current in good agreement with the results of enhanced sampling MD simulations reported above (Table 1). They are more favorable than MD Kd values of 1.5 to 2.3 mM, but the latter were computed for an open state hERG model, whereas in patch-clamp measurements using a voltage protocol in Fig. 5, the channel is mostly in the inactivated state, and sotalol binding is stabilized by channel inactivation [66]. Our experimental and MD data are in agreement with an electrophysiology study of hERG channel interaction with racemic sotalol, where the IC50 for wild-type hERG and S620T hERG mutant with abolished inactivation are 0.516 ± 0.036 mM and 2.22 ± 0.38 mM, respectively [79]. This suggests that our US-MD derived kinetic parameters (Table 1 and Section 3.2) are robust enough to be used for functional scale modeling.

Fig. 5.

d- and l-Sotalol concentration-response data from patch clamp electrophysiology experiments. (Upper panels) Whole-cell currents from HEK-293 cells expressing hERG channels in response to voltage-clamp pulses from −80 mV to +20 mV and then to −60 mV to record pulse-current and tail currents, respectively. Black lines indicate the control currents and blue as well as red lines indicate the decreased currents in the presence of 300 μM of d- or l-sotalol, respectively. (Bottom panel) Concentration-response for block of hERG channels by d- (blue) or l-sotalol (red). Drug effects on tail current peak amplitude are plotted. Data are means ± standard error of means (SEM) from 3 cells each for d- or l-sotalol. Blue and red curves represent the fitted Hill equation.

Functional scale modeling of sotalol – hERG channel interactions and their effect on cardiac rhythm

To develop a functional kinetic model of sotalol – hERG channel interactions, we used the multiscale strategy we recently reported [61]. We used a drug-free hERG Markov model [61] and incorporated open and inactivated hERG channel – drug interactions via corresponding neutral and cationic drug “on” (k and k) and “off” (r and r) transition rates, where x = o or i for open and inactivated state, respectively (see Appendix A Fig. S2). We used US- MD computed drug “on” and “off” rates (Table 1) as our open-state hERG – sotalol(0) and sotalol(+) model estimates (ko_d, kod, ro_d and rod in Appendix A Fig. S2 and Tables S2–S3).

Since we lacked a structural atomistic model of inactivated hERG channel we relied on literature data for relative sotalol affinities for open versus inactivated states with experimental estimates indicating more potent inhibition for the latter [66,79]. We first used the reported 4.3-fold ratio of sotalol IC50 values for a non-inactivated S620T mutant to WT hERG channel from experiments in Chinese hamster ovary (CHO) cells at room temperature [79] as KDo / KDi in our Model 1 (solid light-blue and dotted black curves in Fig.6A) to optimize “on” rates, which provided excellent fit to experimental hERG inhibition data from this study (blue circles in Fig.6A). Notably, however, there is very little hERG inhibition for physiologically relevant drug concentrations from 1.4 to 14.7 μM [34,100] (black box in Fig. 6A): up to ~2% for this model.

Fig. 6.

Concentration dependent block of hERG and QT prolongation by sotalol. (A) Experimentally measured dose dependent inhibition of hERG channel by sotalol (colored symbols) and model 1 optimization based on experiments by us and others from expressed channels in mammalian cell lines for d-sotalol (solid light-blue line) and for l-sotalol (dashed black line). Black box indicates therapeutic plasma concentrations. Experimental data are from: Kramer et. al, 2013 – ref. [117]; Duff et. al., 1995 – ref. [80]; Perrin et. al., 2008 – ref. [79]; and this work – see Fig. 5. (B) Concentration dependent increase in QT intervals by d,l-sotalol with hERG channel block alone (blue circles) compared to clinical data (black diamonds) from ref. [82]. (C) During sympathetic stimulation via concurrent ISO 1 μM application, simulations showed a concentration dependent increase in QT interval by d,l-sotalol dependent hERG block and βAR blockade (blue asterisks) compared to clinical data from ref. [82] (black diamonds). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

To investigate the effect of the hERG channel – sotalol interactions on cardiac electrophysiology, we included hERG – sotalol model 1 (by swapping the drug-free IKr model) into the O’Hara-Rudy human ventricular cardiomyocyte model [81]. We simulated the effect of sotalol on the pseudo-ECGs computed from a 1D strand of O’Hara-Rudy ventricular myocytes and calculated the QT interval [61]. We observed negligible prolongation of the QT interval with model 1 when considering hERG channel block alone by the sotalol racemic mixture (d- and l-sotalol in a 1:1 ratio) as shown in Fig. 6B. In QT prolongation traces shown in Fig. 6C we simulated the sympathetic stimulation of 1D ventricular myocyte strand via the additive impact of 1 μM of potent beta-adrenergic agonist isoproterenol (ISO) and its stereoselective inhibition by sotalol [36-40]. l-Sotalol has a strong beta-blocking activity with Ki = 650 ng/ml vs. 38,000 ng/mol for d-sotalol [38]. Beta adrenergic stimulation did not notably increase the QT prolongation for model 1 (see Appendix A Fig. S4A&C for d- and l-sotalol, respectively). Simulating the stereospecific effect of beta-blockade in addition to hERG inhibition effect by racemic dl-sotalol, model 1 increases (due to l-sotalol, Appendix A Fig. S4C) but still significantly underestimated QT prolongation compared to clinically reported QT data [82]: compare blue and black curves in Fig. 6C.

When the functional model of hERG channel – sotalol interactions was optimized to data from expressed hERG channels in vitro to determine the relative strength of inactivated state block (as shown for model 1 above), we were unable to predict the clinical features of hERG channel block by sotalol. Therefore, we set out to make predictions about what relative affinities of sotalol to the open and inactivated state would allow reproduction of the clinical data. Because the open state affinity is determined by the MD simulations, we performed a “top down” model optimization to predict the inactivated state affinity of sotalol to hERG channel. In other words, we worked backward from the clinically reported effect of sotalol on the QT interval in humans to determine the inactivated state affinity since all the other parameters in the model are constrained by specific measurements or simulations.

The results of the top-down model optimization (model 2) are shown in Fig. 7. The model shown in Fig. 7A for d- and l-sotalol (solid red and dotted black curves, respectively) was identified as best able to reproduce the effect of sotalol on the QT interval. The model predicted that, consistent with some other hERG channel blockers such as dofetilide [79], the relative affinity of sotalol to the hERG inactivated state should be >50-fold the affinity to the hERG open state. We found that a 63-fold open to inactivated state affinity ratio yielded computed IC50 values of 20 μM during a voltage clamp protocol consistent with reference guinea pig ventricular myocyte hERG inhibition data [80] used for fitting this model. This functional kinetic hERG – sotalol interaction model 2 resulted in up to ~40% hERG inhibition for physiologically relevant d- and l-sotalol concentrations from 1.4 to 14.7 μM [34,100] as shown in Fig. 7A.

Fig. 7.

Validation of computational drug screening with human clinical data. (A) Experimentally measured dose dependent inhibition of hERG by sotalol (colored symbols) and model 2 optimization for d-sotalol (solid red line) and for l-sotalol (dashed black line) using data from ref. [80]. See Fig. 6 caption for sources of experimental data. (B) Concentration dependent increase in QT intervals by d,l-sotalol with hERG channel block alone (red circles) compared to clinical data (black diamonds) from ref. [82]. (C) During concurrent ISO 1 μM application to model sympathetic stimulation, model 2 simulations showed a concentration dependent increase in QT interval by d,l-sotalol dependent hERG block and βAR blockade (red asterisks), which compared well to clinical data from ref. [82] (black diamonds).

In Fig. 7B we showed that our predicted model 2 effects of dose-dependent racemic dl-sotalol induced hERG channel inhibition alone in a 1D strand of ventricular myocytes (red curve) underestimates clinical QT prolongation data [82] (black curve). However, after including the combined stereospecific beta-blocking effect of sotalol (dominated by l-sotalol as shown in Appendix A Fig. S4D) in addition to hERG inhibition effects by racemic dl-sotalol, we were able to reproduce clinical QT interval prolongation [82] by model 2 (compare black and red curves in Fig. 7C). Notably, using a 4.3 ratio for hERG open to inactivated state affinities as in model 1 resulted in worse fit to reference hERG inhibition data [80] and smaller QT prolongation compared to clinical data [82] as shown in Appendix A Fig. S5.

Thus, we observed that model 2 accurately predicts clinical QT prolongation by dl-sotalol if we take into account hERG block, beta-adrenergic stimulation and stereospecific beta-blockade. We tested a range of drug concentration effects as a function of time to predict arrhythmogenic impacts of d-, l- and dl-sotalol using model 2 as shown in Fig. 8. We observed that upon beta-adrenergic stimulation, hERG block alone (Fig. 8A–C, top rows) leads to appearance of early afterdepolarizations (EADs) on pseudo-ECGs (left panels for each condition) and variable effects on the cardiac rhythms in time (right panels for each condition indicate the peak positions of R and T waves). Upon beta-blockade pro-arrhythmia markers emerged for d-sotalol at all concentrations similar to the setting of hERG block alone (Fig. 8A, top two rows). For l-sotalol (Fig. 8B, middle rows) and the dl-sotalol racemic mixture (Fig. 8C, bottom rows), a different story emerged: In the setting of hERG block alone, arrhythmogenic patterns emerged for all drug concentrations tested, but when both hERG block and βAR blockade effects were simulated, the computed cardiac rhythms fully normalized. These tests were extended to a wide range of l-sotalol concentrations from 200 to 2400 ng/ml resulting in the same outcome (data not shown). These predictions suggest that a more potent beta-adrenergic blockade by l-sotalol may serve as the key plausible mechanism to reduced arrhythmogenicity with l-sotalol alone and for the racemic mixture.

Fig. 8.

Model prediction of sotalol effects on hERG block and βAR blockade during ISO 1 μM application. The timecourse of pseudo ECGs of model 2 between 820 and 830 s (red bars) at three sotalol plasma concentrations are shown in the left columns of every panel. The peaks of R and T waves are plotted during a 1000 s long simulation timecourse as summary data shown in the right columns of every panel. (A) The predicted additive effect of d-sotalol to block hERG and βAR did not abolish EADs at any tested concentrations. (B) l-sotalol with both hERG and βAR blockade abolished EADs at 1200 ng/ml and higher concentrations. (C) Similar to l-sotalol, racemic d,l-sotalol 1200 ng/ml removed EADs. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Discussion

A crucial link between sotalol stereochemistry and its arrhythmogenicity

Computational modeling and simulation approaches at various length and time scales have demonstrated usefulness to make predictions, suggest new experiments, and to reconcile seemingly disparate existing data [59,60,101,102]. In this study, multi-scale modeling and simulation approaches [57,61,103-105] allowed us to specifically ask questions about the importance of isomerism in determining fundamental mechanisms of arrhythmia proclivity induced by d- sotalol, l-sotalol and the dl-sotalol racemic mixture [36,39,40,106-109]. One benefit of computational modeling and simulation is the ease of performing component dissection to determine which attributes of a system yield emergent impacts [57,103]. Here, we show clearly in MD simulations, which were validated by experimental data, that d-sotalol and l-sotalol have similar affinity and interaction with the hERG channel (see Table 1 and Figs. 3–5), in contrast to observed or proposed stereospecific binding of some other hERG-blocking chiral drugs [4] such as methadone [6], quinidine vs. quinone [8], bupivacaine and ropivacaine [12,13]. Similar stereoisomer hERG channel affinities were previously reported for verapamil [110] and chromanol 293B [111]. Our findings are also in line with a previous experimental observation of similar effects of d-,l- and racemic dl-sotalol on action potential and time-dependent potassium current in isolated cardiac muscle and Purkinje fibers from animal models in the absence of catecholamines [36,112]. Therefore, we conclude that stereo-specificity of drug-induced cardiotoxicity attributed to d-sotalol is unlikely to result from differential binding affinities of sotalol stereoisomers to hERG channel. However, the comprehensive multiscale modeling approach [61] allowed us to go a step further and ask if the stronger beta blocking effect of l-sotalol [36-38,40,112-114] was necessary and sufficient to account for a reduction in pro-arrhythmic risk.

How to resolve disconnect between clinical and experimental data

One of the longstanding conundrums related to hERG block by sotalol is that relatively high concentration of drug required to block the channel when assessed in isolated cloned channels expressed in vitro compared to the relatively low plasma concentrations reported in humans that apparently result in substantial QT interval prolongation [34,82,115]. Human plasma concentrations from patients on sotalol are consistently reported in the 1.4–14.7 μM (500–4000 ng/ml) range [34,100]. Even smaller, 0.5 μM sotalol plasma concentration, results in a 50% likelihood of ≥10 ms increase in human QT interval [115]. This is in contrast to concentrations reported from in vitro measurements that suggest up to 100-fold higher sotalol concentrations required for substantial channel inhibition, although the range of experimentally reported IC50 values is variable [34,79,80,115-123].

We attempted to better understand this disconnect by building computational models that were optimized to experimental cell line data from our laboratory and were similar to previous studies (IC50 values of ~290 μM from this work, 586 ± 179 μM and 4.3-fold inactivated to open state inhibition from [79], 210–450 μM from [118], 103 μM from [115], etc.). Consistent with earlier reports [115,124,125], we also noted that models built on data from expressed channel – drug interactions [79] did not yield clinically observed prolongation of the QT interval [82]. However, when we performed a top-down model optimization, working backward from the clinically reported effects of sotalol on QT interval we found an IC50 value ~20 μM, which reflected ~63-fold higher affinity for inactive to open state. Interestingly, these values are similar to sotalol inhibition data of the reported low- and high-affinity dofetilide binding sites in guinea pig ventricular myocytes [80], although in the study reporting those binding assays, the authors interpret the low affinity site as non-specific binding, not necessarily related to hERG channel current inhibition [80,126]. Additional studies may be needed to provide further clarification.

There is no question that collecting data in the future from a physiologically relevant system with all of the interacting components like the guinea pig ventricular myocyte [80], is likely to yield more reliable data for model optimization. Nevertheless, we were able to use a top-down model approach to constrain model parameters and predict QT interval prolongation as a function of drug concentration that was an excellent agreement with clinically reported data [82] (see Fig. 7C). These results may suggest the critical importance of making drug affinity assessments in a physiologically relevant environment as the differences between the clinical dose and resulting effects cannot currently be reconciled with the dose response curves from expressed channels. Indeed, one possibility is that the induced pluripotent stem cell-derived cardiomyocyte (iPSC-CM) based drug screening technology may create new capacity and potential to do just that [127-129].

Limitations and future directions

In this study we were able to predict sotalol induced QT interval prolongation om the ECG and its stereospecific pro-arrhythmia proclivities with our functional hERG – sotalol models optimized using “on” and “off” drug rates for the open-state hERG channel. Since we lacked a structural atomistic hERG channel model in the inactivated state, we had to use a top-down approach to predict an optimal ratio of sotalol affinities to open and inactivated hERG channels, which turned out to correspond to a ratio of sotalol IC50 values for low and high affinity dofetilide binding sites in ventricular cardiomyocytes [80]. Yet, the low-affinity binding site in the Duff et al. study was interpreted as non-specific binding rather than hERG channel inhibition [80,126]. Such drug binding to other proteins or a lipid membrane itself will alter its effective intracellular concentration, availability for hERG channel inhibition and thus pro-arrhythmia proclivities, which can be investigated in future work. As an alternative, we also tested a ratio of sotalol IC50 values for wild-type and non-inactivated S620T mutant channels from cultured cells [79], with the latter agreeing well with our US-MD K estimates. Functional models optimized using these data provided good fit to our experimental WT hERG inhibition dose response curves but failed to capture dose-dependent clinical QT prolongation [82] in line with previous reports suggesting a significant gaps between drug doses required to cause in vitro hERG channel inhibition and QT prolongation observed in vivo [115].

Such discrepancies can be related to multiple factors not explicitly considered here. For instance, hERG channel isoforms, 1a and 1b, with the same transmembrane but different intracellular domain composition contribute to physiological IKr and have distinct gating kinetics [130], which may modulate drug effect on cardiac electrophysiology. Moreover, mutations in auxiliary hERG binding protein MiRP1 were found to affect drug-induced arrhythmogenicity also due to altered channel gating kinetics [52,131], although WT MiRP1 binding was shown not to affect sotalol hERG channel block [66].

Also, not only hERG channel – drug interaction, but also ventricular myocyte model may have an effect on predicted drug-induced QT prolongation and arrhythmogenicity as was tested and/or discussed previously [132-135]. Here and in previous studies [52,61,134,135] the O’Hara-Rudy (ORd) human ventricular cardiac myocyte model [81] was used providing good agreement with clinical QT prolongation data [82] for our myocyte-based model, whereas comparison with other widely used models such as Grandi – Bers [136], ToR-ORd (Tomek, Rodriguez following ORd) [137], and ten Tusscher – Panfilov [138] models will be performed in subsequent work. We recognize that all models have strengths and weaknesses and that models continue to undergo continuous development.

In this study, MD simulations used an available cryo-EM hERG channel structure [65], which we demonstrated to be in the open conducting state [61,63,64]. However, many hERG blocking drugs including sotalol are suggested to have a higher affinity for the inactivated channel state [48,66,79]. Structural information about potential hERG channel inactivated state just started to emerge from recent cryo-EM [139] and MD simulation [140] studies, but was not yet available for our present work. Interestingly, however, our multi-microsecond long unbiased drug flooding simulations resulted in a hERG channel structure with a distorted SF, pinched in the middle (in 3 out of 4 simulations) and significantly widened at the top, in line with some findings from those new studies [139,140]. A similar SF conformation was also shown to be non-conducting in our previous hERG channel simulations [63] and may potentially represent a transition to the inactivated state of the channel bearing some similarity to C-type inactivated KcsA [141,142] and Shaker family [143,144] channel structures, although there might be multiple structural mechanisms of K+ channel inactivation [144]. Such conformational transition may also explain stable sotalol binding poses, located deep inside the pore in some unrestrained drug flooding MD but not in any restrained US-MD runs, in line with previously suggested sotalol binding to open state and subsequent stabilization by hERG channel inactivation [66].

Finally, we found that beta-blocking properties of sotalol enantiomers are crucial for their pro-arrhythmia proclivities. Yet, in this study we did not explicitly study sotalol binding to beta-adrenergic receptors via MD simulations due to a complexity of this problem, requiring a separate study currently underway in our laboratory using recent human β1 adrenergic receptor (β1AR) structures [145] and state-specific functional models of sotalol – receptor interactions [146], We have also developed a functional model that explicitly couples the effects of the autonomic nervous system (ANS) stimulation on the sinoatrial node (SAN) and ventricular contractile cells and allows for coupled effects of changes in SAN pacing on ventricular myocyte electrophysiology (Clancy et al., manuscript in preparation). Those upcoming studies will allow for a comprehensive exploration of sotalol QT modulation effect at different BCL (RR interval) values [82].

Conclusions

Cardiotoxicity in the form of cardiac arrhythmia is a major reason for multiple drug withdrawals from the market. Such drug-induced arrhythmogenicity has been largely attributed to a blockade of hERG K+ channels, which repolarize cardiac myocytes. However, not all hERG-blocking drugs cause deadly arrhythmias. In the case of the anti-arrhythmic beta-blocker sotalol, a subtle change in drug chemistry due to stereoisomerism leads to different pro-arrhythmia outcomes. In this study, we investigated and predicted molecular determinants of stereospecific d- and l-sotalol pro-arrhythmia proclivities through a combination of atomistic molecular dynamics (MD) simulations with electrophysiological recordings as well as experiment- and MD-informed multi-scale modeling of cardiac electrophysiology at cell and tissue levels. MD simulations predicted that up to 2 sotalol molecules can bind in the hERG channel pore, demonstrated stronger binding of neutral drug forms as well as similar d- and l-sotalol binding poses and affinities for the open channel pore. These predictions are in good agreement with electrophysiology recordings of hERG channel inhibition in HEK cells by both sotalol enantiomers. MD simulation computed “on” and “off” rates for open hERG – sotalol interactions were used as parameters for functional kinetic models, which were optimized to reproduce hERG inhibition dose response curves using cultured cells from this work or cardiomyocyte data from a previous study. We also explored different estimates of sotalol affinities for open and inactivated hERG conformational states using a top-down approach. Optimized models were incorporated into standard ventricular cardiomyocyte cell and tissue simulations to predict dose-dependent d- and l-sotalol effect on heart rhythm. The model based on cardiomyocyte hERG inhibition data was able to predict clinically observed sotalol induced QT prolongation when taking into account both hERG channel and beta-adrenergic receptor inhibition, whereas the model based on a less potent sotalol effect on the hERG channel in cultured cells significantly underestimated clinical QT interval lengthening. The former model also correctly predicted higher pro-arrhythmia risk of d-sotalol compared to l-sotalol and racemic mixture attributing this to stereospecific beta-blocking properties, ameliorating the arrhythmogenic hERG channel inhibition effect.