Introduction

Understanding the dynamic conformational changes of proteins is fundamental to elucidating their functions, interactions, and roles in biological processes. Many proteins, especially membrane proteins that constitute a significant portion of drug targets, exist in multiple functionally distinct states. Capturing these various conformations is crucial for predicting how proteins interact with ligands, designing drugs that selectively target specific states, and uncovering the mechanisms that regulate these interactions. However, experimental techniques like cryo-electron microscopy (cryo-EM) often provide only static snapshots of proteins, typically capturing a single conformational state due to experimental constraints. Computational methods such as molecular dynamics simulations can sample alternative conformations but are limited by timescales and computational resources, often failing to observe meaningful conformational changes that result in functional effects. Enhanced sampling simulation techniques can help extend the timescales, but the biasing factors introduced to accelerate the simulations can sometimes push proteins into non-physiological conformations, potentially skewing the accuracy of the predictions and limiting their biological relevance.

Recent advances in deep learning have revolutionized protein structure prediction, with tools like AlphaFold2 (Jumper et al., 2021a) achieving remarkable success in predicting protein structures based on amino acid sequences. However, conventional applications of these AI-based methods often result in the prediction of a single, static conformation, akin to experimental snapshots. This raises a significant question: Can we harness the capabilities of artificial intelligence to predict different physiologically relevant conformations of proteins, thereby capturing the dynamic spectrum of states essential for understanding protein function and drug interactions?

To address this question, we employed and validated different strategies to guide AlphaFold2 to predict multiple physiologically relevant conformations, surpassing the usual single-state predictions. As a proof of concept, we applied this approach to the human voltage-gated potassium channel K_v11.1, encoded by the KCNH2 or human Ether-à-go-go-Related Gene (hERG) gene, a well-known drug anti-target in pharmacology and cardiology due to its role in drug-induced arrhythmias. hERG is a key player in cardiac electrophysiology, underpinning the rapid component of the delayed rectifier K⁺ current (I_Kr) in cardiac myocytes (Vandenberg et al., 2012). This current plays a crucial role in the repolarization phase of the cardiac action potential (Sanguinetti and Tristani-Firouzi, 2006). Perturbation to hERG channel function, resulting from genetic anomalies or pharmacological interventions, can precipitate multiple arrhythmogenic disorders (Sanguinetti and Tristani-Firouzi, 2006).

The hERG channel is a homotetramer, with each subunit containing six membrane-spanning segments (S1-S6) (Wang and MacKinnon, 2017). The segments S5 and S6 along with the intervening loops and pore helix form the channel pore domain (PD), crucial for potassium ion passage along the central pore, while segments S1-S4 form the voltage sensing domain (VSD), responding to voltage changes across the cell membrane. Notably, the hERG channel also features specialized intracellular regions: the Per-Arnt-Sim (PAS) domain at the N-terminus and the cyclic nucleotide-binding domain (CNBD) at the C-terminus.

The distinctive pharmacological promiscuity of the hERG channel makes it prone to blockade by a diverse array of drugs, creating cardiac safety pharmacology risk in the drug discovery process. Blockade of the hERG channel by drugs can lead to QT interval prolongation known as an acquired long QT syndrome (aLQTS) and escalate the risk of torsades de pointes (TdP), a potentially fatal arrhythmia (Li and Ramos, 2017). This issue has prompted the withdrawal of various drugs from the market and underscored the necessity of incorporating hERG safety evaluations in the drug development pipeline (Ferri et al., 2013; Kocadal et al., 2018; Waldo et al., 1996). The susceptibility for drug blockade is not uniform but varies depending on the channel conformational state, a phenomenon known as state-dependent drug block. Drugs may preferentially bind to and block the channel in specific states (open, closed, or inactivated), which can differentially affect cardiac repolarization and rhythm (Perrin et al., 2008; Priest et al., 2008) and thus confer different risks for aLQTS and TdP arrhythmias as shown in our previous study (Yang et al., 2020).

However, capturing the dynamic spectrum of hERG channel states poses a formidable challenge. While cryo-electron microscopy (cryo-EM) has offered invaluable insights into the putative open state of the channel (Asai et al., 2021; Wang and MacKinnon, 2017), a comprehensive view of the closed and inactivated states has remained elusive. Thus, even as we embark on a scientific era of explosive growth fueled by the convergence of protein structure insights, computational capabilities, and artificial intelligence (AI) based modeling and synthetic data, the next frontier is marked by the need to reveal all relevant conformational states of proteins. The existing knowledge gaps constrain both predictive capabilities regarding drug – protein interactions and the creation of therapies through drug discovery to find specific and selective drugs, or in the case of hERG, to minimize their adverse interactions. For example, our recent study by Yang et al. introduced a multiscale model framework to forecast drug-induced cardiotoxicity at cellular and tissue levels, utilizing atomistic simulations of drug interactions with the hERG channel (Yang et al., 2020). However, the absence of hERG structural models in the inactivated and closed states limited the predictive potential of atomistic scale simulations of state-specific drug binding.

The emergence of AlphaFold2, a protein structure prediction tool driven by machine learning, has brought a paradigm shift in structural biology (Jumper et al., 2021a). AlphaFold2 represents a significant advance over previous methods by using deep learning to predict the three-dimensional structures of proteins based purely on their amino acid sequences (Jumper et al., 2021a). AlphaFold2 primarily requires the amino acid sequence of a protein as its input, but the method utilizes other key elements: in addition to the amino acid sequence, AlphaFold2 can also utilize multiple sequence alignments (MSAs) of similar sequences from different species, templates of related protein structures when available, and/or homologous proteins (Jumper et al., 2021a). Evolutionarily conserved regions over multiple isoforms and species indicated that the sequence is crucial for structural integrity. These additional inputs help the model to identify evolutionary and structural constraints that are crucial for accurate predictions. The output of AlphaFold2 is a predicted 3D structure of the protein that includes inter-residue distance predictions, whereby the model predicts the distances between every pair of amino acid residues in the protein. Predictions about the angles between bonds that connect amino acid residues are also generated as angle predictions that are crucial for determining the precise shape of the protein fold.

AlphaFold2’s limitation, in its default configuration, is that it generates only a single-state structure (Lane, 2023), which for the hERG channel corresponds to the open state. In this study, we introduced an easily replicable and generally applicable approach to guide AlphaFold2 in predicting multiple, physiologically relevant conformations of proteins. By employing multiple structural templates and refining input parameters, we enhanced the predictive capabilities of AlphaFold2, enabling it to generate highly relevant and physically plausible protein conformations beyond the default single-state prediction. We conducted drug docking simulations to predict how specific drugs interact with the hERG channel in different conformational states and performed molecular dynamics simulations to assess ion conduction across these states. Throughout the process, we validated our predictions by comparing them with experimental data, ensuring that both the drug interactions and ion conductive properties aligned with observed experimental outcomes. This method opens new possibilities for in silico studies of protein dynamics, drug design, and safety assessments, allowing researchers to explore the full range of conformational states that proteins may adopt.

Results

Generating hERG channel conformational states

It well known that the hERG channel resides in discrete functional states, minimally comprising closed, open, and inactivated states, which interconnect as a function of time and membrane voltage (Vandenberg et al., 2012). In the open state, the channel conducts K⁺ ions through a central pore. In contrast, the closed and inactivated states are non-conductive due to either a constricted pore at the intracellular gate (closed state) or a distorted selectivity filter (inactivated state) (Vandenberg et al., 2012). So far, published experimental cryo-EM structures resolved the channel in an open state (Asai et al., 2021; Wang and MacKinnon, 2017). Starting with the experimental structure, computational studies explored hERG inactivation by simulating how different membrane voltages can change the selectivity filter and thus affect ion conduction through the channel (Li et al., 2021; Miranda et al., 2020; Yang et al., 2020). These studies are essential but have some limitations, as the high voltages applied can force the channel into unnatural conformations, and the simulations are not long enough to allow observation of state transitions (Shi et al., 2020). To overcome these limitations, we adopted different modeling strategies to guide AlphaFold2 in producing diverse conformations relevant to specific functional states of the hERG channel.

