Simulating PIP2-Induced Gating Transitions in Kir6.2 Channels.

ATP-sensitive potassium (KATP) channels consist of an inwardly rectifying K+ channel (Kir6.2) pore, to which four ATP-sensitive sulfonylurea receptor (SUR) domains are attached, thereby coupling K+ permeation directly to the metabolic state of the cell. Dysfunction is linked to neonatal diabetes and other diseases. K+ flux through these channels is controlled by conformational changes in the helix bundle region, which acts as a physical barrier for K+ permeation. In addition, the G-loop, located in the cytoplasmic domain, and the selectivity filter might contribute to gating, as suggested by different disease-causing mutations. Gating of Kir channels is regulated by different ligands, like Gβγ, H+, Na+, adenosine nucleotides, and the signaling lipid phosphatidyl-inositol 4,5-bisphosphate (PIP2), which is an essential activator for all eukaryotic Kir family members. Although molecular determinants of PIP2 activation of KATP channels have been investigated in functional studies, structural information of the binding site is still lacking as PIP2 could not be resolved in Kir6.2 cryo-EM structures. In this study, we used Molecular Dynamics (MD) simulations to examine the dynamics of residues associated with gating in Kir6.2. By combining this structural information with functional data, we investigated the mechanism underlying Kir6.2 channel regulation by PIP2.

Introduction

Inwardly rectifying K+ (Kir) channels are expressed in diverse tissues and regulate physiological processes by setting the cellular resting membrane potential. K+ efflux is reduced to different degrees due to block by intracellular Mg2+ and polyamines at potentials positive to the K+ equilibrium potential (Ho et al., 1993; Kubo et al., 1993; Nichols and Lee 2018). X-ray and cryo-EM structures of several different Kir family members are available, revealing a remarkably conserved pore architecture, despite widely different ligand gating mechanisms. While all Kir channels require phosphatidylinositol-4,5 bisphosphate (PIP2) binding for channel activation (Hilgemann and Ball 1996; D’Avanzo et al., 2010; Hansen et al., 2011), gating by many additional ligands is unique. For example, gating of Kir1 and Kir4/5 channels is controlled by pH, Kir3 channels are regulated by Gβγ proteins, and Kir6 by ADP/ATP and sulfonylurea receptor subunits (Nichols and Lopatin 1997).

Recent cryo-EM structures of KATP channels (Martin et al., 2017a; Martin et al., 2017b; Lee et al., 2017; Li et al., 2017; Wu et al., 2018; Ding et al., 2019) provide important progress towards understanding the complex gating regulation of this important subfamily, which couples the metabolic state of a cell to its electrical excitability (Hibino et al., 2010; Rorsman and Ashcroft 2018). Structures confirm previous expectations that KATP channels consist of four Kir6.x pore-forming subunits and four regulatory sulfonylurea receptor (SUR) subunits. Like in Kir2 and Kir3 channels, the pore below the selectivity filter (SF) is lined by two main constrictions: the so-called helix bundle crossing (HBC) gate and the G-loop gate. K+ flux through KATP channels is inhibited by direct interactions of the cytoplasmic Kir6 domain with ATP, while Mg-nucleotide binding to SUR also modulates these channels. Furthermore, channel opening is potentiated by PIP2 binding to Kir6, which reduces channel inhibition by ATP (Baukrowitz et al., 1998; Shyng and Nichols 1998). The molecular determinants of PIP2 activation in KATP channels have been investigated in functional studies (Fan and Makielski 1997; Shyng et al., 2000; Cukras et al., 2002; Schulze et al., 2003; Haider et al., 2007; Männikkö et al., 2011; Pipatpolkai et al., 2020). However, structural information of the binding site and the gating transitions leading to channel opening are still missing, since all available cryo-EM structures have been solved in the absence of this activator, even though the molecule was present in some cryo-EM experiments (Lee et al., 2017; Wu et al., 2018) (see Table 1 for details). This significantly limits our understanding of the molecular mechanisms by which PIP2 and other ligands regulate these channels.

TABLE 1

Available Kir6.2 cryo-EM structures in the PDB compared to a pre-open Kir3.2 x-ray structure.

