Upgraded molecular models of the human KCNQ1 potassium channel

Summary

The voltage-gated potassium channel KCNQ1 (K_V7.1) assembles with the KCNE1 accessory protein to generate the slow delayed rectifier current, I_KS, which is critical for membrane repolarization as part of the cardiac action potential. Loss-of-function (LOF) mutations in KCNQ1 are the most common cause of congenital long QT syndrome (LQTS), type 1 LQTS, an inherited genetic predisposition to cardiac arrhythmia and sudden cardiac death. A detailed structural understanding of KCNQ1 is needed to elucidate the molecular basis for KCNQ1 LOF in disease and to enable structure-guided design of new anti-arrhythmic drugs. In this work, advanced structural models of human KCNQ1 in the resting/closed and activated/open states were developed by Rosetta homology modeling guided by newly available experimentally-based templates: X. leavis KCNQ1 and various resting voltage sensor structures. Using molecular dynamics (MD) simulations, the capacity of the models to describe experimentally established channel properties including state-dependent voltage sensor gating charge interactions and pore conformations, PIP2 binding sites, and voltage sensor–pore domain interactions were validated. Rosetta energy calculations were applied to assess the utility of each model in interpreting mutation-evoked KCNQ1 dysfunction by predicting the change in protein thermodynamic stability for 50 experimentally characterized KCNQ1 variants with mutations located in the voltage-sensing domain. Energetic destabilization was successfully predicted for folding-defective KCNQ1 LOF mutants whereas wild type-like mutants exhibited no significant energetic frustrations, which supports growing evidence that mutation-induced protein destabilization is an especially common cause of KCNQ1 dysfunction. The new KCNQ1 Rosetta models provide helpful tools in the study of the structural basis for KCNQ1 function and can be used to generate hypotheses to explain KCNQ1 dysfunction.

Keywords:

Biology and life sciences – Biochemistry – Computational biology – Physical sciences – Chemistry – Research and analysis methods – proteins – Molecular biology – Macromolecular structure analysis – Database and informatics methods – Bioinformatics – Sequence analysis – Sequence alignment – Neuroscience – Simulation and modeling – Medicine and health sciences – Physiology – Lipids – Physics – Biological databases – Electrophysiology – Neurophysiology – Biophysics – Microscopy – Physical chemistry – Electron microscopy – Electron cryo-microscopy – Protein structure databases – Biochemical simulations – Protein structure – Chemical bonding – Hydrogen bonding

Introduction

Voltage-gated ion channels are ubiquitously expressed in human tissues and contribute to diverse physiological phenomena such as generation and modulation of the membrane potential in excitable cells, myocyte contraction, modulation of neurotransmitter and hormone release, and electrolyte transport in epithelia. KCNQ1 (K_V7.1) is a voltage-gated potassium (K_V) channel expressed in the heart and in epithelial cells of the inner ear, stomach, kidney and colon [1]. Like other K_V channels KCNQ1 is comprised of four identical α-subunits, each containing six membrane-spanning segments (S1-S6) and a pore loop (P loop) that contributes to the ion selectivity filter in the homo-tetramers to create the KCNQ1 channel. The central pore domain (PD, S5-P-S6) is surrounded by four voltage-sensing domains (VSDs, S1-S4) that respond to membrane depolarization by undergoing a conformational change that triggers structural rearrangements in the PD (called electromechanical coupling), which opens the channel gate making the channel conductive [2, 3]. VSD activation occurs stepwise and proceeds from an initial resting VSD conformation in which the PD is closed (resting/closed, RC) to an activated VSD with an open PD (activated/open, AO) [4, 5]. In KCNQ1, these transitions also involve an experimentally resolvable intermediate state [6, 7]. KCNQ1 also contains a sizable C-terminus that interacts with both phosphatidyl-4,5-bisphosphate (PIP2) and calmodulin [8, 9]. While modeling its structure is beyond the scope of this study, it plays important roles in regulating KCNQ1 channel function and contains many disease mutation sites [10, 11]. Finally, we note that a hallmark of the KCNQ1 channel is its co-assembly with the KCNE1 auxiliary subunit in the heart to generate the channel complex that is responsible for the slow delayed rectifier current (I_KS) necessary for myocardial repolarization [12, 13].

Heritable mutations in KCNQ1 are associated with several cardiac diseases including LQTS, atrial fibrillation, and short QT syndrome [14]. About 50% of the genetic cases of LQTS, which predisposes children and young adults to sudden cardiac death, are associated with dominant mutations in KCNQ1 (type 1 LQTS) [15]. While progress in the functional characterization of LQTS-associated mutations has been made [16–19], the molecular mechanisms underlying channel dysfunction remain difficult to assess without the availability of high-accuracy structural data. A detailed molecular understanding is needed to improve decision making for new unclassified KCNQ1 mutations and to support the development of new anti-arrhythmic therapeutics.

The lack of an experimentally determined structure for human KCNQ1 has prompted molecular modeling efforts. The first structural model of the KCNQ1 channel containing all membrane-embedded regions S1-S6 was published in 2007 [20] and since then has stimulated numerous structural and functional studies that provided new insights into KCNQ1, e.g. its regulation by KCNE1 [21–24] and other KCNE proteins [25, 26], binding of PIP2 [27, 28], the structural determinants and transitions in VSD activation [5, 29], and mechanisms of VSD-PD electromechanical coupling [30]. The recent determination of a cryo-electron microscopy (cryo-EM) structure of X. leavis KCNQ1 [31] which shares 78% sequence identity to human KCNQ1 provides a new template to upgrade the earlier Rosetta homology models of KCNQ1. The PD of this structure is believed to be in its closed state conformation, while its VSD occupies the activated state. In addition to this cryo-EM structure, recent experimentally determined structures of non-mammalian VSDs in resting conformations [32, 33] provide new templates for modeling the resting/closed channel, which previously was informed by purely computational models built by combining the bacterial KcsA template with de novo prediction methods [20, 34].

Here, using this new set of experimental structures leveraged in Rosetta multiple-template homology modeling, we report the development and structural validation of a second generation of Rosetta models of the human KCNQ1 channel in the closed and open states. These models have improved completeness and accuracy and thus represent useful tools for studying the structure-function relationships in KCNQ1 and the molecular mechanisms of mutation-induced changes in KCNQ1 phenotypes. We showcase the utility of the new models by (1) applying molecular dynamics (MD) simulation to our KCNQ1 models which generated useful insight into the VSD and PD conformations, channel-PIP2 interaction properties and VSD-PD coupling contacts, and (2) by applying Rosetta energy calculations to KCNQ1 VSD mutants to estimate the mutation-induced thermodynamic changes in protein stability, with results supporting the notion that protein destabilization is a common cause of mutation-induced KCNQ1 dysfunction. We anticipate that these Rosetta models will stimulate a new series of integrated structure-electrophysiology studies and personalized medicine approaches in KCNQ1 research.

Results

Model building and quality check of human KCNQ1 models