The first modeling strategy involves using a structural fragment from an experimental structure of a homologous protein that exhibits the desired characteristics of the target state we aim to model in our channel. This fragment serves as a structural template, and AlphaFold2 is used to rebuild the rest of the channel while adhering to the constraints of the template. For example, to model the closed state of the hERG channel, it is known that channel closure requires the voltage sensor in the voltage-sensing domain to be in a deactivated conformation. To achieve this, we used the deactivated voltage-sensing domain from the closed-state rat EAG channel cryo-EM structure (PDB 8EP1, residues H208 – H343) (Mandala and MacKinnon, 2022) and combined it with the selectivity filter from the open-state hERG cryo-EM structure (PDB 5VA2, residues I607 – T634) (Wang and MacKinnon, 2017). This hybrid structure was used as the template for AlphaFold2 predictions, as illustrated in Figure 1a. Using these discrete structural fragments, AlphaFold2 was then applied to generate 100 models, specifically configured to encourage diverse prediction outcomes for further analysis.

Figure 1.

For modeling the open hERG channel state, we utilized the existing cryo-EM structure of hERG (PDB 5VA2) (Wang and MacKinnon, 2017) and rebuilt the missing extracellular loops using the Rosetta (Fleishman et al., 2011) with the results shown in Figure 1b. This reconstructed model served as a basis for molecular dynamics and drug docking simulations.

The second modeling strategy addresses situations where structural information about state transitions is either limited or inconsistent. In this approach, we erase regions – expected to undergo changes during state transition – from an existing structure and use AlphaFold2 to sample potential conformations for these regions. The predictions are evaluated using standard accuracy metrics to ensure their physical plausibility. This allows AlphaFold2 to identify possible substates, which are then grouped into clusters of structurally similar models for further analysis. For example, during hERG inactivation, the selectivity filter shifts from an open, conductive to a distorted, non-conductive conformation, as shown by numerous studies on hERG and other K⁺ channels (Butler et al., 2020; Cuello et al., 2010; Fan et al., 1999; Li et al., 2021; Pettini et al., 2023; Schönherr and Heinemann, 1996; Tan et al., 2022; Wu et al., 2023). Moreover, there are a number of studies that do not uniformly suggest a single discrete structure of the inactivated state selectivity filter but propose several alternative conformations (Lau et al., 2024; Li et al., 2021).

To model the inactivated state of hERG, first we configured AlphaFold2 to introduce more uncertainty into the sampling process. As illustrated in Figure 1c, starting with the open state cryo-EM structure (PDB 5VA2) (Wang and MacKinnon, 2017), we removed everything except for the cytosolic domain (S660 – R863), then let AlphaFold2 reconstruct the transmembrane domain. In half of the resulting predictions, including the top-ranked model by prediction confidence, the selectivity filter showed a distinct lateral flip of the backbone carbonyl oxygens at residue V625 compared to the open-state structure. This flip created a potential barrier that could prevent K⁺ ions from crossing between the S3 and S2 ion-binding sites. To further investigate this conformation, we extracted the predicted selectivity filter region (residues Y607 – T634) and merged them with the activated VSDs (W398 – V549) and the cytosolic domain (S660 – R863) from the open-state hERG structure (PDB 5VA2) (Wang and MacKinnon, 2017) to create a new structural template. Then, we generated 100 more models for further analysis.

Clustering of AlphaFold2-generated hERG models reveals predominant substates

A key distinction of our approach is that, rather than relying solely on single-model predictions, we generated a diverse population of models to better explore the conformational landscape. By clustering these models, we identified predominant substates, represented as clusters of structurally similar models. To determine which of these substates are likely to be physiologically relevant, we quantitatively assessed the structural reliability within each cluster using the predicted Local Distance Difference Test (pLDDT). Higher pLDDT scores indicate more reliable and accurate structural predictions (Jumper et al., 2021a). This clustering approach helps to capture a range of conformations that might represent stable states of the protein.

For each predicted conformational state of hERG, we clustered 100 predicted structural models based by their degree of similarity, quantified by the root-mean-square deviation (RMSD), as shown in Figure 2.

Figure 2.

Closed-state clusters

The analysis of closed-state clusters showed only minor differences in RMSD and pLDDT values between them. When comparing the top-ranked models from each cluster, the all-atom RMSD between cluster 1 (the cluster with the highest confidence) and cluster 2 was just 0.36 Å, while the RMSD between cluster 1 and the outlier cluster was 0.95 Å. (Outlier clusters are those with fewer than three members.) This indicates small structural differences among the models. Aside from the outlier cluster, the average pLDDT scores for clusters 1 and 2 were also very similar. The low RMSD values suggest that the predictions are converging on a similar overall conformation, with the minor differences likely due to slight variations in the positioning of the intracellular loop regions. As a result, the top-ranked model from cluster 1 (Figure 2b) was selected for further simulations.

Inactivated-state clusters

As inactivation is known to affect the selectivity filter (SF), we grouped the models by focusing exclusively on similarity of the selectivity filter (S624 – G628) conformations. We ranked the clusters according to the average pLDDT of these specific residues. This method led to the identification of four main clusters and one outlier.

Cluster 1, which has the highest certainty score (pLDDT = 90.3 ± 2.3), consists of models where the majority of the SF carbonyl oxygens point inward toward the central axis, resembling the open-state SF. In contrast, cluster 2 (pLDDT = 85.2 ± 1.6) predominantly shows the carbonyl of SF residue V625 flipped outward from the central axis, along with a noticeable narrowing between the G626 carbonyls. Cluster 3 (pLDDT = 77.5 ± 2.5) is characterized by the reorientation of the G626 carbonyls and occasionally F627. In most models from cluster 2 and beyond, the S6 helix rotates to varying degrees, causing the side chains of the pore-lining, drug-binding residues Y652 and F656 to rotate and extend further into the hERG inner cavity. Both cluster 4 and the outlier cluster have low average pLDDTs (<70), making them unreliable for further analysis.

Interestingly, the inactivated state SF conformations predicted by AlphaFold coincide with proposed hERG C-type inactivation mechanisms as highlighted in other experimental and computational studies (Lau et al., 2024; Li et al., 2021). Specifically, the flipping of V625 and the constriction at G626 carbonyls in cluster 2, was previously reported in a recent study (Lau et al., 2024). Moreover, Li et al.’s computational work revealed an asymmetric SF conformation, where two opposing subunits exhibited similar V625 flipping and G626 narrowing characteristics, while the other two subunits displayed the G626 and F627 carbonyl reorientation characteristic of our cluster 3 (Li et al., 2021). Remarkably, AlphaFold2 was able to independently predict these conformations, despite the fact that they were not part of its training dataset, which had a cutoff year of 2021 (Jumper et al., 2021a) and did not include simulated models. Since cluster 1 SF closely resembles the open-state SF, we picked the best model from cluster 2 for further analysis and simulations, which consistently shows flipping of V625 carbonyls and pinching between G626 carbonyls across all subunits.

Open-state clusters

We combined parts of the presumed open-state cryo-EM hERG structure (PDB 5VA2) (Wang and MacKinnon, 2017), specifically the conductive selectivity filter, the activated VSDs, and the cytosolic domain, as the structural template for AlphaFold2 to test whether it would predict changes to the pore region similar to those observed in other states’ predicted clusters. Post-prediction, all 100 generated models are nearly identical, converging almost uniformly into a single cluster. The highest scoring model closely mirrors the experimental open-state structure (Wang and MacKinnon, 2017), showing no major changes in the pore region.

Comparison of hERG channel states models reveals structural differences in the selectivity filter and channel pore

After further structural refinement in Rosetta (Fleishman et al., 2011; Leman et al., 2020) to resolve steric clashes, the resulting models are compared in Figure 3. In Figure 3a, b, the closed-state model displayed the most constricted channel pore, followed by the inactivated-state and then the open-state model. In the pore-lining S6 helix, the canonical drug-binding residue Y652 (Vandenberg et al., 2012) retains a relatively consistent position with minor variation across all channel state models. The rotation and shift of the S6 helix in the inactivated and closed states affect the position of another canonical drug-binding residue, F656 (Vandenberg et al., 2012). The adjustment caused the F656 side chain to extend more into the hERG inner cavity in both the closed and inactivated states, compared to the open state.