PDB	State	Experimental construct	Ligands in model	Ligands in experiment	Resolution (Å)	r(L164) (Å)	r(F168) (Å)	r(I296) (Å)	Reference
5TWV	Closed	Hamster SUR1, rat Kir6.2	ATP	Glibenclamide, ATP	6.30	1.27	1.66	2.65	Martin et al. (2017b)
6C3P	Closed (propeller)	Human SUR1, human Kir6.2 fusion construct	ATP, ADP, Mg2+	ATP, ADP, Mg2+, diC8-PIP2	5.60	0.97	1.24	2.36	Lee et al. (2017)
6C3O	Closed (quatrefoil)	Human SUR1, human Kir6.2 fusion construct	ATP, ADP, K+, Mg2+	ATP, ADP, Mg2+, diC8-PIP2	3.90	0.82	1.06	1.59	Lee et al. (2017)
5YW8	Closed	Hamster SUR1, mouse Kir6.2 fusion construct	ATPγS	ATPγS, glibenclamide	4.40	0.84	0.76	2.7	Wu et al. (2018)
5YW9	Closed, (T state, tense)	Hamster SUR1, mouse Kir6.2 fusion construct	ATPγS	ATPγS	5.00	0.84	0.76	2.7	Wu et al. (2018)
5YWC	Closed	Hamster SUR1, mouse Kir6.2 fusion construct	ADP, Mg2+	Mg-ADP, VO4 3-, diC8-PIP2, NN414 (a KCO)	4.30	0.87	0.74	2.64	Wu et al. (2018)
5YWA	Closed, (R state, relaxed)	Hamster SUR1, mouse Kir6.2 fusion construct	ATPγS	ATPγS	6.10	0.89	0.71	2.91	Wu et al. (2018)
5YKF	Closed, (T state, tense)	Hamster SUR1, mouse Kir6.2 fusion construct	ATPγS, glibenclamide	ATPγS, glibenclamide	4.33	0.82	0.69	2.58	Wu et al. (2018)
5YWB	Closed	Hamster SUR1, mouse Kir6.2 fusion construct	ADP, Mg2+	Mg-ADP, VO4 3-, diC8-PIP2, NN414 (a KCO)	5.20	0.99	0.62	2.7	Wu et al. (2018)
5YKG	Closed, (R state, relaxed)	Hamster SUR1, mouse Kir6.2 fusion construct	ATPγS, glibenclamide	ATPγS, glibenclamide	4.57	0.86	0.59	2.8	Wu et al. (2018)
5YKE	Closed, focus refined TMD + SUR, no CTD	Hamster SUR1, mouse Kir6.2 fusion construct	Glibenclamide	ATPγS, glibenclamide	4.11	0.81	0.5	N/A	Wu et al. (2018)
6BAA	Closed	Hamster SUR1, rat Kir6.2	Glibenclamide, ATP	Glibenclamide, ATP	3.63	1.1	0.46	2.21	(Martinet al., 2017a)
6JB1	Closed, (T state, tense)	Hamster SUR1 (Y1209S), mouse Kir6.2 fusion construct	ATPγS, repaglinide, lipids (POPC, PE), digitonin	ATPγS, repaglinide, digitonin, KNtp peptide	3.30	0.96	0.38	2.41	Ding et al. (2019)
5WUA	Closed	Hamster SUR1 (Q608K), mouse Kir6.2	−	Glibenclamide	5.60	1.05	Too narrow	1.5	Li et al. (2017)
3SYA (Kir3.2 ref)	Closed, pre-open	Mouse Kir3.2 crystal structure	diC8-PIP2, K+, Na+	diC8-PIP2, Na+	2.98	2.95	2.45	2.96	Whorton and MacKinnon (2011)

The red-white-green color scale of the Resolution column ranges from green (higher resolutions) to red (lower resolutions). Additionally, among the columns for the three pore constriction sites r(L164), r(F168), and r(I296), the colors range from red (narrow pore radii) to green (wider radii).

Functional studies and modeling strongly support a common structural basis for PIP2 regulation in different Kir family members (Fan and Makielski 1997; Shyng et al., 2000; Cukras et al., 2002; Schulze et al., 2003; Haider et al., 2007; Männikkö et al., 2011; Fürst et al., 2014; Pipatpolkai et al., 2020). Since PIP2 binding sites are well resolved in Kir2 (Hansen et al., 2011; Lee S. J. et al., 2016; Zangerl-Plessl et al., 2019) and Kir3 channels (Whorton and MacKinnon 2011; Whorton and MacKinnon 2013; Niu et al., 2020), we used this structural information, together with functional data, to investigate the possible structural basis for KATP channel regulation by PIP2. Specifically, the aim of the present study was to examine PIP2-induced gating transitions (after unbinding of ATP from the cytoplasmic domain) using Molecular Dynamics (MD) simulations. Furthermore, we investigated structural changes of the Permanent Neonatal Diabetes mutation (PNDM) L164P that dramatically alters the open state stability of the channel.