Structural models of human KCNQ1 (helix S0-S6) in the RC and AO state (Fig 1) were generated by Rosetta homology modeling guided by the multiple sequence alignments shown in S1 and S2 Figs. The template with the highest sequence identity (78%) to human KCNQ1 is the recently published cryo-EM structure of X. leavis KCNQ1 (PDB 5VMS) [31], which appears to represent a decoupled state with an activated VSD and a closed pore (referred to as the AC state). Thus, the two domains served as separate templates for the AO and RC states, respectively. The X. leavis KCNQ1 VSD provided a template for the human KCNQ1 AO model whereas its PD was used in RC state modeling. An additional template for the open pore conformation in the AO model was derived from the crystal structure of the chimeric K_V1.2–2.1 channel (PDB 2R9R) [35]. Because no experimental resting state VSD structure of a voltage-gated potassium channel was available, templates of the VSD in the RC model were obtained from VSD structures of related non-mammalian proteins and another structural model: the resting VSD of the C. intestinalis voltage-sensing phosphatase (Ci-VSP) (PDB 4G7Y) [32], VSD2 in the A. thaliana two pore calcium channel protein 1 (TPC1) (PDB 5DQQ) [33], and the resting VSD conformation C3 in a model of the Shaker channel [36].

Rosetta homology modeling yielded an ensemble of low energy models that were assessed with Procheck [37] and Molprobity [38]. The RC and AO models with the best Molprobity statistics (compare with S1 Table) were selected as final models. They exhibited favorable Molprobity and clash scores that were equal or better than the 90^th percentile of structures in the Molprobity database. No violations for backbone bonds or Cβ deviations were observed, and only four bond angles in both models fell outside of four times the standard deviation of dictionary values. In the AO model, all residues fell within favored (95.9%) or allowed (4.1%) regions of the Ramachandran plot. In the RC model, only 1% of the residues fell in disallowed regions with the remaining being in favored (92.8%) and allowed (6.2%) regions. Analysis of the side-chain conformations revealed 100% (AO model) and 99.5% (RC model) favored rotamers and no bad rotamers. Taken together, these values are well within the range of high-resolution structural models and a positive indicator of the quality of our KCNQ1 models. PDB files of these models are provided in the Supporting Information (S1 PDB and S2 PDB). Morphs between the RC and AO model are shown in supporting videos S1 Movie, S2 Movie and S3 Movie.

Structural validation of human KCNQ1 RC and AO models

Our KCNQ1 RC and AO models are displayed in Fig 1 and their topology is compared to that of X. leavis KCNQ1. KCNQ1 has a domain-swapped homotetrameric topology with the VSD of one subunit interacting with helix S5’ from a neighboring subunit. The VSD and PD in the AO model had a Cα-atom RMSD of 1.9 Å and 2.3 Å relative to the X. leavis KCNQ1 AC structure, respectively. For the RC model the Cα-RMSD for the VSD and PD were 4.1 Å and 1.0 Å, respectively. These values are considerably lower than those calculated for the first generation of human KCNQ1 models [20] which had RMSD values for the VSD and PD of 4.8 Å (3.5 Å when excluding the S2-S3 linker) / 3.0 Å (AO model) and 4.3 Å (3.6 Å) / 3.1 Å (RC model), respectively. Thus, the new KCNQ1 models are closer to the X. leavis KCNQ1 structure (i.e. the most closely related structural homolog available), which is a positive indicator for their improvement over the Rosetta models by Smith et al. [20] (see Discussion).

The KCNQ1 RC and AO models were validated rigorously against known experimental and structural data to assess model confidence. First, we calculated the channel pore radius and compared it to the closed PD conformation in X. leavis KCNQ1 (Fig 2). In the RC model, the positions with the narrowest constriction along the ion conduction pathway are at residues G345, S349 and L353 which correspond to G335, S339 and L343 in X. leavis KCNQ1. The pore radius at those residues (~0.5 Å) is clearly below the ionic radius of a K⁺ ion (1.38 Å) confirming that the channel gate in the RC model is closed. In contrast, the channel gate is open in the AO model. The radius at the narrowest point of the pore at A344 is not smaller than 2.0 Å, which will allow a K⁺ ion to pass through. Beyond that point, the pore radius increases and remains larger than the K⁺ ionic radius. This transition point in the KCNQ1 pore coincides with the functionally important PAG motif, which corresponds to PXP in other K_V channels. The PAG motif serves as a hinge point allowing the S6 helix to bend and swing outward, which opens the channel gate (S4 Movie). Mutations of the PAG motif abolish current but not channel surface expression [39] consistent with the idea that this region is important for movement of the activation gate. These observations for our Rosetta models were reproduced by pore radius measurements made for these models when employed in MD simulations (S3 Fig, see explanations below), further confirming the validity and stability of the PD structure in our KCNQ1 models.

As a second validation, we measured state-dependent distances in our models between amino acid residues suggested to be responsible for voltage-dependent activation of the VSD and compared them to expectations from known residue pairings inferred from electrophysiological studies on charge-reversal VSD mutants of KCNQ1 [5, 6]. Positively charged residues located along the S4 helix, commonly referred to as gating charges and labeled R1 through R6, undergo a series of stepwise transitions and successively pair with acidic residues on S2 and S3 that creates movement of S4 and confers voltage-sensitivity to the ion channel [40] (S5 Movie). Compared to the Shaker [36] and K_V1.2 channel [41], it is currently assumed that KCNQ1 adopts a limited number of VSD states (resting, intermediate and activated) [42], possibly due to the lower number of positively charged S4 residues. Residues R3 and K5 in K_V1.2 are replaced by Q3 and H5 in KCNQ1. Critical residues that coordinate S4 movement include E1 (E160/E150 in human/X. leavis KCNQ1), E2 (E170/E160) and D202 (D192 in X. leavis KCNQ1). In the resting VSD, E1 likely pairs with R1 (R228 in human KCNQ1) whereas the activated conformation is stabilized by an interaction between E1 and R4 (R237) [5, 6, 43]. In addition, residue Q3 (Q234/Q224 in human/X. leavis KCNQ1) is found close to E2 in the cryo-EM structure of KCNQ1, which likely represents an activated state of the VSD. The described state-dependent residue pairings are correctly captured by our KCNQ1 Rosetta models and persistently detected in MD simulations of the RC and AO model (Fig 3; see also S5 Movie for a morph between the resting and activated VSD structure). The structural models also suggest additional interactions occur between helix S4 and residues of the charge transfer center: in the RC model R2 and Q3 interact with E2 and D202 whereas in the AO model H5 and R6 are close to those two acidic residues. Moreover, another source of structural stabilization of the VSD in its resting and activated conformation may come from interactions of the gating charges with membrane phospholipids as observed for the VSD of other K_V channels [44–46]. This idea is supported by the observation of frequent hydrogen bond contacts of R6 (RC) and R1 (AO) with phospholipids in the lower and upper bilayer leaflet during MD simulation (Fig 3).

Overall, the measured channel pore dimensions and the good agreement of residue pairings in the VSD with experimental data lend confidence to the reliability of our KCNQ1 structural models. Further, they provide insight into the structural dynamics of VSD activation and channel opening and closing.

Refinement of KCNQ1 models by MD simulation and study of KCNQ1-lipid interactions