Figure 3.

Shown in Figure S1a, b, the SFs of the open- and closed-state models display similar conformations, with carbonyl oxygens along the ion path all oriented toward the central axis as in other K⁺ channel structures, e.g., KcsA and K_v1.2, enabling efficient knock-on K⁺ conduction (Doyle et al., 1998; Long et al., 2005). In contrast, the inactivated-state model SF is distinct, marked by the lateral rotation of the V625 backbone carbonyls away from the central axis (Figure S1c), thereby creating a potential barrier preventing ion crossing from the S3 to S2 binding site. Additionally, we noted a constriction between the G626 backbone carbonyls and a repositioning of the S624 sidechain hydroxyl oxygens. In the model representing the inactivated state, the carbonyl oxygens of G628 and F627 exhibit an upward shift relative to their positions in the open-state model. Figure S1d-f presents the SF of all three models from an extracellular perspective. In both the closed- and inactivated-state models, the F627 side chain undergoes a clockwise rotation when contrasted with its orientation in the open state. This rotational behavior aligns with findings from prior simulation study (Miranda et al., 2020) where it was noted in a metastable non-conductive state. The loop that links the upper SF to the S6 helix rotates anti-clockwise relative to its position in the open-state model, consequently narrowing the upper part of the SF.

For the VSD, we measured the distances between the backbone C_α atoms of the gating charge residues (K525, R528, R531, R534, R537, K538) on the S4 helix and the gating charge transfer center residue F463 on the S2 helix, as shown in Figure 3c. Although we observed an increased distance between the gating charge residues and the charge transfer center residue in the closed state, this separation was not due to a straight downward movement of the S4 helix. Instead, the closed-state model S4 exhibited a minor kink around residue R531 and lateral movement toward the channel center, impacting the S4 – S5 linker and consequently nudging the S5 helix inward, effectively narrowing the pore. The predicted closed-state model exhibits lower confidence levels for the S4 helix and S4-S5 linker residues (pLDDT ≤ 75) when compared to their counterparts in models of other states, necessitating caution in interpreting the physiological implications of this observation. Conversely, the pore region, which demonstrates closure, is characterized by a higher prediction confidence (pLDDT ≥ 75), suggesting a more robust and reliable structural representation. Supporting Movie 1 shows an animation of state changes of the hERG channel models.

We aimed to further investigate the molecular interactions that contributed to channel inactivation through modulation of the selectivity filter conformation. In Figure S2, we analyzed interactions of extracellular S5-P linker (I583 – Q592) along with SF and surrounding SF residues (S620 – N633, Figure S2a) through heatmaps detailing hydrogen bonding, π stacking, cation-π, and salt bridge formation similarities and differences (Table 1 shows detection criteria; distance-based contact maps (Noel et al., 2012) for all residues are shown in Figure S3). Distinct interaction patterns between open- and inactivated-state models were observed in these regions (Figure S2b, c). In the open-state model, N633 atop the SF forms hydrogen bonds with S5-P linker G584 and Q592 from an adjacent subunit, while N629 forms hydrogen bonds with an adjacent-subunit S631. Additionally, SF G626 forms an intra-subunit hydrogen bond with S620 behind the SF. Within the S5-P linker, N588 and D591 also display hydrogen bonding. However, these stabilizing interactions in the open-state model SF region are absent in the inactivated-state model, where only intra-subunit hydrogen bonds between I583 (S5-P linker) and N633 (SF) occur, along with V630 hydrogen-bonding with the same-subunit N629 atop the SF. To corroborate our findings, mutations involving the residues discussed above have been shown to impact hERG inactivation as evidenced in numerous clinical and experimental studies (Butler et al., 2018; Clarke et al., 2006; Cordeiro et al., 2005; Dun et al., 1999; Fan et al., 1999; Ficker et al., 1998; Miranda et al., 2020; Nakajima et al., 1998; Satler et al., 1998) (see Table 2 for more details).

Table 1.

Table 2.

There are fewer differences between the open- and closed-state models (Figure S2b, d). In the open-state model, D591 from the S5-P linker forms intra-subunit hydrogen bond with N588, and G584 hydrogen-bonds with N633 at the top of the SF. These interactions are absent in the closed-state model, where H587 (instead of N588) hydrogen-bonds with D591 within the S5-P linker, and I583 (replacing G584) interacts with N633 at the SF top. Additionally, G628 from an adjacent subunit forms a hydrogen bond with N629 atop the SF. Analyzing the differences between the inactivated- and closed-state model (Figure S2e), the inactivated-state model uniquely features an intra-subunit V630-N629 hydrogen bond, whereas the closed-state model exhibits intersubunit hydrogen bonds between N633 and Q592, and between S631/G628 and N629. Furthermore, in the closed-state model, S620 forms an intra-subunit hydrogen bond with G626, stabilizing the SF conformation.

In Figure S4, we compared the S6 pore-lining helix orientation across various models. The closed-state model features a straight S6 helix with a minor kink. On the contrary, both the open- and inactivated-state models exhibit a pronounced kink around I655, as identified in a prior study (Thouta et al., 2014), which facilitates pore opening and distinguishes the inactivated-state from the closed-state model. Notably, a slight rotation differentiates the S6 helix in the open- and inactivated-state models, altering the conformation of drug-binding residues Y652 and to a greater extent, F656. The interaction network analysis results from Figure S2b, c suggest that alterations in the hydrogen bond network around the SF region, during the transition from open to inactivated state, might pull on the S6 helix and influence its orientation (Figure S1e, f) – a subtle yet potentially impactful change for drug binding. In agreement with our observations, a study by Helliwell et al. (Helliwell et al., 2018a) also suggested that a slight clockwise rotation of the S6 helix in the hERG open-state cryo-EM structure (Wang and MacKinnon, 2017) could align the S6 aromatic side chains, particularly F656, into a configuration enabling interactions with inactivation-dependent blockers that more accurately reflects experimental data (Helliwell et al., 2018b).

Molecular dynamics simulations show K⁺ ion conduction in the open-state model but not in the inactivated state

We performed all-atom molecular dynamics (MD) simulations on two hERG channel models described above, one in the open state and the other in the predicted inactivated state, to evaluate their ion conduction capabilities. Unlike the closed-state model, both open- and inactivated-state models should allow ions and water to enter and traverse the channel pore reaching the SF region. However, only the open state-model is expected to facilitate ion conduction through its SF.

To investigate ion conduction in the selectivity filter, we considered two conditions, as shown in Figure S5a: one in which the SF initially contained only K⁺ ions and another in which both ions and water were present to test previously proposed direct (or Coulombic) and water-mediated K⁺ conduction knock-on mechanisms (Lam and de Groot, 2023; Roux, 2017) as in a previous study (Miranda et al., 2020). In the direct knock-on (ions-only) scenario, we manually positioned K⁺ atoms in the putative K⁺ binding sites of S0, S2, S3, and S4 within the SF. For water-mediated knock-on (the alternating ions and water molecules) scenario, K⁺ ions were placed in the S1, S3, and S_cav positions, while water molecules were inserted into the S0, S2, and S4 positions. These models were incorporated into phospholipid bilayers consisting of 1-palmitoyl-2-oleoylphosphatidylcholine (POPC) molecules and hydrated by 0.30 M KCl, as depicted in Figure S5b. Subsequently, we conducted MD simulations for each case under three membrane voltage conditions: 0 mV, 500 mV, and 750 mV, each lasting 1 μs. This resulted in a total of six MD simulations for each model.

In all instances where a non-zero membrane voltage was applied after equilibration, we observed K⁺ conduction for the open-state model (Figure 4a, c and Figure S6a, c), whereas such conduction was not observed for the inactivated-state model (Figure 4b, d and Figure S6b, d). For ions-only initial SF arrangement we observed that all K⁺ ions initially located in the SF went across during 1 μs MD runs under applied 750 and 500 mV membrane voltages (Figure 4a and Figure S6a), whereas for the alternating water-ion initial SF configuration we observed conduction of SF ions as well as additional K⁺ ions moving all the way across the channel pore (Figure 4c and Figure S6c). In both cases we saw combination of direct and water-mediated knock-on mechanisms as in our previous hERG channel MD simulations (Miranda et al., 2020; Yang et al., 2020). Control MD simulations conducted under zero voltage conditions revealed a single K⁺ SF conduction event for the open-state model (Figure S6g) when the SF was initially filled with water molecules and ions, while no conduction events were observed in the remaining cases (Figure S6e, f, h).