Results and Discussion

Currently Available KATP Structures Contain Three Constriction Sites

Fourteen cryo-EM structures were resolved to date, with resolutions ranging from 3.3 to 6.3 Å. A complete list of currently available Kir6.2 models with included ligands is given in Table 1. To obtain a comprehensive overview of pore dimensions and possible constriction sites, we analyzed pore profiles using the program HOLE (Smart et al., 1996). As shown in Figure 1 and Table 1, three main constriction sites could be identified. Constriction site 1 is located at L164, one helix-turn above the so-called helix bundle crossing (HBC). Constriction site 2, formed by F168 side chains, frames the canonical HBC gate, while I296, a residue associated with the G-loop gate, constitutes the third constriction site. The pore diameter at constriction site one ranges from 0.81 to 1.27 Å, in the respective cryo-EM structures. This is in line with previous studies which identified L164 as a narrow, pore-lining site in Kir6.2 (Loussouarn et al., 2000; Loussouarn et al., 2001; Kurata et al., 2004; Proks et al., 2005; Kurata et al., 2006; Walczewska-Szewc and Nowak 2020). Distances at constriction site 2 range from 0.81 to 1.27 Å, while the pore was slightly wider at I296 (constriction site 3), with distances ranging from 1.5 to 2.9 Å. Overall, none of the structures showed pore radii large enough to enable hydrated K+ flux through a continuous pore. Even though recent studies suggest that partially dehydrated K+ ions can pass the HBC gate formed by aromatic side chains (Bernsteiner et al., 2019; Black et al., 2020), the constrictions at sites 1 and 3 are formed by hydrophobic residues, rendering K+ passage (especially at site 1) very unlikely.

FIGURE 1

Comparison of the pore dimensions of Kir6.2 structures with Kir3.2. Slices through the pore-forming Kir channel surfaces show three Kir models in the closed state, whereas the pore of Kir3.2 (3SYA) has been described as ‘pre-open’. Kir surfaces are colored according to the electrostatic potential. On the right side of each Kir channel, a HOLE profile shows the pore diameter. The three narrowest constriction sites in Kir6.2 below the SF are annotated. In the HOLE profile, red indicates a pore radius too tight for a water molecule. Green suggests space for a single water molecule (0.6 Å < r < 1.15 Å), while blue indicates at least twice the radius of a single water molecule.

Dynamics of Constriction Sites in MD Simulations

Given our recent success in using MD simulations to provide functional interpretation of conductive states (Bernsteiner et al., 2019; Lee S. J. et al., 2016; Zangerl-Plessl et al., 2019), or mutant-induced gating transitions (Linder et al., 2015), we chose to simulate the pore domain of two Kir6.2 structures. Therefore, we selected PDB accession no. 6C3O and 6C3P (highest-resolution structure pair, solved under identical conditions), acquired in the presence of the Kir6.2 activators Mg-ADP and PIP2 during the experimental setup (Lee et al., 2017) which have quite different radii at the constriction sites as shown in Figure 1. Our recent MD simulations on a Kir3 structure (PDB accession no. 3SYA) revealed that spontaneous wetting of hydrophobic gates could be obtained in classical, atomistic simulations within relatively short time scales (Bernsteiner et al., 2019). Thus, we performed similar MD simulations with the ATP-unbound, closed pore domains of these Kir6.2 structures after placing short-chain PIP2 molecules in the respective binding sites, using the 3SYA structure as a template. For details about the setup see Supplementary Figure S1. Both, simulating Kir6.2 in the absence of SUR1 (Pipatpolkai et al., 2020), as well as cross-talk analysis between ATP and PIP2, revealed that removal of SUR1 or ATP (Pipatpolkai et al., 2021) from the cryo-EM structure have a marginal effect on the PIP2 binding site. Specifically, Pipatpolkai et al. had shown that the presence of ATP primarily influences the dynamics of residue K39 (Pipatpolkai et al., 2021). Comparison of the root-mean-square fluctuations (RMSF) of the PIP2 binding residues of our simulations (Supplementary Figure S2) with those from Pipatpolkai et al. revealed a similar behavior: Simulation setups containing only PIP2 versus setups with PIP2 plus ATP suggest no critical bias of our simulations caused by the removal of ATP.