Computational simulations have become an indispensable tool in the study of voltage-gated ion channels and have provided insights into different aspects of KCNQ1 function and regulation, such as KCNE1 binding [21, 49, 50], drug [51] and lipid binding [52]. Here, we assessed the structural stability of our KCNQ1 AO and RC homology models in all-atom MD simulations with Amber [48]. We further validated that the simulations can reproduce known interactions of KCNQ1 with PIP2 previously revealed by experimental data. PIP2 is an important second messenger for cell signaling that binds to and regulates a wide variety of ion channels, including KCNQ1 [53]. Depletion of PIP2 from the membrane suppresses VSD-PD coupling in KCNQ1 [27] but leaves the VSD activation intact.

A total of four RC state and four AO state simulations starting from the selected KCNQ1 models as well as three additional conformations from the ensemble of Rosetta homology models were conducted, each running for 400 ns. Representative snapshots from the trajectory of AO and RC models are provided as supporting files S1 MDtraj and S2 MDtraj, respectively. Simulations were performed under constant temperature and pressure conditions in explicit POPC bilayers that contained 10 mol% (~28 molecules) of PIP2 in the cytosolic leaflet (Fig 4A). Overall, the KCNQ1 models remained stable in MD and deviated by not more than a Cα-RMSD of 4–5 Å from their respective starting conformations (Fig 4B and S4 Fig). The RC and AO models maintained their closed and open pore structure (Fig 2 and S3 Fig), respectively, as well as their characteristic gating charge pairings (Fig 3) as discussed above.

The number of hydrogen bonds with POPC and PIP2 was analyzed after letting the system equilibrate for 100 ns. Because the mean square displacement values for lateral diffusion of POPC and PIP2 in Amber lipid simulations are on the order of several ten Ų per 10 ns [54], this time span was long enough to allow the initially randomly placed PIP2 molecules to diffuse towards KCNQ1 and engage in close protein contacts (Fig 4C). PIP2 was observed to bind preferentially in the cleft between two neighboring VSDs (Fig 4C). It made numerous hydrogen bond interactions with residues on the S2-S3 linker (S2-S3L) (R181, K183, R190, R192), the S4-S5 linker (S4-S5L) (R249, R259) and the cytoplasmic end of S6 (S6_C) (K354, R360, K362) (Fig 4D) in agreement with the experimentally reported PIP2 binding sites [27, 28]. Previous studies suggested that interactions of PIP2 with these positively charged residues strengthen the coupling between the VSD and PD and are critical for pore opening [55]. Consistent with its high functional importance, a series of LQTS-associated mutations affect the PIP2 binding region: K183R [56], K183M [57], R190Q [58], R190L [59], R190W [60], R192H [61], R192P [60], R259C [62], R259H [61], R259L [63], R360M [59], R360T [60].

Comparison of the locations of POPC and PIP2 interaction sites on the cytosolic juxtamembrane side of KCNQ1 (Fig 4D) revealed some significant differences: POPC was unable to participate in hydrogen bonding with residues R181 and K183 in the short loop between S2 and the S2-S3 linker and with residues on S6c which are located within or below the water-interface region of the membrane. As a consequence of the small PC head group those residues are outside of the range of POPC but can still be spanned by PIP2, with its 1-phosphoinositol-4,5-diphospate headgroup.

In conclusion, all-atom MD simulations of our KCNQ1 models correctly identified experimentally established PIP2 binding sites in S2-S3L, S4, S4-S5L and S6_C, and revealed the molecular details of KCNQ1-PIP2 interaction. Owing to its increased head group size, PIP2 can reach and bridge sites in the VSD (S2-S3L, S4) with sites in the PD (S4-S5L, S6_C) to act as a molecular coupling element for KCNQ1.

High impact mutation sites map intra -⁠ and inter-subunit contact regions between the VSD and PD in KCNQ1 models

Using MD simulation, we further characterized the interactions between the VSD and PD in the RC and AO model and compared them to reported experimental data [30, 64–66]. VSD-PD interactions translate the voltage-dependent movement of S4 to channel opening and closing. Inference from studies of other K_V channels [67–69] and mutational scanning experiments on KCNQ1 [30, 64] have pointed to the interface between the S4-S5 linker (S4-S5L) and the cytosolic end of S6 (S6_C) as one important region for electromechanical coupling. In those studies, one group of mutations in S4-S5L [30] and S6_C [64] reduced the channel opening rate and shifted channel activation to more depolarized voltages whereas another group of mutations, specifically at V254 in S4-S5L and at L353 in S6_C, prevented the channel from closing and led to a constitutively open channel. These results were taken to suggest that S4-S5L binds to S6_C to stabilize the closed state of the channel and that S4-S5L relocation in voltage-dependent activation releases tension on S6 that allows it to kink, promoting channel opening.

In our RC and AO models, we identified extended regions of intra -⁠ and inter-subunit contacts (displayed in Fig 5A) that tie the VSD and S4-S5L to the neighboring PD. Intra-subunit contacts were found between S4-S5L and S6_C and between S5 and the N-terminal end of S6 (S6_N). Inter-subunit contacts were observed for S4 and S4-S5L with S5’ from an adjacent KCNQ1 chain and between two S6 helices from neighboring chains that cross each other near the channel gate. Strikingly, mutation sites in S4-S5L, S5 and S6 with high impact on channel opening and closing [30, 64, 66] are clustered in those contact regions and have consistently high contact numbers (Fig 5B and 5D) indicating that those regions are important for channel function. For clarity, mutation sites in the selectivity filter and P helix were not marked in Fig 5B. Those elements share only a small contact area with the VSD and many mutations cause complete channel loss of function (LOF) because of the detrimental effect that these mutations have on the structural integrity of the selectivity filter and extracellular pore entrance gate.

In the RC model, the S4-S5L –S6_C contact surface is formed by residues L250, L251, V254, V255, H258 on S4-S5L and I346, L347, G350, F351, L353 and K354 on S6_C. The side-chains of V254 and L353 were found to be proximal to each other in the model consistent with the results of combinations of mutations at those two sites [30] that were able to rescue channel function. Interestingly, L353 also engages in side-chain–backbone interactions with one of its symmetry mates in the neighboring subunit. As viewed from the cytosolic side (S4 Movie), the four symmetric copies of L353 are arranged circularly with the orientation of their side-chains directed clockwise. L353 is located one helical turn below the helix-helix crossing point at S349, which forms the narrowest constriction of the channel pore (compare with Fig 2). Hence, our structural model can provide a conclusive explanation for the constitutively open channel phenotype of KCNQ1 mutants with charged amino acid substitutions (K/E) at L353 [64]: steric repulsion between close local charges at the position of L353 pushes the S6 helices apart and prevents the channel from closing.

The contact region between S4-S5L and S6_C is maintained in the AO model as well and comprises the same group of residues but is reoriented slightly when S4-S5L and S6_C slide against each other during pore opening. The L353 sites, however, become separated from each other when S6 is kinked (S4 Movie) and are not in contact in the AO model.

The inter-subunit contacts span similar regions in the RC and AO model, respectively, but differ in their local structure reflecting the conformational changes that occur during channel activation as described below. Inter-subunit contacts are probably also the reason for the global motions in KCNQ1 that we observed, as characteristic regular patterns in the dynamic residue cross-correlation matrix (DCCM) (S5 Fig), and by performing principal component analysis (PCA) on the intra-membrane segments in KCNQ1 (S5 Fig). The global motions of KCNQ1 in MD simulation feature a rigid-body-like swing movement of the VSDs that includes both an ‘in-membrane-plane’ as well as an ‘out-of-membrane-plane’ component. Two neighboring VSDs were observed moving antiparallel to each other whereas the two domains on opposite sides of the tetrameric channel moved in the same direction (see S5 Fig for further details).