Figure 4.

Subsequently, we conducted an analysis of pore radius changes throughout the MD simulations (see full results in Figure S7). In under zero voltage conditions, we observed consistent and distinct pore radius profiles across all simulations within their respective models (left panel in Figure S8a). Specifically, MD simulations featuring the inactivated-state model consistently displayed a narrower pore radius when compared to simulations involving the open-state model. However, when subjected to high voltage conditions, the open-state model exhibited a shift towards an inactivated-like state, leading to a reduction in the pore width (right panel in Figure S8a), which is consistent with an increased hERG channel inactivation propensity at more depolarized voltages (Vandenberg et al., 2012). Although we did not observe the outward flipping of the V625 backbone carbonyl oxygens in the SF during the 1 μs simulations of the open-state model, we did observe the flipping of the F627 backbone carbonyl oxygens as shown in Figure S8b. Interestingly, this SF conformation with flipped F627 but not V625 is also observed in cluster 3 of the AlphaFold-predicted models shown in Figure 2. In the inactivated-state model MD simulations, increasing membrane voltage resulted in an increased probability of flipping of the V625 backbone carbonyl oxygens into ion conducting conformation. Notably, a single outward oriented carbonyl oxygen still prevented K⁺ conduction through the SF.

Comparison with previously reported K⁺ channel C-type inactivation mechanisms

Cuello et al. in their study of KcsA channel identified a similar constriction at G77 within the SF and a corresponding reorientation of the V76 carbonyl, resulting in a dilation in the SF at this location and corresponding loss of the S2 and S3 ion binding sites (Cuello et al., 2010). They suggested this backbone rearrangement as a fundamental molecular mechanism underlying C-type inactivation in K⁺ channels (Cuello et al., 2010). In other studies on Shaker and K_v1.3 channels, dilation in the upper SF that disrupts the S1 and S2 K⁺ binding sites has been proposed to be a potential C-type inactivation mechanism (Chandy et al., 2023; Selvakumar et al., 2022; Tan et al., 2022; Tyagi et al., 2022). Similar dilations in the SF are also predicted by AlphaFold2, particularly within cluster 3 of the predicted inactivated state hERG channel clusters shown in Figure 2a. Although these models were not simulated under our study, such dilated conformations of the SF also emerged during our MD simulations of the open-state model under applied voltage.

We further analyzed predicted inactivated-state conformations by plotting the cross-subunit distances between the carbonyl oxygen atoms of SF residues at 750 mV and 500 mV applied voltages, as shown in Figure S9. The dilation observed in the hERG channel, which also occurs in the upper SF, differs from that in the aforementioned K⁺ channels. In Shaker-family channels, the most considerable widening occurs at the SF tyrosine residue (Y445 in Shaker/Y377 in K_v1.2) immediately below the topmost SF residue (G446 in Shaker/G378 in K_v1.2) (Tan et al., 2022; Wu et al., 2023). Conversely, in the hERG channel, the topmost SF residue (G628) exhibits the most significant widening, followed by the residue immediately below it (F627). Our MD simulations of the hERG channel reveal that its dilation process involves two sequential steps: SF near residues F627 dilates first, followed by SF near topmost G628 residues. The latter step occurs faster at higher voltages (750 mV) compared to lower voltages (500 mV). We present these steps in Figure S10. Notably, despite the dilation of the hERG SF, ion conduction is still observed across all replicas, in contrast to the Shaker channel (Chandy et al., 2023; Selvakumar et al., 2022; Tan et al., 2022; Tyagi et al., 2022).

Computational drug docking reveals state-specific differences in drug binding

We utilized Rosetta GALigandDock software (Park et al., 2021) to dock 19 drugs from different classes, considering their multiple protonation states, into our hERG state-specific channel models. This process aimed to evaluate and corroborate state-dependent binding interactions with experimental studies, specifically in terms of relative binding affinities. Figure S11 presents these findings in the form of Rosetta GALigandDock (Park et al., 2021) binding energies (lower values mean more favorable binding). Consistent with published studies, most drugs showed stronger binding to the inactivated-state hERG channel model, including astemizole, terfenadine, cisapride, d/l-sotalol, dofetilide (Perrin et al., 2008), and E-4031 (Numaguchi et al., 2000; Wang et al., 1997). Drugs like moxifloxacin (Alexandrou et al., 2006), quinidine (Perrin et al., 2008), and perhexiline (Perrin et al., 2008) did not show strong preference for the inactivated-state model, aligning with findings from hERG experimental studies using inactivation-deficient mutants (Perrin et al., 2008) or “step-ramp” voltage protocol (Alexandrou et al., 2006).

In our GALigandDock docking results, most drugs exhibited increased binding affinity to the closed-state hERG model compared to the open-state hERG model. Drugs are unable to bind to the closed state from the intracellular space because the pore is closed. However, they can become trapped if they are already bound when the channel transitions from open to closed state, as shown in experiments for dofetilide (Windley et al., 2017), cisapride (Windley et al., 2017), terfenadine (Windisch et al., 2011; Windley et al., 2017), E-4031 (Windisch et al., 2011), and nifekalant (Kamiya et al., 2006).

To model drug trapping, we placed the drug in a pocket beneath the selectivity filter in the closed pore configuration before docking. However, this method does not consider how the conformational shift from the open to the closed state might influence drug binding. Under physiological conditions, the pore gating motion from open to closed might expel drugs from the pore instead of pushing them deeper. This limitation might account for some inconsistencies noted in our docking study, particularly regarding the apparent trapping of drugs such as amiodarone and haloperidol, which is at odds with experimental results (Stork et al., 2007). However, these preliminary results could pave the way for more thorough investigations, employing advanced computational techniques to delve deeper into the dynamics of drug trapping (Branduardi and Faraldo-Gómez, 2013; Miao et al., 2020).

As an example, Figure 5 focuses closely on the binding of astemizole, dofetilide, and quinidine in their cationic forms (the predominant species at the physiological pH 7.4) to different hERG channel models. Experimental evidence (Perrin et al., 2008) suggests that dofetilide and astemizole bind more potently to the inactivated state, while quinidine does not exhibit this preference.

Figure 5.

Astemizole

In the open-state model, astemizole engages in π stacking with the Y652 aromatic side chain via its benzimidazole group and forms hydrophobic contacts with adjacent residues (A653 and F656) and Y652 from another subunit. Its binding occurs just below S6 helix residue Y652, not deeply within the pore, resulting in a binding pocket that remains partially open to external solvent. In contrast, in the inactivated-state model, astemizole demonstrates an increase in interactions and becomes fully enclosed within the binding site deeper into the pore, engaging in hydrophobic interactions with S6 residues S649, Y652, F656, and S5 residue F557. The astemizole methoxyphenyl group engages in the π stacking with hERG Y652 and F656 side chains from adjacent subunits. At the other end, the aromatic ring of Y652 hydroxyphenyl side chain also engages in π stacking, while its oxygen atom acts as a weak hydrogen bond acceptor (i.e., donor angle of ≤ 45° and bonding distance of ≤ 3.8 Å) to astemizole secondary amine group. Furthermore, the SF residue S624 side chain participates as a hydrogen bond acceptor with the cationic ammonium group of astemizole piperidine ring. The binding pattern in the closed-state model is similar to that in the inactivated-state model, but it includes additional hydrophobic interaction with G648 and I655.

In 2021, a study by Asai et al. solved the cryo-EM hERG structure in complex with astemizole (Asai et al., 2021). Although the structure deposited to the PDB (7CN1) (Asai et al., 2021) did not include ligand coordinates, a possible astemizole binding pose was suggested based on residual EM density. The density map showed a π-π interaction of astemizole with Y652 and a hydrogen bond with S624 (Asai et al., 2021). In addition, electrophysiology study by Saxena et al. has identified F557 and Y652 as crucial for astemizole binding, as determined through mutagenesis (Saxena et al., 2016).