To differentiate the effect induced by PIP2 from inherent protein dynamics, we compared PIP2-bound simulations to apo control runs. For simplicity, we will hereinafter refer to these simulations systems as “6C3O” and “6C3P” for the PIP2-bound WT simulations, and “6C3O apo” and “6C3P apo” for the control runs without PIP2. See Table 2 for an overview of simulations.

TABLE 2

Overview of MD simulations.

Name	n	PIP2	Length (ns)	Electric field (mV/nm)
6C3O	5	+	1,000	−
6C3O apo	5	−	200	−
6C3P	5	+	1,000	−
6C3P apo	5	−	200	−
6C3P L164P	1	+	1,000	−
6C3P L164P	1	+	1,000	40 (outward)

The stability of the different protein systems, measured as the root-mean-square deviation (RMSD), is shown in Supplementary Figure S3. As can be seen in the RMSD plots, all systems reached stable values below 4 Å, which agrees well with previous simulations and a similar system setup (Pipatpolkai et al., 2021).

To assess the pore’s lumen at the constriction sites, we monitored the minimum distance of each residue between opposing subunits of 5 × 1 µs simulations over time, as shown in Figures 2A–C. L164, forming constriction site 1, sampled minimum distances between fully closed and slightly widened apertures, with average minimum distances of 3.4 and 4.4 Å in 6C3O and 6C3P, respectively. Still, both systems remained not only too narrow for hydrated K+ ions (6.6 Å (Conway 1981)), but also for water molecules to pass. The pore diameters at the canonical HBC gate were wider in both systems. As shown in Figure 2B, a bimodal distribution at F168 caused an asymmetrical pore geometry in both systems. Minimum distances between two opposing subunits frequently sampled 2.6 and 5.8 Å in the 6C3P structure, while the wider pair exhibited distances around 6.7 Å in 6C3O. Despite starting from a narrower HBC gate (2.12 Å diameter), pore diameters increased up to 10 Å in 6C3O, leading to sporadic wetting of the lower HBC gate. However, no persistent wetting in this region could be observed within 1 µs simulation time, as shown in Supplementary Movies S1, S2. Summarizing, PIP2 influences the gate diameter of the HBC gate. This is particularly evident for the 6C3O structure, which occasionally permitted distances wider than the first K+ solvation shell. The G-loop gate diameter increased in both structures, peaking at 9.3 Å in 6C3O and 10.2 Å in 6C3P, which led to continuous solvation in all simulations (Figure 2C, and Supplementary Movies S1, S2). Interestingly, no changes could be observed between PIP2-holo and apo simulations at constriction site 1 (Supplementary Figure S4), whereas the bimodal distributions of both systems were shifted towards narrower apertures without bound PIP2 in constriction sites 2 and 3.

FIGURE 2

Histograms of pore diameters, sampled in two Kir6.2 systems (6C3O, 6C3P). Minimum distances for three major constriction sites in Kir6.2 were measured in 5 × 1 μs MD simulations for both systems between two opposing subunits and averaged over the number of simulations, subsequently. Crosses mark the corresponding distances of the initial cryo-EM state before equilibration and production run, measured with the HOLE program. A vertical line is drawn at 6.6 Å, indicating the time-averaged hydration diameter of K+ (Conway 1981). Analogous average minimum-distance values, measured in the G-loop gate (narrowest passage in 3SYA), were associated with conduction in previously published Kir3.2 3SYA simulations (Bernsteiner et al., 2019) (*, gray curve).

Monitoring PIP2-Induced Gating Transitions