We further noticed that many mutation sites that have been associated with LQTS fall within intra -⁠ and inter-subunit contact regions in KCNQ1 (see Fig 5C). A simple classifier that distinguishes LQTS from non-LQTS sites based on residue contact number yields a receiver operating characteristic (ROC) with an area under the curve (AUC) of 59% and 58% for the RC and AO model, respectively (Fig 5E). This shows that LQTS sites are enriched in areas with high contact number and that mutations at those sites are critical to ensure channel function.

Taken together, we show that the VSD and PD in KCNQ1 share extended intra -⁠ and inter-subunit contact regions between S4, S4-S5L, S5 and S6. Experimentally tested high-impact mutation sites and sites associated with LQTS fall in those channel regions, which highlights their importance for channel function and renders them possible candidates for roles in electromechanical coupling.

Transition from the RC to the AO state triggers changes in inter-subunit contacts between the VSD and pore

Comparison of the RC and AO model revealed conformational changes at the VSD-PD interface and in the PD (Fig 6) (in addition to S4 in the VSD), which likely reflect the changes happening during depolarization-activated channel gating. Triggered by the gating charge movement in the VSD, S4 slides along S5’ and moves upwards by about two helical turns. Hydrophobic residues on the S5’-facing side of S4 swap their positions and come into contact with new residues on S5’ while S4 is moving (compare with left plot in Fig 6A). The S4-S5L is pulled upwards and sideward (Fig 6B and 6C) and engages in new interactions with residues located two to three turns higher on S5’. For instance, in the RC model, T247 makes contacts with the first residues in S5’ (Q260, E261) but is packed against T264, Y267 and I268 in the AO model. W248 makes few interactions with S5’ in the RC model but makes many putative stabilizing contacts with Y267, I268 and L271 in the AO model (Fig 6B). Moreover, movement of S4-S5L allows S6 to kink and move sideward cantilever-like. This movement breaks contacts between the ends of S6_C and induces a crosswise sliding of the S6 helices (compare with right plot in Fig 6A) that widens the channel pore diameter. Specifically, residue F351 slides upwards by about two helical turns and becomes tightly packed in the AO model where it appears to bridge sites on S4-S5L with S6_C of the same and S6’_C of the neighboring subunit (Fig 6B). The observation that F351 is at a critical spot in KCNQ1 with many contacts to other residues is consistent with the fact that F351 is highly conserved in K_V channels, and that it is a high impact mutation site [64]. Experimental studies [64] showed that the F351A mutant featured a rightward shifted current-voltage relationship and longer activation and deactivation times that mimicked the I_KS-phenotype of KCNQ1, and was ascribed to a destabilization of the KCNQ1 intermediate/open state [6] similar to the effects of KCNE1.

Hou and coworkers [7] demonstrated that KCNQ1 when expressed alone can also undergo an opening transition once the VSD is in the intermediate or activated conformation. To characterize the structural changes associated with channel opening from an activated VSD, we further compared the AO model with an additional model for the KCNQ1 AC state. The AC model was generated by comparative modeling based on the X leavis cryo-EM structure [31], and is available as supporting file S3 PDB. Residue contacts were analyzed by conducting four independent 350 ns all-atom MD simulations in POPC bilayers. The AC model maintained its activated VSD and closed PD conformation during MD (S3 and S4 Figs). The S4 helix is in an upward position as in the AO model, but the intracellular channel gate is closed with the narrowest constrictions formed by residues S349 and L353 (S3 and S6 Figs). At its N-terminal site, S4 makes fewer contacts with S5’ in the AC model (S6A and S6B Fig) suggesting that the VSD-PD interface is less tightly packed in the AC state. As seen from within the membrane plane, S4-S5L adopts a similar position in the AC and AO model, but is oriented slightly different and tilted inwards and closer to S6 in the AC model (S6B and S6C Fig). Inward tilting of S4-S5L likely stabilizes the closed PD conformation in the AC state by holding neighboring S6 helices together and increasing the number of inter-subunit contacts (S6A and S6B Fig). This observation suggests that in the AC state the channel is primed to open but that an additional outward tilting of S4-S5L is needed before S6 can bend to facilitate widening of the channel pore. We speculate that the separation of S4 and S4-S5L movement as inferred from our modeling of the AC and AO states and supported by the X. leavis KCNQ1 structure [31] is the result of decreased electromechanical coupling in KCNQ1 which is strongly dependent on PIP2 binding [6, 27, 55].

In conclusion, we characterized VSD-PD interactions and their state-dependent changes in KCNQ1 models when the channel transitions from the closed (RC or AC) to the open (AO) state. Our observations provide a stepping stone for systematically testing protein-contact-mediated coupling mechanisms in the activated channel state.

Rosetta-predicted stability changes in KCNQ1 models correlate with experimental channel loss of function data

While progress in the functional characterization of LQTS-associated mutations has been made [16–19, 70–72] the mechanistic molecular basis of channel dysfunction for most KCNQ1 mutations is still unclear. The dearth of direct experimental data that can help to determine how mutations alter KCNQ1 structure and function has prompted computational modeling approaches [20, 73].

Here, we assessed the capability of our structural models to predict the change in thermodynamic stability (ΔΔG) in KCNQ1 variants and whether the predicted energy changes can account for alterations in channel function. We used a Rosetta ΔΔG prediction protocol [74] allowing for flexibility in the protein backbone and side-chains around the mutation site. Protein backbone flexibility was modeled by a series of backrub moves [75] that rotate short main-chain segments (three to twelve residues) as rigid bodies about an axis defined by the starting and ending Cα-atom of the segment. This aims at alleviating the prediction for mutations with a large change in amino acid size. Side-chain flexibility was modeled using discrete rotamer conformations and simulated annealing (i.e. side chain packing). The ΔΔG was calculated as mean difference between the three top-scoring models of wild type (WT) and mutant KCNQ1 and averaged over ten different starting conformations taken from the ensemble of Rosetta homology models as described in Methods.

ΔΔG predictions were made for a set of single-site KCNQ1 variants with the affected sites located in the VSD (shown in Fig 7A, and divided by functional classes in S7 Fig). These KCNQ1 VSD variants were previously biochemically, biophysically, and functionally characterized [19], which provided the mechanistic basis for their pathogenicity. Comparison with this variant dataset allowed us to validate our computational stability predictions and assess the capabilities and limitations of a structure-based approach to predict KCNQ1 dysfunction. Not included from this original set of 51 KCNQ1 variants was the deletion mutant ΔF167 because modeling of deletions and insertions is currently not supported by the Rosetta ΔΔG prediction protocol. Furthermore, ΔΔG predictions for mutations to proline were found to be off-scale due to incompatible backbone torsions yielding Ramachandran and ring closure penalties as noted earlier [76]. Therefore, their ΔΔG values were considered less reliable and also not included in the analysis leaving ΔΔG data for 44 single-site KCNQ1 variants.