Dofetilide

In the open-state hERG model, dofetilide interacts with multiple residues of the pore lining S6 helices. For instance, it forms hydrogen bonds with two S660 residues from opposite subunits and two A653 residues from adjacent subunits. Additionally, there are hydrophobic interactions with Y652 and G657. When bound to the inactivated-state hERG model, dofetilide settles deeper into the pore, extending horizontally above the Y652 residue, where it forms π stacking and cation-π interactions, as well as hydrogen bonds with the Y652 residues of all subunits. The backbone carbonyl oxygen of residue S649 also engages in hydrogen bonding with the nitrogen in the sulfonamide group of dofetilide. There are also hydrophobic contacts with hERG S6 residues S649, F656, and I655 as well as residue F557 on the S5 helix. The closed-state model pore provides a substantially narrower space, leading dofetilide to bind vertically along the pore axis, forming hydrophobic contacts with S624, S649, Y652, F656, and S660. Here, the oxygens and nitrogen of its sulfonamide group create a series of hydrogen bonds with hERG residues T623 and Y652. The central cationic tertiary ammonium group of the ligand forms cation-π interactions with two F656 residues from adjacent subunits. Lastly, the oxygens in the remaining sulfonamide group are involved in hydrogen bonding with the S660 residues from 3 subunits. In the inactivated-state and closed-state models, the binding pocket is completely enclosed, unlike the open-state model.

In line with our results, mutations of A653 affected hERG channel blocking by dofetilide (Stepanovic et al., 2009). Additionally, electrophysiology study by Saxena et al. revealed that mutants Y652A and F557A altered dofetilide-induced hERG inhibition (Lees-Miller et al., 2000; Saxena et al., 2016). Another experimental mutagenesis study by Kamiya et al. have highlighted possible dofetilide interactions with hERG S6 residues G648, Y652, F656, and V659 as well as its pore helix residues T623, S624, and V625 (Kamiya et al., 2006) in general agreement with our study.

Quinidine

Unlike the previous two drugs, quinidine does not exhibit significant state-specific binding based on our docking predictions (Figure S11). In the open state hERG model, its small size allows for complete enclosure within the binding pocket, where it forms distorted hydrogen bond with residue Y652, hydrophobic contacts with another Y652 from an adjacent chain, as well as other S6 residues A653 and S660. In the inactivated-state model, quinidine binds lower in the pore, resulting in a binding pocket partially open to solvent. Here, its quinoline rings engage in π-stacking interactions with the F656 side chain above and hydrophobic contacts with another F656 from the opposite subunit. The quinidine hydroxyl oxygen forms a weak hydrogen bond with residue S660. In the closed-state model, quinidine becomes trapped deeper in the pore and is fully enclosed. Its quinoline rings π stack with multiple F656 side chains below and the Y652 side chain above drug binding position. The quinidine hydroxyl oxygen now establishes a weak hydrogen bond with Y652, with additional hydrophobic contacts observed with other Y652 and F656 side chains and the T623 backbone. Corroborating our results, Sănchez-Chapula et al. previously found that mutating key aromatic residues (Y652, F656) in the S6 segment reduced the efficacy of quinidine by more than 100-fold (Sănchez-Chapula et al., 2003).

Validation of state-dependent drug block with experimental data using hERG Markov model

There are several complications that make it difficult to directly compare experimental binding affinities with predicted affinities from simulations. During electrophysiological recordings of hERG inhibition by various drugs, hERG has been shown to adopt various functional states, presumably corresponding to different protein conformation states. These states can be bound by drugs with varying affinities, with drug ionization state also being a contributing factor. Additionally, the variability in experimental protocols affects the measured affinities (Gomis-Tena et al., 2020). In general, electrophysiological measurements report the IC₅₀, the drug concentration required for 50% inhibition of current. However, the IC₅₀ value is not directly comparable to computed affinities from drug docking.

To address these challenges, we developed a novel computational approach that combines modeling and simulation to predict hERG channel conformational state probabilities (open, closed, inactivated) over time. First, we collected a comprehensive set of experimental data and employed a hERG functional model with five functional states, which was extensively validated in our earlier study (Romero et al., 2015). For each drug, we ran in silico electrophysiological experiments under the same conditions as the experimental studies, allowing us to calculate the relative probabilities of the various hERG channel states specific to the drug and protocol. These state probabilities were then used to refine the computed binding affinities from docking simulations. We adjusted the affinities for both neutral and charged forms of each drug according to their prevalence in each conformational state. This method allowed us to scale the binding predictions based on the likelihood of each channel state occurring during the experimental protocols. Finally, we compared the simulated binding affinities with experimental hERG binding potencies (Table 3), offering a new validation technique that enhances the accuracy of our predictions and helps reconcile the differences between experimental IC₅₀ measurements and computed affinities.

Table 3.

We compared the experimental drug binding potencies with the simulated binding affinities, starting with the traditional approach of using only open-state drug binding simulations (Figure 6a, commonly used in ion channel pharmacology due to the absence of multi-state models) and then expanding to simulations that account for drug binding to different states (Figure 6b). The comparison revealed that while the conventional open-state drug binding model showed a moderate correlation with experimental drug potencies (R² = 0.45, Pearson’s r = 0.67, p-value = 0.00046), including multi-state binding improved how well the model explained the experimental data by 35.6% (R² increased to 0.61) and strengthened the correlation (Pearson’s r increased to 0.78, p-value = 0.00001).

Figure 6.

This improvement highlights the importance of including multiple conformational states in drug docking simulations, offering a more complete view of how the drug interacts with the ion channel. It also further supports potential physiological relevance of our predicted inactivated and closed state hERG channel models. Notably, this stronger correlation was achieved without considering how drugs may influence ion channel gating dynamics — an aspect that could further refine our predictions and will be addressed in future studies. Despite these limitations, the results suggest that multi-state drug binding simulations may provide a more reliable approach than conventional single-state ion channel models, showing promise for more accurate drug potency predictions in ion channel pharmacology.

Discussion

In this study, we introduced a methodology to extend AlphaFold2’s predictive capabilities by guiding it to sample multiple protein conformations reminiscent of different ion channel states, using the hERG potassium channel as a proof of concept. While other studies have focused on generating diverse conformations, we aimed to take this further by validating the physiological relevance of these predicted states through computational simulations and experimental data. By incorporating multiple structural templates and refined input parameters, we directed AlphaFold2 to predict distinct functional states, including the closed and inactivated forms of hERG. This validation is particularly important, as understanding the full conformational spectrum of a protein can provide valuable insights for pharmacological development, especially for an anti-target like hERG, which exhibits state-specific drug interactions (Li and Ramos, 2017; Perrin et al., 2008; Sanguinetti and Tristani-Firouzi, 2006; Vandenberg et al., 2012; Yang et al., 2020). Such knowledge can aid in designing safer drugs that avoid undesirable interactions, such as drug-induced cardiotoxicity, by selectively targeting or avoiding particular hERG channel states. This proof of concept illustrates the potential of our methodology for hERG and offers a framework that could be extended to other membrane proteins, contributing to advancements in in silico drug safety assessments and development.

The first strategy involved using AlphaFold2 as a generative tool, where we provided structural fragments indicative of a specific state of the channel as templates and instructed AlphaFold2 to reconstruct the remaining structure from these fragments. This approach was particularly effective in modeling the hERG channel in its closed state with a non-conductive constricted pore. The second strategy focused on using AlphaFold2 to investigate a wide range of potential conformational states, especially in scenarios where there is limited knowledge. Remarkably, AlphaFold2 managed to generate conformations that had been previously documented in the literature (Lau et al., 2024; Li et al., 2021) but not included in the training dataset, showcasing its broad predictive capacity. We further validated the physiological significance of these states using computer simulations and comparing them to experimental results as discussed below. This methodology has provided new insights into the structural and functional behavior of the hERG channel, including state-specific drug interaction mechanisms.