Given the different behavior of the PIP2-bound structures in the simulations, we analyzed the structural changes that led to the short-lived wetting and opening transitions at the HBC gate observed in the 6C3O structure. As previously described in structural and functional studies (Hansen et al., 2011; Fürst et al., 2014; Li et al., 2015; Poveda et al., 2017; Niu et al., 2020), we observed a key structural change in the C-linker, characterized by a PIP2-dependent conversion of the loop into a helix. Thus, we monitored changes in the secondary structure of the C-linker over simulation time in both systems, as shown in Figures 3C,E, and subplots A, C, E, G of Supplementary Figures S5–S7. Interestingly, there is a ∼30% higher helical content (310 helix, corresponding to gray lines in the qualitative plots) in the 6C3O C-linker, compared to 6C3P, suggesting a possible connection between the conformation of the C-linker region and the HBC gate. What might cause this early PIP2-dependent step towards opening? When comparing residues that render the PIP2-binding site in 6C3O and 6C3P, a critical difference becomes apparent. R176 was determined to be a key residue for PIP2-dependent activation in Kir6.2 (Fan and Makielski 1997; Baukrowitz et al., 1998; Shyng et al., 2000; Haider et al., 2007; Pipatpolkai et al., 2020; Pipatpolkai et al., 2021) and other Kir families (Soom et al., 2001; Lopes et al., 2002; Xie et al., 2005; Lee C. H. et al., 2016; Lacin et al., 2017). As shown in Figures 3D,F and Supplementary Figures S5, S6, S8, R176 formed stable h-bonds with PIP2 in all five replicas of the 6C3O system, while h-bonds were missing in one to three subunits of 6C3P simulations. Quantification of PIP2 occupancies of the two Kir6.2 systems revealed a significant difference for this residue: While R176 was in close distance and formed h-bonds with the PIP2 head group almost 100% of the time in the 6C3O system, the same interaction was observed only ∼50% of the time in the 6C3P system (Supplementary Figure S8).

FIGURE 3

PIP2 binding site and PIP2-induced gating changes. Snapshots of polar Kir6.2-PIP2 interactions after 1 µs of a representative (A) 6C3O and (B) 6C3P simulation. Relevant residues are shown as sticks, where colors correspond to a protein subunit (SU). The C-linker is colored in teal. (C,E) Plots correspond to four Kir subunits of a representative simulation system. The secondary structure of the (C) 6C3O and (E) 6C3P C-linkers is shown over time, where the y-axis ranges between residues Q173 and T180. The color legend is at the bottom of the figure. (D,F) The number of hydrogen bonds is shown between PIP2 molecules and Kir6.2 residues over time for (D) 6C3O and (F) 6C3P. Subunits in (C,D) and (E,F) correspond.

Analysis of the helical content of apo vs. holo runs (Supplementary Figure S7) showed a trend for a higher helical content of the C-linker region between holo vs. apo systems in 6C3O. However, only the difference between 6C3O and 6C3P holo systems was statistically significant. Nevertheless, the helical content of the 6C3O holo systems increased slightly with longer sampling times (5 × 1 µs), albeit not statistically significant. This suggests that much longer simulations will be required to sample a more complete picture of the conformational changes in this region.

In agreement with recent simulations (Pipatpolkai et al., 2021), our analyses revealed that PIP2 binding reduces the dynamics of all binding residues (see Supplementary Figure S2). The effect is statistically significant for most residues in 6C3P, whereas the effect on 6C3O was less pronounced. Furthermore, it has previously been shown that K39 plays a critical role for ATP and particularly for PIP2-dependent interactions. We saw similar effects in our PIP2 simulations, as shown in Figures 3D,F and Supplementary Figures S2, S5, S6, S8. However, in contrast to R176, residue K39 seems to play a more general role in stabilizing PIP2 interactions, independent of HBC gate changes or the helical content of the C-linker, since there is no statistically significant difference of occupancies between the different Kir6.2 systems for this residue. In contrast, significant differences of PIP2 contacts between 6C3O and 6C3P were observed for residues R54, L66, W68, P69, K170, and R716 (Supplementary Figure S8). Combined with the previous observation that a salt bridge between ATP and K39 is broken in the presence of PIP2, which subsequently leads to a rotation of the K39 amine group towards the PIP2 head group (Pipatpolkai et al., 2021), this might suggest the following order of gating transition events: In a first step, PIP2 binding competes with ATP for coordination of the amine group of K39. In a second step, the K39 side chain reorients towards more favorable interactions with PIP2, thereby reducing ATP inhibition (Hilgemann and Ball 1996; Baukrowitz et al., 1998; Shyng and Nichols 1998; Fan and Makielski 1999). Thirdly, after the release of ATP, h-bonds between R176 and PIP2 are gradually strengthened, combined with a conformational change of the C-linker and rotational changes of the CTD (see below for discussion), which in turn lead to widening of the gate regions and subsequent ion flux. It will probably require many further simulations, including enhanced sampling techniques in the future to address the exact order of events and kinetics of these processes.