Fig 7B shows the distribution of ΔΔG values for those variants calculated with the KCNQ1 AO model as a function of variant class. Variants were previously assigned one of six classes based on their channel phenotype: WT-like and classes I-V for variants that have perturbations in the expression level, surface trafficking, and/or channel current density (see also the caption for Fig 7 and the footnote of S2 Table for a description of the class phenotype). Roughly, the severity of the channel LOF increases from class I to V with variants in class IV and V having the most severe perturbations. The median ΔΔG value of these two classes was significantly higher (ΔΔG˜IV=10.7REU (p_IV < 0.05), ΔΔG˜V=26.9REU (p_V < 0.001)) than that one of WT-like variants (ΔΔG˜WT‐like=0.72REU) indicating the mutant protein is thermodynamically destabilized. It must be mentioned that Rosetta ΔΔG values (in Rosetta Energy Units, REU) do not necessarily have the same scale as experimentally measured protein stability values (in kcal/mol), but that their relative order (e.g. between different variant classes) is more instructive.

The predicted ΔΔG values were correlated with the expression, trafficking and current density data (Fig 7C) suggesting protein stability is an important factor that contributes to these channel properties, and protein destabilization is a main cause of channel LOF. This idea is supported by NMR measurements on the KCNQ1 voltage sensor [19], which for two of seven class IV variants and ten of 14 class V variants had revealed spectra that were characteristic for a folding-destabilized protein. Another four class V variants failed to express the voltage sensor in E. coli also indicating protein folding defects [19]. The channel state of the structural model (RC or AO) had only minor influence on the ΔΔG predictions. Similar results were also obtained with the RC model (S8 and S9 Figs), indicating that most of the investigated mutations affected both channel states to a similar degree. Examples of mutants with varying ΔΔG predictions in the RC and AO model are shown in supporting S10 Fig and further discussed in the caption to S10 Fig.

The correlation plots in Fig 7C seem to have a complex, non-linear dependency, specifically the plot for peak current density, suggesting additional factors other than protein stability affect the named properties and/or that some variants had incorrect ΔΔG predictions. Variants that fell outside the 95% confidence interval for a linear model included four class I variants (C122Y, A150T, K218E, I227L), three class V variants (G179S, W176R, H126L), A102S and A149V (both WT-like) as well as P197L (class II) and R109L (class III). Structural models for these variants are shown in S11 Fig in the supporting information and their experimental and ΔΔG data are listed in S2 Table. In particular, it was noticeable that class I variants had ΔΔG values (ΔΔG˜I=3.4REU or ΔΔG˜I=1.4REU with the AO or RC model, respectively) not significantly different from WT-like variants (ΔΔG˜WT‐like=0.7REU), which failed to explain why these channels have low current density. For example, C122Y was predicted a negative ΔΔG (i.e. more stabilized than WT) but this channel was nonconductive. Two other variants for which channel function was not correlated with thermodynamic stability were A102S and R109L which had normal (0.2 REU) or increased ΔΔG values (10.3 REU), respectively, but exhibited current densities more than two times larger than WT (i.e. these were gain-of-function mutations).

The other few cases where the stability prediction was inconsistent with the functional data were likely due to incorrect ΔΔG calculations, which was supported by visual inspection of the respective models (see S11 Fig). For instance, A149V was predicted to have a high positive ΔΔG but had only moderately reduced expression levels and a trafficking efficiency higher than WT. Visual inspection of the structural model of this variant (S11 Fig) indicated remaining steric clashes between the valine side-chain and the nearby protein backbone causing a high repulsive score. This suggests that the backrub moves that are performed in the Flex ddG protocol prior to amino acid mutation had optimized the packing of the native alanine residue with the neighboring S3/S4 loop and generated a backbone conformational ensemble that was not sufficiently diverse to accommodate the larger valine side-chain upon mutation. Mutants G179S, W176R and H126L had ΔΔG values that were lower than other class V variants, mostly due to score balancing interactions made by the substituted amino acid residue in the mutant model, inconsistent with their low expression levels, trafficking and conductance.

In summary, Rosetta energy calculations correctly predicted destabilizations for most of the folding-defective KCNQ1 LOF variants: 13/18 variants with low expression (<65%WT) had increased ΔΔG values above 10 REU. Our results support growing evidence [19] that mutation-induced protein destabilization is a common cause of KCNQ1 dysfunction. They further suggest other so far unknown mutational effects exist (e.g. in class I variants), which are not detected in protein stability calculations, but have an impact on KCNQ1 function.

Discussion

Improvements of upgraded KCNQ1 Rosetta models

Our ability to improve homology models of human KCNQ1 was largely informed by the new X. leavis KCNQ1 template, the closest structural homolog available. Consequently, the Rosetta models presented here bear high similarity to the X. leavis KCNQ1 cryo-EM structure and are closer to it than the Rosetta models by Smith et al. [20]. The backbone RMSD for the VSD / PD compared to X. leavis KCNQ1 is 1.9 Å / 2.3 Å (this study) and 4.8 Å / 2.9 Å (Smith et al. [20]), respectively. For the closed state, the RMSD values for the VSD / PD are 4.1 Å / 1.0 Å (this study) and 4.3 Å / 3.1 Å (Smith et al. [20]), respectively. This indicates considerable remodeling changes were applied to the Rosetta models: comparison of the models presented herein and in Smith el al. shows an RMSD for the VSD / PD of 5.3 Å / 3.2 Å (closed) and 5.2 Å / 2.4 Å (open). It is important to mention that the Rosetta models of this work also compare favorably with their other structural templates: the PD in the AO model with the PD in K_V1.2–2.1 (1.9 Å) and the VSD in the RC model with the VSD in Ci-VSP (3.6 Å for sequence-aligned regions; 2.6 Å for S1+S2) and TPC1 (4.4 Å for sequence-aligned regions; 2.1 Å for S1+S2). These results, together with the favorable stereochemical properties of the new human KCNQ1 Rosetta models and their valid VSD and PD conformations as illustrated in Results, lends confidence that they represent a substantial improvement over the preceding models. It is also important to point out that the Rosetta models by Smith et al. did not include the S0 helix in the VSD, which is now appreciated to represent a crucial element for structural stability (see reference [19] and Results on KCNQ1 VSD mutations). The previous models further fail to satisfy the E1-R1 (E160-R228) gating charge interaction in the RC state which is an experimentally established contact stabilizing the resting state [5, 6], however, was not known for KCNQ1 at the time of the first model generation.

The improved Rosetta models should be regarded as working models that are subject to revision as more experimental structural data become available. Specifically, we need to point out a few limitations. Of lower confidence are those regions for which no coordinate information in the template structures was available, requiring de novo modeling. Furthermore, higher uncertainty applies to transition regions where templates from different sources were hybridized. Lower confidence regions include the loop between S4 and S4-S5L in the RC and AO models, the loop between S3 and S4 in the AO model and the loops between S0 and S1 and between S2-S3L and S3 in the RC model.

Implications for KCNQ1 structural biology questions

Given improved completeness and accuracy, the Rosetta models provide helpful tools to refine and advance our understanding of the structure-function relationships in KCNQ1. A hallmark of the KCNQ1 channel is its distinctive functional properties observed upon association with different members of the family of KCNE accessory proteins. While KCNE1 increases current density and delays channel opening and closing [12, 13], KCNE2 constitutively activates KCNQ1 and reduces current levels to approximately 50% of KCNQ1 alone [77]. KCNE3 also makes KCNQ1 constitutively open but increases current by ~10 fold [26, 78].