The predicted models of the closed, open, and inactivated states of the hERG channel, have revealed detailed structural aspects that were previously elusive, particularly in the pore region where drug binding occurs. The closed-state hERG channel model demonstrates constricted pore and deactivated voltage-sensing domains, suggesting a structural basis for channel deactivation. Conversely, the open state hERG channel model indicates a pore conformation that allows ion conduction and aligns with the physiological role of the channel in cardiac repolarization. The inactivated-state hERG channel model suggests a non-conductive selectivity filter conformation with notable changes to the pore region that underpin new insights into state-specific drug-induced channel blockade. The structural models similar to ones predicted in our study have not been published in any publicly available databases and were not used for training AlphaFold2, to our knowledge. Our approach thus underscores the substantial yet untapped potential of AlphaFold2 to predict de novo channel conformations with physiological relevance.

One of the most compelling predictions here is that AlphaFold2 predicted models that directly corresponded to the previously reported “inactivated state” conformations from prior experimental and simulation studies. These structures were identifiable in discrete protein clusters (Figure 2a) with high prediction confidence metrics and likely together constitute components of the hERG inactivation mechanism. In this way, AlphaFold2 provided a strong method to reconcile apparently disparate previous experimental and computational data.

Predicted by AlphaFold2 and validated in our study, the non-conductive SF conformation is characterized by the flipping of SF V625 carbonyl oxygens away from the central axis and minor constriction at the level of G626. Echoing our observations, Lau et al. suggested that such flipping may constitute a C-type inactivation mechanism for the hERG channel (Lau et al., 2024). Their study involved altering the external K⁺ concentration to predispose the channel towards a non-conductive, presumed inactivated state (Lau et al., 2024). Li et al. made a similar observation in their computational study, demonstrating the constriction between G626 carbonyls which occurs asymmetrically in only two opposing subunits and not in the others to be a potential C-type inactivation mechanism (Li et al., 2021). Indeed, this constriction was accompanied by the flipping of V625 carbonyls in the same subunits where the constriction was observed (Li et al., 2021). Although these studies have proposed potential molecular mechanisms for hERG C-type inactivation (Lau et al., 2024; Li et al., 2021), they have primarily concentrated on conformational shifts leading to a non-conductive selectivity filter. Our study takes a step further and sheds light on how these alterations extend to the pore region and subsequently impact drug binding as seen in experiments, which is an outstanding issue in safety pharmacology and drug development as hERG channel is a notorious drug anti-target (Perrin et al., 2008; Yang et al., 2020).

For other K⁺ channels, dilation in the SF has been proposed to be potential C-type inactivation mechanisms (Chandy et al., 2023; Cuello et al., 2010; Selvakumar et al., 2022; Tan et al., 2022; Tyagi et al., 2022). While inactivation processes across various K⁺ channels may share some similarities, the associated conformational changes can adopt distinct differences due to small variations in the SF sequences, which could explain the observed variability in inactivation rates among hERG and other ion channels (Vandenberg et al., 2012). As an example, a study on Shaker-family channels suggested that a two-step widening process in the upper SF could be a mechanism for C-type inactivation (Wu et al., 2023). The SF of Shaker-family channel has the sequence “TVGYGD”. The first step involves a partial twist of the P-loop backbone, particularly involving the upmost SF residue D30’ (D379 in K_v1.2), which originally stabilizes itself by interacting with W17’ (W336 in K_v1.2). The second step is the reorientation of the upper-middle SF residue Y28’ (Y377 in K_v1.2) upward, which normally participates in hydrogen bonds with nearby pore-helix residues, to fill some of the original volume occupied by D30’.

In contrast, the hERG SF is characterized by the sequence “SVGFGN”. In the upper SF of the hERG channel, phenylalanine (F) replaces the tyrosine (Y) found in many K⁺ channels including Shaker and K_v1.2. Phenylalanine has a non-polar benzyl side chain, whereas tyrosine has a polar hydroxyl group attached to its benzene ring. The hydroxyl group in tyrosine can form additional hydrogen bonds, which may stabilize different SF conformations in other K⁺ channels compared to hERG. Similarly, at the outermost end of the selectivity filter of hERG, asparagine (N) replaces the aspartate (D) found in many other K⁺ channels. Asparagine is uncharged, while aspartate introduces a negative charge through its carboxylate group, which could explain the differences in ion coordination and gating dynamics. These structural differences may explain why hERG channel adopts a similar but distinct SF rearrangement compared to other K⁺ channels (such as Shaker and K_v1.2) and can have a slightly different structural mechanism of the C-type inactivation.

Molecular dynamics simulations were crucial in confirming conformational state functional properties. The observed differences in ion conduction between the predicted open and inactivated hERG channel states illustrate the impact of structural changes on ion channel activity. In addition, the drug docking analysis matched expectations, with significant preference for the inactivated-state model across a variety of hERG blockers. In most cases, the increased drug binding affinity to the inactivated state hERG model can be attributed to the protrusion of residue F656 aromatic side chain into hERG central cavity (Figure 3a). This protrusion has the potential to strengthen drug binding through providing extra π stacking or hydrophobic interactions. Additionally, the side chains of F656 may act as a barrier, blocking drugs that are bound deeper within the pore from escaping even when they do not have direct interactions with the drugs. Similar observations were made for the closed-state model, where the pore closure also positions the F656 side chains into the central cavity. In addition, the slanted orientation of the Y652 aromatic rings and repositioning of the S624 sidechain hydroxyl oxygens in the inactivated-state model enabled the formation of a binding pocket between S624 and Y652 (Figure 3b, at Z = 21 Å), enabling drugs to reach deeper into the pore.

Furthermore, we carried out molecular docking of multiple drugs to test their state-dependent binding preferences and relative affinities. An additional novelty of our study was the development of a new computational modeling and simulation approach that allowed us to use the predicted affinities of drug binding from docking simulations and compare them directly to measured hERG channel inhibition potencies from electrophysiology experiments. We employed a simulated hERG functional model comprising of five functional states that we have extensively validated in our earlier studies (Romero et al., 2015) and performed in silico electrophysiological experiments under the same conditions as reported in the experimental papers. In doing so, we could compute the relative hERG channel state probabilities during the experimental protocol. The channel state probabilities were then used to scale the computed affinities in each state from the docking simulations, allowing comparison of the overall computed relative affinity with experimentally reported relative potencies. Such an analysis created an opportunity for an “apples-to-apples” comparison between structurally derived affinity predictions and abundant functional measurements that have been conducted for half a century. This novel linkage can be readily expanded to other protein targets and any variety of drugs. The intersection of structural modeling, molecular docking and functional simulations, and supporting experimental data offers a comprehensive approach to understanding protein structures and their links to biological functions.

Despite the promising results, our study is not without limitations. While AlphaFold has demonstrated remarkable accuracy in numerous instances (AlQuraishi, 2019; Jumper et al., 2021b; Varadi and Velankar, 2023), it is important to note that the predicted models may not always be reliably accurate to assess drug binding (Karelina et al., 2023). Moreover, hERG channel inactivation and closure might encompass a range of states as was shown for other ion channels (Catterall et al., 2020; Goldschen-Ohm et al., 2013; Hite and MacKinnon, 2017; Li et al., 2024; Yao et al., 2023), and the conformations we have identified could potentially represent just a few possibilities within this broad spectrum. Recently, AlphaFold3 was released, incorporating a diffusion model that enables the prediction of proteins in complex with other proteins, small molecules, nucleic acids, and ions (Abramson et al., 2024). Although AlphaFold3 demonstrated high accuracy in various predictive tasks and is expected be widely adopted, the source code is not publicly available. Consequently, the approaches employed in this study with AlphaFold2 cannot yet be applied to AlphaFold3. The results from GALigandDock drug docking, while providing valuable insights, should be regarded as provisional at this juncture (Maly et al., 2022). Correlating simulated drug binding affinities with experimental results is challenging because multiple studies have shown that drug binding potency is highly dependent on the measurement technique used, resulting in different IC₅₀ values being reported for the same channel – drug pairing (Alexandrou et al., 2006; Asai et al., 2021; Cameron et al., 2021; Chiu et al., 2004; Ficker et al., 1998; Hansen et al., 2006; Johnson and Trudeau, 2023; Kushida et al., 2002; Orvos et al., 2019; Paul et al., 2002; Perrin et al., 2008; Rampe et al., 1997; Suessbrich et al., 1997; Tanaka et al., 2014; Walker et al., 1999; Zhang et al., 1999; Zhou et al., 1998) as was explored in detail in our recent study (Gomis-Tena et al., 2020). Additionally, knowing the binding free energies of a drug is not the complete story; binding rates such as K_v(association rate) and K_v(dissociation rate) are also crucial for a quantitative evaluation of drug binding to the channel. The results herein primarily serve as a basis for further detailed investigation using methods like molecular dynamics simulations linked to multiscale functional modeling, as demonstrated in our previous studies (DeMarco et al., 2021; Yang et al., 2020).