Another important conformational change, reported in several previous studies, concerns the dynamics of the cytoplasmic domain including a rotation of the whole domain (Clarke et al., 2010; Whorton and MacKinnon 2011; Whorton and MacKinnon 2013; Bavro et al., 2012; Li et al., 2015; Li et al., 2017; Linder et al., 2015; Wu et al., 2018; Zangerl-Plessl et al., 2019). We, therefore, monitored the rotation of the CTD over simulation time, as shown in Supplementary Figure S9. Indeed, we observed a stronger relative clockwise rotation in the 6C3O simulations (9.2° ± 4.3° std. after 1 µs, n = 5, viewed from the extracellular side), compared to 6C3P. Absolute rotation angles of PIP2-bound systems converged at ∼70° (measured between TMD and CTD of a single protein subunit) in both systems. Nevertheless, fluctuations of the rotation angles render the interpretation of this conformational change rather speculative. Future experiments are necessary to address this question in more detail. Furthermore, an important limitation of our study concerns the lack of the SUR domain, which was not included in our simulation systems, due to the large system size, poor resolution, and many missing residues of these domains. It has to be noted, however, that the N-terminal transmembrane domain of SUR1 modulates PIP2 sensitivity in Kir6.2 channels (Pratt et al., 2011; Walczewska-Szewc and Nowak 2020). This might explain why no wetting or widening at the uppermost constriction site was observed in our simulations, rendering the channels essentially non-conductive.

Disease Mutation L164P Strongly Influences Pore Geometry

L164P, a mutation targeting constriction site 1, exhibits very high open-state stability and is insensitive to ATP regulation (Enkvetchakul et al., 2000; Loussouarn et al., 2000; Loussouarn et al., 2001). Heterozygous L164P mutations cause permanent neonatal diabetes by significantly increasing intrinsic open probability (Po) compared to WT channels, while retaining single-channel current amplitudes (Flanagan et al., 2006; Tammaro et al., 2008). To investigate proline-induced structural changes on the pore we performed MD simulations of the mutant, using the 6C3P system after 1 µs simulation time (Figure 4A). As the mutation targets constriction site 1, which stayed closed and desolvated in our previous WT simulations, we expected pronounced effects on this gating region. Indeed, Figures 4B,C and Supplementary Movie S3 already show changes in pore solvation during a 1 µs simulation. Due to the lack of the hydrophobic side chain, the minimum distance at constriction site 1 increased to ∼10 Å, leading to rapid solvation of the cavity above residue F168 (Figure 4D). Figure 4E shows a shift of the bimodal distribution towards wider average minimum distances (∼6.2 Å). Subsequently, constriction site 2 wetted spontaneously with intermittent desolvated periods in both mutant replicas (see Supplementary Movie S3). Furthermore, the G-loop gate stabilized in an open and fully solvated conformation, peaking at ∼10 Å (Figure 4F). To examine whether the effect on the HBC gate was already enough to facilitate ion conduction, we set up a control replica with a 40 mV nm−1 outward driving force. While the simulation revealed a marked effect on the HBC gate, 1 µs sampling time was still not enough to allow continuous K+ permeation through the whole pore, probably requiring longer sampling times in future studies.

FIGURE 4

L164P leads to widening of the gates and pore solvation. Snapshots of (A) a WT 6C3P system and the equilibrated mutant system 6C3P L164P (B) before and (C) after 1 µs simulation show the solvation state of the HBC gate. The G-loop gate stabilized fully in an open conformation. (D–F) As described in Figure 2, histograms of pore diameters show minimum distances for the three major constriction sites of 6C3P and 6C3P L164P. Constriction sites one and three sampled distances wider than the time-averaged hydration diameter of K+ in the mutant system. The bimodal distribution at constriction site two shifts towards wider apertures, leading to intermittent solvation of the HBC gate. Distances were measured in 2 × 1 μs MD simulation for both systems between two opposing subunits and averaged over the number of simulations, subsequently.