The mechanism by which KCNE1 modulates KCNQ1 is a longstanding yet still unsettled debate and different models have been put forward such as alteration of the VSD movement [72, 79–83], perturbation of pore opening [84] and changes in VSD-PD coupling interactions [6, 81, 82]. A mapping onto our KCNQ1 Rosetta models of putative KCNE1 binding sites identified in cysteine crosslinking [22, 23, 84–86], mutational scanning and double mutant cycle [25, 87] experiments (see S12 Fig) suggests very similar binding pockets for KCNE1 in the RC and AO state. Preliminary docking calculations with our Rosetta models also suggest that KCNE1 makes similar interactions with KCNQ1 in the RC and AO states and adopts similar orientations (S12 Fig), which is an unexpected result and was not observed in the original Rosetta models [21] indicating that a revision of the old KCNQ1/KCNE1 models is well-merited.

Of the many previous structure-function studies of KCNE1, we want to highlight here the model proposed by Zaydman et al. [6] who suggested that KCNE1 modulates KCNQ1 function by changing the VSD-PD coupling interactions. This idea stems from observations of KCNE1 inhibiting channel opening when the VSD is in the intermediate conformation between resting and activated state, as well as studies of F351A, a mutation that impairs coupling by weakening the contact between S4-S5L and S6_C, mimicking the KCNE1-phenotype of KCNQ1. Combining electrophysiology and kinetic modeling, Zaydman et al. concluded that changes in state-dependent VSD-PD interactions may explain many of the effects of KCNE1 including its delay of channel opening, increased channel conductance and pharmacological properties. Consistent with this idea, we identified extended intra -⁠ and inter-subunit contact regions between the VSD and PD in our Rosetta models (Fig 5) that potentially contribute to electromechanical coupling in addition to the established S4-S5L –S6_C contact. Several other studies [6, 7, 27, 55, 81] hinted at the role of protein-protein and protein-PIP2-dependent mechanisms in KCNQ1 electromechanical coupling. It appears to be intrinsically weak and requires PIP2 binding to strengthen the interaction between the VSD and PD [6, 27, 55]. This is evidenced by the cryo-EM structure of X. leavis KCNQ1 [31] which was determined in the absence of PIP2 and represents a decoupled state. These data led to the conclusion that PIP2 acts as important coupling element which is especially prominent in KCNQ1. Recently, Hou and coauthors [7] also proposed different protein-protein coupling mechanisms may exist in KCNQ1 for channel opening from the intermediate and activated state of the VSD, and that these differences underlie the inactivation phenotype of KCNQ1. We speculate that KCNE1 modifies these VSD-PD interactions, changing the properties of the open pore conformation and affecting the channel open probability. Previous docking studies [21, 49, 87] converged on a model in which KCNE1 is bound in a crevice between the VSD and PD of two neighboring subunits. The intimate interactions that KCNE1 makes in these models with several sites on S4, S4-S5L, S5 and S6 seems well suited to modulate electromechanical coupling interactions.

Explicit structure-function knowledge concerning KCNE2 binding and regulation of KCNQ1 is limited. Results from cysteine scanning [88], double mutant cycle analysis and computational docking [25] suggested a different interaction mode for KCNE2 than for KCNE1. However, a widely recognized molecular model for the KCNQ1-KCNE2 interaction has yet to be established. The KCNQ1 Rosetta models provide a means to devise new experimentally-testable hypotheses on KCNE2’s mode of action that should trigger further progress in this field.

For KCNE3, a series of functional studies [89, 90] have converged on the appreciation that it acts primarily by stabilizing the VSD, especially S4, in its fully activated state, thereby increasing the open probability at resting potential and making KCNQ1 a constitutively open channel. The stabilization of S4 arises from salt bridge interactions between R228 on S4 and residues D54 and D55 on KCNE3 [89] and possibly by another polar contact between Q244 on KCNQ1 and R83 on KCNE3 [26]. Our Rosetta model of the AO state supports this mechanistic explanation by showing R228 and Q244 at surface-accessible sites at the extracellular and cytosolic ends of S4 where those residues extend into the inter-subunit crevice formed between the VSD and PD. Binding of KCNE3 to R228 on the extracellular side and to Q244 on the cytosolic side is suggested to lock the S4 helix in place and stabilize the activated VSD conformation. Recent studies also suggest an additional effect of KCNE3 on the PD [81] which can be investigated with the new KCNQ1 models e.g. through docking and MD simulations.

Potential use of KCNQ1 models for personalized medicine

For a variety of genetic human disorders, substantial evidence suggests that a majority of disease-causing mutations act by destabilizing the folded conformation of the encoded protein [91, 92] thereby reducing the population of the functional form of the affected protein in cells. Predicting the phenotype of a mutation and interpreting variants of unknown significance discovered by genome sequencing data can greatly benefit from methods to quantitatively compare the stabilities of WT and mutant forms of proteins. For water-soluble proteins, computational algorithms are able to predict mutation-induced stability changes with reasonable accuracy [93]. However, these methods perform significantly worse when applied to membrane proteins [94], highlighting the pressing need to develop methods with improved performance for membrane proteins.

In light of personalized treatment strategies for LQTS, the KCNQ1 Rosetta models will be broadly useful for interpreting the functional phenotype of KCNQ1 variants that we demonstrated here by predicting the energetic change in protein stability upon mutation compared to WT KCNQ1. Of the 44 single-site missense mutations studied here, 25 were LOF mutations (channel current <65%WT). Huang et al. [19] found 18 of these LOF mutations exhibited lower (<65%WT) cell surface expression than WT and 13 were severely folding-destabilized as indicated by their NMR spectrum or the inability of the VSD to express in E. coli. For these subsets of variants, we found the median ΔΔG was 17.7 REU and 24.2 REU, respectively, which translates to 13 of 18 and 11 of 13 variants having largely increased ΔΔG values above 10 REU (Fig 7C). In contrast, the median ΔΔG of the 19 non-LOF variants was significantly lower (0.7 REU, p<0.001), and for the 26 variants with WT-like cell surface expression the ΔΔG was 1.1 REU (p<0.001). Thus, the Rosetta ΔΔG predictions support the conclusion that the KCNQ1 LOF variants were energetically destabilized and that protein destabilization represents the likely cause for protein misfolding and low cell surface expression. We previously noticed [19] that many class IV and V variants, which are folding-destabilized, have mutation sites located in S0 and S0-contacting regions (see S7 Fig). The Rosetta models and MD simulations revealed a network of interactions of residues in S0 with residues in S2 and S4 that appear to glue the VSD helix bundle together. This observation illuminates the role of S0 as an important scaffolding element in KCNQ1 and potentially also in other K_V channels.

Energetic destabilization and trafficking defects is a common disease mechanism associated with gene mutations that impact membrane proteins [95]. Examples include the HERG potassium channel and cardiac arrhythmia [96], the human peripheral membrane myelin protein 22 (PMP22) and Charcot-Marie Tooth disease [97], and a subset of mutations in the cystic fibrosis transmembrane conductance regulator (CFTR) protein [98]. It is likely that mutation-induced destabilization is also a trigger for the disease-causative intracellular retention and/or degradation of other membrane proteins e.g. the myelin proteolipid protein that causes Pelizaeus-Merzbacher disease [99].