In conclusion, this study not only advances our knowledge of the hERG channel structural dynamics and state-dependent drug binding but also sets a precedent for using AlphaFold2 or similar computational methods, such as RoseTTAFold (Baek et al., 2021) and ESMFold (Lin et al., 2023), to study other ion channels and membrane proteins. The insights gained from this study are expected to have a lasting impact on the fields of ion channel physiology and pharmacology, including cardiac drug safety, providing a foundation for future research aimed at developing safer and more effective therapeutic agents.

Materials and methods

hERG channel model generation and refinement with AlphaFold2 and Rosetta

AlphaFold2 employs a deep learning architecture that integrates several innovative components including, residue pair representation comprising an architecture module that represents each possible pair of amino acid residues in the sequence. This pairwise representation captures the interactions between residues that determine the protein folding. AlphaFold2 also applies an attention mechanism, which constitutes a transformer-based model (similar to the architecture used in natural language processing) to weigh the influence of different parts of the input data differently (Jumper et al., 2021a). In AlphaFold2, the effect is to emphasize interactions between certain amino acid residues more than others based on how they might impact folding. There is also a so called Evoformer block within the learning model that specifically processes the evolutionary data from multiple sequence alignments, enabling the model to effectively incorporate evolutionary information (Jumper et al., 2021a). After processing through the Evoformer, intermediate representations are used to predict the distances and angles between residues as part of an iterative feedback process. One of the most critical pieces of the AlphaFold2 architecture is the recycling layers that promote an iterative refinement of predictions. The output from each iteration of the model feeds into the next one as input. This process, known as “recycling”, helps refine the accuracy of the predicted structure with each pass.

Due to the time-intensive process of creating multiple sequence alignments (MSA) for AlphaFold2, the ColabFold (Mirdita et al., 2022) software was made to streamline protein structure prediction by combining MMseqs2 (Mirdita et al., 2019) sequence search toolkit with AlphaFold2, enhancing runtime efficiency while preserving high prediction accuracy. First, we provided this hybrid structure as a template for ColabFold. ColabFold utilized only a portion of the MSA in the modeling process to reduce memory usage. We adjusted the limit to 256 cluster centers and 512 additional sequences (max_msa = “256:512”), a reduction from the standard values of 512 and 1024, respectively. Altering the random seed can lead to variance in cluster centers and, consequently, in structure predictions. Therefore, we set AlphaFold2 to cycle through 20 random seeds (each seed was used to produce 5 models) to foster diversity in the predicted conformational states (num_seeds = 20), resulting in 100 total models per run. An exception to this approach was the initial phase of inactivated-state hERG model creation, where a single seed number was used to generate 5 models. Additionally, we enabled the model stochastic component (use_dropout = True) during the inference pipeline to sample from an ensemble of models for the more ambiguous segments of the structure prediction. The number of iteration cycles was set to 20, with an early stopping tolerance of 0.5 based on per-residue confidence metric pLDDT deviation from the previous cycle. Structures were converted from the PDB file format to CIF prior to being used as templates in ColabFold. Previous studies indicated that these settings may permit the sampling of alternative protein conformations (Brown et al., 2024; del Alamo et al., 2022; Mirdita et al., 2022; Sala et al., 2023). The N-terminal PAS domain (residues M1 – R397) was not included in the modeling due to graphics card memory limitation making the resultant model (W398 – R863) resemble hERG 1b isoform (Phartiyal et al., 2007).

Structural data for the hERG channel in its open state (PDB 5VA2) (Wang and MacKinnon, 2017) were sourced from the Protein Data Bank (PDB) (Berman et al., 2000), and missing segments of the loops were modeled using Rosetta (Leman et al., 2020) software suite LoopRemodel (Huang et al., 2011). This reconstructed open-state model served as a basis for molecular dynamics and drug docking simulations. However, we also wanted to test the potential for AlphaFold2 to emulate the open state model. To do so, specific regions of the putative open state hERG model (PDB 5VA2) (Wang and MacKinnon, 2017) — namely the VSD (W398 – V549), the SF with adjacent pore helix (I607 – T634), a part of S6 helix and the cytosolic domain (S660 – R863) — were provided to ColabFold as structural templates. Then, 100 diverse models were generated for further analysis.

Post-prediction, the resulting 100 models for each structural state were categorized into clusters based on all-atom RMSD between those models. Closed-state models were clustered with a threshold of 0.75 Å across the entire channel, whereas inactivated and open states (non-empirical, AlphaFold-predicted) models focused on the SF (residues S624 – G628), with a more stringent threshold of 0.35 Å. We ranked the models in the cluster by their average per-residue confidence metric (predicted local distance difference test, or pLDDT), which assesses the likelihood that the predicted structure aligns with an experimentally determined structure (Jumper et al., 2021a). pLDDT above 90 are considered to be highly reliable and those between 70 and 90 as reliable with generally good protein backbone structure prediction (Jumper et al., 2021a). Lower scores indicate regions of lower confidence and may be unstructured. Cluster 1 is defined in this study to be the cluster with highest average pLDDT among all the clusters.

Afterward, we refined the preliminary atomistic structural models putatively representing each functional state of the hERG channel (open, inactivated and closed) using Rosetta FastRelax protocol (Fleishman et al., 2011; Leman et al., 2020) with an implicit membrane to optimize each residue conformation and resolve any steric clashes. The protocol was set to repeat 15 times and included an implicit POPC membrane environment. For each final model, 10 separate relaxation runs were executed, and the highest-scoring model from these runs was selected for further simulations and analyses.

Atomistic MD simulations of hERG channel conduction

The CHARMM-GUI web server (Jo et al., 2008) was employed to embed hERG channel structural models within tetragonal patches of phospholipid bilayers, each comprising approximately 230 to 240 POPC lipid molecules. The resulting assemblies were immersed in a 0.3 M KCl aqueous solution, culminating in molecular systems with an approximate total of 138,000 to 144,000 atoms. Residue protonation states were determined assuming a neutral pH=7.4 environment, and each channel subunit was terminated with standard charged N- and C-termini.

MD simulations were conducted using the Amber22 (Case et al., 2005) software suite. The simulations utilized standard all-atom Chemistry at Harvard Macromolecular Mechanics (CHARMM) force fields such as CHARMM36m (Huang et al., 2017) for proteins, C36 (Klauda et al., 2010) for lipids, and standard ion parameters (Beglov and Roux, 1994), in conjunction with the TIP3P water model (Jorgensen et al., 1983). The systems were maintained at 310.15 K and 1 atm pressure in the isobaric-isothermal (NPT) ensemble, facilitated by Langevin thermostat and the Berendsen barostat. Non-bonded interactions were calculated with a cutoff of 12 Å. Long-range electrostatic forces were computed using the Particle Mesh Ewald (PME) method (Darden et al., 1993), and van der Waals interactions were not subjected to long-range correction as per recommendations for the C36 lipid force field (Klauda et al., 2010). All hydrogen-connected covalent bonds were constrained using the SHAKE algorithm to enable a 2-fs MD simulation time step (Ryckaert et al., 1977).

Equilibration began with harmonic restraints imposed on all protein atoms and lipid tail dihedral angles, initially set at 20 kcal/mol/Ų and reduced to 2.5 kcal/mol/Ų over 6 ns. A subsequent 90 ns equilibration phase further decreased the restraints to 0.1 kcal/mol/Ų, initially encompassing all atoms in the protein and eventually focusing solely on the backbone atoms of pore-domain residues (G546 to F720). In selected MD simulations, an electric field was applied along the Z direction to mimic membrane voltage (Gumbart et al., 2012), increasing linearly over the final 10 ns of equilibration to reach either 500 or 750 mV. This setup prefaced a production phase lasting 910 ns, totaling 1 μs of total simulation time per each case.