Since proline is known to induce kinks in helices (Cordes et al., 2002), we analyzed the maximum helix angle along the length of the M2 helix over time. Supplementary Figure S10 shows a significant increase in the maximum M2 helix angle from ∼10° to more than 15° on average between 6C3P WT and L164P mutant simulations. This kinking of M2 helices and the lack of the hydrophobic side chain at constriction site 1 led to a widening of the pore. Given the fact that previous studies reported that mutating L164 to cysteine, alanine, valine, threonine or glycine also leads to a large increase in Po (Enkvetchakul et al., 2000; Enkvetchakul et al., 2001; Loussouarn et al., 2000; Loussouarn et al., 2001), this suggests that pore geometry and hydrophobicity at these sites are critical parameters for keeping the pore closed, independent of allosteric ligand effects, such as PIP2. Unfortunately, mutations at this site can reduce block by drugs such as verapamil (L164C) (Ninomiya et al., 2003) or chloroquine (L164A) (Ponce-Balbuena et al., 2012). Moreover, patients carrying the L164P mutation could not be converted to sulfonylurea therapy. The fact that mutant channels are insensitive to regulation by these drugs, SUR, and ATP (Tammaro et al., 2008) corroborates the relevance of L164 in Kir6.2 gating. Our simulations suggest a link between increased open probability and changed pore geometry. Nevertheless, future studies, including docking and drug screening, will be necessary to address this issue in more detail.

Conclusion

Our study provides mechanistic insights into how PIP2 influences the narrow constriction sites observed in Kir6.2 cryo-EM structures and unravels first structural changes induced by the permanent neonatal diabetes mutation L164P on an atomistic level.

Based on our simulations and recently published simulations on ATP-PIP2 interactions (Pipatpolkai et al., 2021), we propose the following order of gating events in the Kir6.2 pore region:

In a first step, a salt bridge between ATP and K39 is broken in the presence of PIP2, which subsequently leads to a rotation of the K39 amine group towards the PIP2 head group (Pipatpolkai et al., 2021). This is followed by a gradual strengthening of the h-bonds between R176 and PIP2 upon release of ATP (our current study). In addition, a conformational change of the C-linker and rotational changes of the CTD will lead to widening of the gate regions and will promote subsequent K+ ion flux (not yet seen in simulations).

Future studies are necessary to obtain a more complete picture, particularly including SUR subunits. This will, however, require higher resolution structures of these regions and much longer sampling times which are beyond the scope of the current paper.

Materials and Methods

Pore Analysis With HOLE

Available molecular assemblies of Kir6.2 from the PDB (Research Collaboratory for Structural Bioinformatics Protein Data Bank (RCSB PDB), RRID:SCR_012820) (see Table 1), were aligned with the Swiss-PdbViewer (Swiss-PdbViewerDeepViewv4.0, RRID:SCR_013295) (Guex and Peitsch 1997) at the selectivity filter motif TTIGYG. Hydrogens were added to the pore-only files with the help of the APBS-PDB2PQR web server (Adaptive Poisson-Boltzmann Solver, RRID:SCR_008387) (Baker, 2001; Dolinsky et al., 2004) with standard settings (pH 7.0, PARSE FF, internal naming scheme, neutral N and C termini). Slices through the pore-forming Kir channel surfaces were generated with PyMOL (PyMOL, RRID:SCR_000305) (Schrödinger 2015). The surface is colored according to the APBS map with a grid spacing of 0.5. The pqr files were used to calculate the pore dimensions of the Kir channels with the HOLE program (version 2.0) (Smart et al., 1996). Standard settings were used, with pore ending radii (endrad) of 6 Å in channels. The HOLE pore geometry was analyzed with Matplotlib (MatPlotLib, RRID:SCR_008624) (Hunter 2007) and visualized with VMD (Visual Molecular Dynamics, RRID:SCR_001820) (Humphrey et al., 1996).

MD Simulations

Four short-chain C8-PIP2 molecules were inserted in both Kir6.2 systems (PDB accession no. 6C3O and 6C3P) based on the respective binding conformations, using the 3SYA structure (Whorton and MacKinnon, 2011) as a template. For details about the process see Supplementary Figure S1. PIP2 parameters were taken from our previous work (Lee S. J. et al., 2016). For our MD simulations, we used the amber99sb force field (Hornak et al., 2006) with Berger lipid parameters (Berger et al., 1997). A cut-off of 1.0 nm was employed for Lennard-Jones and electrostatic interactions and long-range interactions calculated with the Particle-Mesh Ewald algorithm (Essmann et al., 1995). The LINCS algorithm (Hess et al., 1997) was used for bond constraints, allowing for an integration time step of 2 fs. V-rescale was used to couple the temperature to 310 K, with a coupling constant of 0.1 ps. The pressure was kept constant semi-isotropically at 1 bar, using the Parrinello-Rahman barostat with τ = 2 ps (Parrinello and Rahman 1981).