Noticeable exceptions of variants for which Rosetta-predicted changes in protein thermodynamic stability were not correlated with KCNQ1 loss -⁠ or gain-of-function were variants in class I (e.g. C122Y) as well as A102S (WT-like) and R109L (class III). C122Y is a LOF variant and is associated with LQTS [100]. For this case, the reason for LOF is unknown. In our KCNQ1 models, C122 is close to the putative KCNE1 binding site in KCNQ1 between S1 and S5’ of two neighboring subunits. One may speculate mutation-induced changes of the local structure around C122 could affect the presentation of side-chains on S1 towards KCNE1 and therefore affect KCNE1 binding. Alternatively, the C122Y mutation could affect S4 movement through an allosteric mechanism. Five other mutations in class I (with normal or higher expression and trafficking levels but low or no conductance) are located in regions that are critical for KCNE1 binding and channel gating: on S1 (L134P), S4 (I227L, Q234P) and on S4-contacting regions (A150, K218). L134 is lipid-exposed and its side-chain is pointed towards the putative KCNE1 binding site. I227 is one position before the first gating charge R1 while Q234 occupies the third gating charge position (R3 in other channels, Q3 in KCNQ1). A150 is located in the S1-S2 loop and K218 in the S3-S4 loop, both of which undergo varying interactions with S4 during VSD activation. Thus, the unusual behavior of these class I variants may be related to the location-specific role these residues play in voltage-dependent channel activation mechanism.

The reason for increased conductance of A102S and R109L is unknown. Both sites are located in the first half of S0 and are lipid exposed. A structural mechanistic explanation for them is currently hampered by the lack of structural data for the N-terminal domain before S0, which is not part of the Rosetta models.

Knowledge of KCNQ1 LOF mechanisms pertaining to a given patient may soon facilitate the design of personalized treatment strategies for LQTS. The observation that many LOF variants are energetically destabilized and prone to protein misfolding motivates the idea of developing targeted therapies that assist in folding of destabilized KCNQ1 variants and/or rescue trafficking-defective variants to reach the plasma membrane. The development of such tailored therapies for varying classes of channel defects has become reality for the cystic fibrosis transmembrane regulator (CFTR) protein. Many patient-associated CFTR mutations are trafficking-defective, some of which have been shown to be rescuable by a pharmacological approach [95, 101]. The structural models developed in this study provide a step toward realizing such tailored pharmacological approaches for KCNQ1.

Conclusion

The structural models of the human KCNQ1 channel in closed and open states were developed by comparative modeling based on the high-grade template structure from X. leavis KCNQ1 and other experimental template structures that together provided nearly complete coverage of all channel regions (S0-S6). The KCNQ1 models are experimentally validated and offer improved completeness and accuracy over earlier KCNQ1 Rosetta models. MD simulations demonstrated the capability of the models to describe experimentally established channel properties such as state-dependent VSD gating charge interactions and pore conformations, PIP2 binding, and VSD-PD coupling interactions. Rosetta energy calculations on KCNQ1 models successfully predicted energetic destabilizations in folding-defective KCNQ1 LOF mutants and distinguished them from mutants with normal phenotype. The results in this work highlight the utility of the KCNQ1 structural models in addressing open questions about KCNQ1 structure-function mechanisms and assisting in the diagnosis of disease-causative KCNQ1 mutations towards the development of new therapeutic approaches.

Materials and methods

Homology modeling of KCNQ1 RC state and AO state structures

Structural models of human KCNQ1 (helix S0-S6, residues 100–369) in the RC, AC and AO state were generated by homology modeling using the protein structure prediction software Rosetta (version 3.8) [102]. The X. leavis KCNQ1 cryo-EM structure (PDB 5VMS) [31] provided a template for the AC state as well as separate templates for the VSD in the AO model and for the pore domain in the RC model. Additional templates for the open pore conformation in AO state KCNQ1 were derived from the crystal structure of rat K_V1.2–2.1 channel (PDB 2R9R) [35]. Templates for the VSD in RC state KCNQ1 were obtained from related VSD structures: the resting VSD in C. intestinalis voltage-sensing phosphatase (Ci-VSP) (PDB 4G7Y) [32], VSD2 in A. thaliana two pore calcium channel protein 1 (TPC1) (PDB 5DQQ) [33], and the resting VSD conformation C3 in a model of the Shaker K channel [36]. Hybridization and regularization of multiple templates was accomplished with the Rosetta comparative modeling (RosettaCM) protocol [103] guided by the RosettaMembrane energy function [104, 105]. Missing loop regions and connecting residues between the VSD and pore domain were modeled de novo through fragment insertion. Template threading in RosettaCM was guided by multiple sequence alignments created with ClustalW [106] (see S1 and S2 Figs). Fragment libraries were created with the Rosetta fragment picker [107] incorporating PSIPRED [108] secondary structure prediction. Transmembrane helix regions were predicted by OCTOPUS [109] and used to impose membrane-specific Rosetta energy terms on residues within the theoretical membrane bilayer. During homology modeling, C4 symmetry of the KCNQ1 tetramer was enforced by means of Rosetta symmetry definition files as described in DiMaio et al. [110]. Amino acid side chain positions were optimized by a simulated annealing protocol, referred to as rotamer packing in Rosetta, and models were refined in internal and Cartesian coordinate space by gradient-based minimization while applying a low harmonic restraint to the initial Cα-atom coordinates. Models were clustered based on Cα-atom RMSD, and the top-scoring models of the ten largest clusters as well as another 10–20 models with low Rosetta score and valid helix and loop conformations were selected as final ensemble. Model quality was additionally checked by MolProbity [38] analysis (see supporting S1 Table), and the model with the lowest Molprobity score was considered the final model.

MD simulations of KCNQ1 RC and AO models

MD simulations of KCNQ1 RC and AO models were performed in explicit phospholipid membranes at 310 K using Amber16 [48] and the Lipid17 force field (Gould, I.R., Skjevik A.A., Dickson, C.J., Madej, B.D., Walker, R.C.,"Lipid17: A Comprehensive AMBER Force Field for the Simulation of Zwitterionic and Anionic Lipids", 2018, manuscript in preparation). Representative starting conformations for MD were selected from the ensemble of KCNQ1 homology models and by hierarchical full-linkage clustering.

Models were aligned to the membrane normal using the PPM webserver [111] and embedded into bilayers of POPC (palmitoyloleoyl-phosphatidylcholine) and PIP2 (phosphatidyl-4,5-bisphosphate) (~280 lipids per leaflet) using the membrane builder tool of the CHARMM-GUI website [47]. A TIP3P water layer with 22.5 Å thickness containing 150 mM of KCl was added on either side of the membrane. In addition, four K+ ions were placed in the channel selectivity filter at positions inferred from the X-ray structure coordinates of PDB 2R9R. Bilayers contained 10 mol% of PIP2 in the inner leaflet which comprised equal numbers of C4-PO₄ -⁠ and C5-PO₄-mono-protonated PIP2 molecules with stearoyl and arachidonoyl conjugations at the sn-1 and sn-2 position. The geometry of PIP2 was optimized with Gaussian 09 (Gaussian, Inc., Wallingford CT) on the B3LYP/6-31G* level of theory, and assignment of Amber atom types and calculation of RESP charges was done with Antechamber [112]. Bond and angle parameters of the protonated C4-PO₄ or C5-PO₄ group in PIP2 were adjusted to values previously reported for phosphorylated amino acids [113] to avoid simulation instabilities. No PIP2 was added to the MD system of the AC model because the decoupled state of KCNQ1 in the cryo-EM structure is assumed to be due to the absence of PIP2 and no conformational change should be induced by the membrane. SHAKE [114] bond length constraints were applied to all bonds involving hydrogen. Nonbonded interactions were evaluated with a 10 Å cutoff, and electrostatic interactions were calculated by the particle-mesh Ewald method [115].