Docking of small-molecule drugs to hERG channel models

Ligand structures (i.e., pharmaceutical compounds) were retrieved from the PubChem (Kim et al., 2016) and ZINC (Tingle et al., 2023) databases. In this study, we considered the protonation states of the top two most dominant species at the physiological pH of 7.4, computed using the Henderson–Hasselbalch equation. After these initial modifications, each ligand partial atomic charges, as well as atom and bond types, were refined using AM1-BCC correction via the Antechamber module in AmberTools22 (Case et al., 2023, 2005).

Each ligand was individually positioned within the pore of the hERG channel, specifically between the key drug-binding residues, Y652 and F656, located on the pore-lining S6 helices. This manual placement was a preparatory step for the docking process. Docking was executed using the GALigandDock (Park et al., 2021) Rosetta mover, a component of the Rosetta software suite (Leman et al., 2020). The DockFlex mode was utilized for this purpose with a spatial padding of 5 Å. For every individual hERG channel – ligand pair, a total of 20 parallel runs were conducted. Each run generated 1,250 unique docking poses, cumulatively providing a substantial collection of 25,000 poses for each ligand – hERG channel complex. This extensive array of poses was intended to ensure a comprehensive exploration of potential ligand binding configurations and orientations within the hERG channel pore.

Utilizing in-house Python scripts, we selected the top 50 docking poses for each ligand based on the lowest free energy of binding as calculated by GALigandDock. These selected poses were then clustered based on ligand all-atom RMSDs, factoring in the symmetry of the hERG channel to ensure accurate grouping. The clustering was executed with an RMSD threshold of 2.5 Å, and only clusters containing a minimum of 8 poses were deemed significant; those with fewer were classified as outliers. Within the cluster with the most favorable average binding free energy, we selected the top scoring pose as the representative pose for further analysis. This pose was then subjected to a detailed analysis of protein-ligand interactions utilizing the Grapheme Toolkit from the OpenEye software suite (www.eyesopen.com). The criteria for detecting interactions are outlined in the OEChem Toolkit manual (https://docs.eyesopen.com/toolkits/python/oechemtk/OEBioClasses/OEPerceiveInteractionOptions.html), with two modifications to minimize clutter: the MaxContactFraction is set to 1 (default: 1.2), and the MaxCationPiAngle is adjusted to 30° (default: 40°). The interaction patterns and binding sites were subsequently rendered as a two-dimensional image for comprehensive visual interpretation. Additionally, for a more detailed understanding of the spatial arrangement, three-dimensional visualization of the protein-ligand complexes was conducted using the ChimeraX (Pettersen et al., 2021) software.

Calculation of state probabilities using a five-state hERG Markov model

Over the time course of experimental recordings of hERG inhibition by various drugs, the channel can be in different functional states, each bound by drugs with different ionization states, making it difficult to compare experimental and simulated binding affinities. Moreover, different studies utilize different electrophysiological protocols to measure state-dependent ligand binding, further complicating comparisons. To address this complication, we used a five-state hERG Markov model to predict the probabilities of the channel in each state (open, closed, inactivated) during a given experimental protocol (Romero et al., 2015). Transition rate constants are provided in Table 4. The protocols that were used in the model to calculate each state probabilities (closed states: C3 + C2 + C1, inactivated state: I, and open state: O) are shown in Table 5. To simulate the inhibitory effects of drug on the hERG channel current, I_kr, we decreased the peak conductance, G_Kr in a concentration dependent fashion using a concentration response relationship with a Hill coefficient of 1 (n = 1) as follows:

Table 4.

Table 5.

where G_Kr,max is the nominal conductance value obtained from each ventricular myocyte model, [Drug] is a molar drug concentration, and the IC₅₀ is the concentration of drug that produces a 50% inhibition of the targeted transmembrane current, i.e. I_kr in this case (see Table 5).

To obtain the simulated binding affinity of a drug, ΔG_bind,sim, we multiplied the binding affinity for each state by the probability of the channel being in that state and combined these values for both neutral and cationic forms of the drug, as represented by the following equation:

Here, ΔG_{bind, O}, ΔG_{bind, I}, ΔG_{bind, C}represents simulated binding affinities for the open, inactivated, and closed state, respectively, to either the neutral or the cationic form of the drug. P_{hERG, O}, P_{hERG, I}, P_{hERG, G} represent the fraction of channels that are in the open, inactivated, and closed state, at the time when the tail current was observed in electrophysiological recordings to calculate drug fractional block, determined for the specific voltage protocol employed. P_{drug, neutral}and P_{drug, cationic} represent the fraction of neutral and cationic species of each drug at physiological pH 7.4 as calculated using the Henderson-Hasselbalch equation using drug pK_a values from Chemaxon. The zwitterionic species of moxifloxacin instead of cationic was included. These data are recorded in Table 3. Experimental IC₅₀ values (in units of M) were converted to equivalent binding free energies using the equation ΔG_{bind, exp} = –RTln(1/IC₅₀) where R = 0.0019872036 kcal K^-1 mol^-1 is the gas constant and T is the experimental temperature in K.

Molecular graphics and interaction analysis

Molecular graphics visualization was performed using ChimeraX (Pettersen et al., 2021). MD trajectory and simulation images were visualized using VMD (Humphrey et al., 1996). Interaction network analysis was performed using the Protein-Ligand Interaction Profiler (PLIP) (Salentin et al., 2015) with criteria outlined in Table 1.

Acknowledgements

We would like to thank all members of the I.V., C.E.C. and V.Y.-Y. laboratories and KN’s cats, Momo and Mimi, for helpful discussions and support. This work was supported by National Institutes of Health Common Fund Grant OT2OD026580 (to C.E.C. and I.V.), National Heart, Lung, and Blood Institute (NHLBI) grants R01HL128537, R01HL174001, R01HL085844, R01HL152681, and U01HL126273 (to C.E.C., V.Y.-Y. and I. V.), American Heart Association Career Development Award grant 19CDA34770101 (to I.V.), National Science Foundation travel grant 2032486 (to I.V.), UC Davis Department of Physiology and Membrane Biology Research Partnership Fund (to C.E.C. and I.V.) as well as UC Davis T32 Predoctoral Training in Basic and Translational Cardiovascular Medicine fellowship supported in part by NHLBI Institutional Training Grant T32HL086350 (to K.N.). Computer allocations were provided through Extreme Science and Engineering Discovery Environment (XSEDE) grant MCB170095 (to I.V., C.E.C., and V.Y.-Y.), Texas Advanced Computing Center (TACC) Leadership Resource and Pathways Allocations MCB20010 (I.V., C.E.C., and V.Y.-Y.), Oracle for Research fellowship and cloud credits award (to I.V., C.E.C.), Pittsburgh Supercomputing Center (PSC) Anton 2 allocations PSCA17085P, MCB160089P, PSCA18077P, PSCA17085P, PSCA16108P (to I.V., C.E.C., and V.Y.-Y.). Anton 2 computer time was provided by the Pittsburgh Supercomputing Center (PSC) through Grant R01GM116961 from the National Institutes of Health. The Anton 2 machine at PSC was generously made available by D.E. Shaw Research. OpenEye academic license was provided by OpenEye Scientific.

Data availability

All final study data are included in the article and/or SI Appendix with key molecular modeling, docking as well as molecular dynamics simulation and analysis data files and scripts available to download from Dryad digital repository at https://doi.org/10.5061/dryad.18931zd5x. (Note for the editor and reviewers: the link will be active after the manuscript acceptance for publication).

Additional information

Ethics declarations

Competing interests

The authors declare no competing interests.

Author Contributions

K.N. designed the research study, conducted simulations, acquired data, analyzed data, and wrote the manuscript. P.-C.Y., V.Y.-Y., C.E.C., and I.V. designed the research study, analyzed data, and wrote the manuscript.

Competing Interest Statement

The authors declare no competing interests.

Supplementary figures

Figure S1.

Figure S2.

Figure S3.

Figure S4.

Figure S5.

Figure S6.

Figure S7.

Figure S8.

Figure S9.

Figure S10.

Figure S11.