Protein topologies were created with the Gromacs module pdb2gmx (Abraham et al., 2018). The PIP2-bound or apo 6C3O structure was placed in the pre-equilibrated POPC (1-palmitoyl-2-oleoyl-sn-glycero-3- phosphocholine) membrane of our Kir3.2 system (Bernsteiner et al., 2019), containing 588 lipids. K+ ions were placed at selectivity filter positions (0,2,4), separated by water molecules. After solvation with 82,131 water molecules (SPC/E) (Berendsen et al., 1987) and neutralization with K+ ions, 150 mM KCl was added to the simulation system. The energy was minimized (steepest descent) in all setups, followed by 1 ns NVT and at least 10 ns NPT equilibration until convergence of temperature and pressure, with restraints on protein atoms to their starting positions (force constant of 1,000 kJ mol−1 nm−2). 6C3P systems were set up in the same way, while the L164P mutation was introduced into Kir6.2 with the Swiss PdbViewer (Guex and Peitsch 1997). We chose 6C3P for the mutation, as the model exhibited a wider pore geometry at the cryo-EM state. Gromacs version 2018.8 (GROMACS, RRID:SCR_014565) (Abraham et al., 2015; Abraham et al., 2018) was used to perform all-atomistic MD simulations (see Table 2 for an overview). In a single 6C3P L164P simulation an electric field of 40 mV nm−1along the z-axis of the simulation box was applied to test K+ conductivity. With a simulation system size of ∼17.3 nm in the z-direction, this resulted in an electric transmembrane potential of ∼700 mV (Treptow and Tarek 2006; Roux 2008). The overall protein stability was assessed by calculating Cα RMSDs with Gromacs after backbone alignment and exclusion of the highly flexible extracellular loop (r97-110) and intracellular C-terminus (r346-359). Figures were rendered with VMD (Humphrey et al., 1996) and PyMol (Humphrey et al., 1996) and movies were produced using Molywood (Wieczór et al., 2020). Maximum helix angles were calculated with the VMD Bendix plugin (Dahl et al., 2012). All data were plotted with Matplotlib (Hunter 2007) and statistical significance was determined using SciPy (Virtanen et al., 2020).

Minimum Distance Analysis

Minimum distances between opposing subunits of pore-forming Kir subunits were calculated with the Gromacs module gmx mindist (Abraham et al., 2018) for residues L164, F168, and I296 in Kir6.2, and corresponding residues in Kir3.2 (V188, F192, M319). The minimum distances of all simulations of a simulation system and both opposing subunits were combined and averaged over the number of simulations.

PIP2 Binding Site Analysis

We screened the MD trajectories with PyLipid (Song et al., 2021) for residues that were in close proximity (single cutoff: 4 Å) to the PIP2 headgroup. Residues that formed frequent h-bonds (cut-off angle 30°, cut-off radius: 0.35 nm) between PIP2 and Kir6.2 were analyzed with the VMD Hydrogen Bonds extension. Selected h-bonds were further analyzed with Gromacs, using the hbond module. The dynamics of residues in the PIP2 binding site were assessed with the Gromacs module rmsf after aligning the protein along the backbone.

Relative TMD-CTD Rotation

The relative rotation of a single subunit’s cytoplasmic domain (CTD) with respect to its transmembrane domain (TMD) was calculated with an in-house script, as described previously (Bernsteiner et al., 2019). The torsional angle was defined by four points: The center of mass (COM) of the TMD’s subunit 1 (point 1), the COM of the whole TMD (point 2), the COM of the whole CTD (point 3), and the COM of CTD’s subunit 1 (point 4). Kir6.2 residues 53–98 and 116–172 were defined as TMD, while residues 32–52 and 173–346 contributed to the CTD. Thus, the highly flexible extracellular loops in the TMD and the C-termini were excluded for the COM calculations (see Supplementary Figure S9E).

Secondary Structure of the C-Linker

The average helical content of the C-linker (= Tether Helix, residues 173–180) was calculated with the Gromacs module do_dssp (Kabsch and Sander 1983; Touw et al., 2015) and was defined as a combination of 310-helices, α-helices, and π-helices for quantitative analysis.