Each MD system was first minimized for 15,000 steps using steepest descent followed by 15,000 steps of conjugate gradient minimization. With protein and ions restrained to their initial coordinates, the lipid and water were heated to 50 K over 1000 steps with a step size of 1 fs in the NVT ensemble using Langevin dynamics with a rapid collision frequency of 10,000 ps^-1. The system was then heated to 100 K over 50,000 steps with a collision frequency of 1000 ps^-1 and finally to 310 K over 200,000 steps and a collision frequency of 100 ps^-1. After changing to the NPT ensemble, restraints on ions were gradually removed over 500 ps and the system was equilibrated for another 5 ns at 310 K with weak positional restraints (with a force constant of 1 kcal mol^-1 Å^-2) applied to protein Cα atoms. The protein restraints were then gradually removed over 10 ns, and production MD was conducted for 350 ns (AC model) or 400 ns (RC and AO model) using a step size of 2 fs, constant pressure periodic boundary conditions, anisotropic pressure scaling and Langevin dynamics. Four independent simulations were carried out for each channel model yielding 1.4–1.6 μs of MD data for each state.

Analysis of KCNQ1 MD simulations

Analysis of MD trajectories with CPPTRAJ (version 18.0) [116] included calculation of Cα-atom root-mean-square deviations (Cα-RMSD) and root-mean-square fluctuations (Cα-RMSF), enumeration of protein-protein and protein-lipid hydrogen bonds and KCNQ1 residue contacts as well as principal component analysis. Residue contact numbers were calculated by counting the number of heteroatoms within a 4 Å distance of any heteroatom in a given KCNQ1 residue averaged over the last 300 ns of production MD as well as over all four RC and AO state MD simulations, respectively. For calculating intra-subunit contact numbers, only residues within the same chain were used whereas inter-subunit contact numbers were calculated considering only residues located on two different, directly adjacent chains. Intra -⁠ and inter-subunit contacts numbers were averaged over all four chains (i.e. A–D) or pairs of adjacent chains (i.e. A+B, B+C, C+D, A+D) in the KCNQ1 tetramer, respectively. Pairs of diagonally opposite chains (i.e. A+C and B+D) were not included in the analysis because they make no contacts between the VSD and PD. Furthermore, only long-range contacts (|i-j| > 5 residues) were used in the analysis of intra-subunit interactions.

Measurement of the channel pore radius was carried out with the HOLE program [117] and using snapshots of KCNQ1 taken at 1 ns intervals during the last 300 ns of MD.

Stability calculations of KCNQ1 VSD variants

Free energy changes in 50 functionally characterized KCNQ1 variants [19] with mutations in the VSD were modeled using the Rosetta Flex ddG protocol [74] and the RosettaMembrane all-atom energy function [105]. In short, the Flex ddG protocol models mutation-induced conformational and energetic changes by a series of “backrub” moves (15,000 steps in this study) of the protein backbone together with side-chain repacking within an 8 Å shell around the mutation site, followed by global minimization of all protein backbone and side-chain torsion angles. Minimization is performed with harmonic Cα-atom pair distance restraints to avoid any large structural deviations from the input model. A mutation was modeled in each of the four channel subunits simultaneously by enforcing C4 symmetry by means of Rosetta symmetry definition files. 50 independent trajectories were carried out for each mutant model and the WT model, and the Rosetta energy change (ΔΔG) was calculated as average score difference between the three top-scoring mutant and WT models. Stability calculations were repeated on ten different input models of the RC as well as AO state taken from the selected ensemble of KCNQ1 homology models to calculate the average ΔΔG value (± S.E.M.) for each KCNQ1 variant.

Supporting information

S1 Table [docx]
Molprobity statistic of KCNQ1 AO, RC and AC homology models.

S2 Table [docx]
Comparison of Rosetta predicted stability changes (ΔΔG) with expression levels, trafficking efficiencies and peak current densities of KCNQ1 VSD mutants.

S1 Fig [tpc1]
Multiple sequence alignment of human KCNQ1 for RC homology modeling.

S2 Fig [tiff]
Multiple sequence alignment of human KCNQ1 for AO homology modeling.

S3 Fig [tiff]
Pore radius of KCNQ1 in MD simulations of RC, AC and AO models.

S4 Fig [ac]
Backbone Cα-RMSD of KCNQ1 RC, AC and AO state MD trajectories.

S5 Fig [a]
Global domain motions in KCNQ1 channel models revealed by residue cross-correlation and principal component analysis.

S6 Fig [a]
Conformational changes in KCNQ1 during transition from the AC to AO state inferred from modeling.

S7 Fig [tiff]
Location of mutation sites in the KCNQ1 VSD by variant class.

S8 Fig [a]
Predicted stability changes of KCNQ1 VSD variants calculated with the RC model.

S9 Fig [tiff]
Comparison of Rosetta ΔΔG values (mean ± S.E.M.) for KCNQ1 variants calculated with the RC and AO model.

S10 Fig [a]
Comparison of Rosetta RC and AO models of selected KCNQ1 mutants.

S11 Fig [a]
Structural models of selected KCNQ1 variants.

S12 Fig [a]
Experimentally restrained docking of KCNE1 to KCNQ1 closed and open channel models.

S1 PDB [pdb]
Coordinates of the Rosetta KCNQ1 AO model in PDB format.

S2 PDB [pdb]
Coordinates of the Rosetta KCNQ1 RC model in PDB format.

S3 PDB [pdb]
Coordinates of the Rosetta KCNQ1 AC model in PDB format.

S1 Movie [mp4]
Morph between the Rosetta KCNQ1 RC and AO channel models as viewed parallel to the membrane plane.

S2 Movie [mp4]
Morph between the Rosetta KCNQ1 RC and AO channel models as viewed from the extracellular side.

S3 Movie [mp4]
Morph between the Rosetta KCNQ1 RC and AO channel models as viewed from the cytosolic side.

S4 Movie [mp4]
Morph between the closed and open pore domain of KCNQ1.

S5 Movie [mp4]
Morph between the resting and activated VSD of KCNQ1.

S1 Input Files and Scripts []
This directory contains input data and scripts that were used to create the results in this manuscript.

S1 MDtraj [nc]
Snapshots from 400ns MD trajectory of KCNQ1 AO model in Amber netcdf format.

S2 MDtraj [nc]
Snapshots from 400ns MD trajectory of KCNQ1 RC model in Amber netcdf format.

S3 MDtraj [nc]
Snapshots from 350ns MD trajectory of KCNQ1 AC model in Amber netcdf format.