PIP₂-driven cytoplasmic domain motions are coupled to Kir2 channel gating

Inwardly rectifying potassium (Kir) channels play critical roles in controlling membrane potentials of both excitable and non-excitable cells (Nichols and Lopatin, 1997; Hibino et al., 2010). They are implicated in multiple physiological phenomena and underlie multiple pathologies, including cardiac arrhythmias (Nichols et al., 1996), neonatal diabetes (Koster et al., 2000), SeSAME syndrome (Sicca et al., 2011), and Cantu syndrome (Cooper et al., 2014), as well as altered neuronal excitability leading to epilepsy, Down syndrome, and Parkinson’s syndrome (Zhao et al., 2021). Kir channels are formed as homo- or pseudo-homotetrameric proteins, each subunit consisting of a transmembrane pore domain (TMD) and a cytoplasmic domain (CTD) (Fig. 1, a and b). Three “gates” are implicated in control of K⁺ permeation through Kir channels (Fig. 1 b): a selectivity filter (SF) gate at the extracellular side of the TMD, the helix bundle crossing (HBC) gate at the intracellular side of the TMD, and a G-loop gate at the top of the CTD (Fig. 1, b and c). Ligands that bind to the CTD are major regulators of channel activity, controlling the HBC and G-loop constrictions (Hibino et al., 2010). While each Kir subtype is sensitive to different and diverse regulators, the binding of phosphoinositide-4,5-bisphosphate (PIP₂) at the intracellular leaflet is commonly required for the opening of all Kir channels (Hibino et al., 2010) (Fig. 1 b). Therefore, the molecular mechanisms of PIP₂ gating are of keen interest.

High-resolution crystallographic structures have provided exquisite molecular details of the PIP₂-channel interaction, but the consequent structural rearrangements associated with channel opening are less clear. The very first Kir channel structure in complex with PIP₂ determined the primary PIP₂-binding site at a CTD–TMD interface that is generated by vertical upward movement of the CTD into a so-called “compact” conformation (Tao et al., 2009; Hansen et al., 2011). CTD rotations about the central pore axis have also been observed in Kir channel structures and proposed to be important for channel gating (Bavro et al., 2012). However, mechanistic linkage between PIP₂ binding and such conformational changes has also been questioned for two main reasons: first, there are multiple structures in complex with PIP₂ that remain closed and exhibit different degrees of CTD rotation (Clarke et al., 2010; Bavro et al., 2012; Zhao and MacKinnon, 2021), and second, similar conformational changes have been observed in the absence of PIP₂ (Whorton and MacKinnon, 2011; Lee et al., 2016; Niu et al., 2020; Zhao and MacKinnon, 2021).

Conformational dynamics of proteins, encoded in their structure (topology), are tightly linked to protein function (Demirel and Keskin, 2005; Bahar et al., 2010b). Coarse-grained in silico modeling approaches, such as normal mode analysis (NMA) (Brooks and Karplus, 1983; Go et al., 1983) or even simpler elastic network model (ENM) analysis (Tirion, 1996; Bahar et al., 1997; Frappier and Najmanovich, 2014), have successfully been employed to elucidate intrinsic protein motions without the need for expensive all-atom MD simulations (Chennubhotla et al., 2008). These methods have been applied to subdomain structures (Haider et al., 2005), prokaryotic Kir channels (Shrivastava and Bahar, 2006), and recently to a full-length Kir2.1 channel (Fernandes et al., 2022), but in a conformation with the CTD disengaged, leaving unclear what the intrinsic dynamics after adopting the compact conformation are.

We previously utilized the crystallographic structure of a forced-open mutant (K62W/G178D) channel in complex with PIP₂, to assess K⁺ ion conduction through open Kir2.2 channels in atomistic MD simulations (Zangerl-Plessl et al., 2020). Opening of the HBC gate was clear, and K⁺ conductance through the open channels was of a similar order of magnitude to those measured in single channel recordings. While thorough analyses were performed to characterize K⁺ ion movement through the open channels, we did not pursue the details of the structural rearrangements of the CTD that underlie the PIP₂-driven opening and closing of the channel. In the present study, we have reprocessed these simulations and generated new MD trajectories after removal of PIP₂ molecules from the open conductive channel to elucidate PIP₂-dependent conformational changes in the CTD that actually control opening and closing. We identify key intrinsic motions of the channel that are variably affected by the binding of PIP₂ and specific PIP₂-dependent conformational changes that occur after a compact conformation is adopted, leading to a coherent structural basis for PIP₂-driven cooperative gating.

MD simulations and subsequent analyses were performed using the software package Gromacs v.5.1.2 (Abraham and Gready, 2011) and MATLAB, using in-house scripts and built-in functions. A Kir2.2 K62W/G178D(G) (KW/GD(G)) + PIP₂ 1 μs (“with PIP₂”) simulation was produced in our previous study (Zangerl-Plessl et al., 2020). Four additional replica simulations were additionally performed. The initial structure of KW/GD(G) “without PIP₂” simulations was obtained by deleting 4 PIP₂ molecules from the snapshot of KW/GD(G) + PIP₂ simulation at 50 ns, and 15 replica simulations were performed (Table 1).

Berger lipid parameters were employed for palmitoyloleoylphosphatidylcholine lipids (Berger et al., 1997; Cordomí et al., 2012), the sole constituent of the bilayer, which was surrounded by explicit water molecules (SPCE water model (Berendsen et al., 1987)) and 0.2 M potassium chloride (Joung and Cheatham, 2008). The amber99sb force field was used for the protein (Hornak et al., 2006). For the ions, we used corrected monovalent Lennard–Jones parameters for the amber force field (Joung and Cheatham, 2008) with cutoff set to 1.0 nm. The cutoff for electrostatic interactions was set to 1.0 nm; long-range electrostatic interactions were treated by the particle-mesh Ewald method (Darden et al., 1993). The LINCS algorithm was used to constrain bonds (Hess et al., 1998). The simulation temperature was kept constant using a V-rescale thermostat (Bussi et al., 2007); protein and lipids, as well as the solvent, together with ions were coupled (τ = 0.1 ps) separately to a temperature bath of 310 K. Similarly, the pressure was kept constant at 1 bar by the Parrinello−Rahman barostat algorithm (coupling constant = 2 ps) (Bussi et al., 2007). Steepest descent energy minimization was performed prior to all simulations. To equilibrate the system, all heavy atoms of the protein were position restrained with a force constant (fc) of 1,000 kJmol⁻¹ nm⁻² for 2 ns NPT simulation, and then another 1 ns NPT simulation was followed without the position restraints. An electric field (40 mV·nm⁻¹) was applied along the channel pore from the intracellular side to the external side. With a box length of ∼14.6 nm, this amounts to a transmembrane potential of ∼580 mV. The Y146 side chain atoms were restrained throughout the simulations with a fc of 1,000 kJmol⁻¹ nm⁻² to prevent the SF from collapsing. This practice was necessary for simulating K⁺ ion conduction through Kir2.2 channels, as shown in our previous study (Zangerl-Plessl et al., 2020) and others (Jogini et al., 2023).

In the condition without PIP₂, analysis of 15 replica simulations revealed that removing PIP₂ molecules led to the loss of the α-helical conformation in their TM2 helices (Fig. S1), while the normal conformation was essentially maintained in simulations with PIP₂ (Fig. S2). In particular, loss of α-helical conformation in the middle of TM2 caused substantial structural deformation (Fig. S3) and kept the HBC gate wide and the channel highly conductive (Fig. S4). Our simulation setup caused these artificial states upon removal of PIP₂, which do not seem to be real and are discussed later. Therefore, α-helical content was used as a criterion to discern trajectories that are to be considered realistic. Only four simulations (Traj_1, 3, 11, and 13) showed channel closure with <10% loss of α-helical content in any residues of the TM2 helix. These four simulations were used to elucidate the conformational dynamics of Kir2 channel closure triggered by the removal of PIP₂ molecules.

To confirm that the protein did not undergo artificial deformation during simulations, we assessed the secondary structure of each residue of the TM2 helices throughout the simulations (every 200 ps), and the secondary structure was assessed using the VMD plugin tool, Timeline (Humphrey et al., 1996). An in-house MATLAB code was used to compute the fractional population of the α-helical content for each residue of each TM2 helix (Figs. S2 and S1).

PIP₂ interactions with neighboring residues were sampled every 1 ns, and the distances between phosphorus atoms and the terminal atoms of basic residue side chains were computed to determine distance changes (ΔDist) from the reference crystal structure (PDB accession no. 5KUM). These changes were averaged over the four subunits, concatenated from five trajectories, and presented as a violin plot (Fig. 2 b).

Diagonal Cα distance deviation of each residue (ΔDist) during MD simulations from the reference distance of Kir2.2 crystal structure (PDB accession no. 5KUM) was computed over each simulation. ΔDist of the HBC gate-forming residues was assessed for the two diagonal pairs (AC and BD) (Fig. 3 a; and Figs. S5 and S6). To examine potential PIP₂ effects on the protein structure, differential diagonal Cα distance (ΔΔDist) was computed by subtracting the averaged distance of each residue in the without PIP₂ simulation from those in the with PIP₂ simulation (Fig. 4 a). While RMSD measures the deviation of individual atom from a reference structure in Cartesian coordinates, ΔΔDist measures the changes in the diagonal distance of every residue independent of reference structures. These metrics are potentially more informative than RMSD for the following reasons: First, the diagonal distance changes of the gate-forming residues can be directly interpretable for the functional states of the channel. Second, it may be more intuitive to infer gross conformational changes from diagonal distance changes of spatially clustered residues. Third, the analysis does not require structural alignment, which prevents arbitrary results due to erroneous alignments. Conformational transitions were inferred after mapping ΔΔDist of every residue onto the Kir2 tetramer structures (Fig. 4, b and c). This strategy is commonly used in FRET (Wang et al., 2012) and EPR (Hubbell et al., 2000) approaches to predict conformational changes from the signal changes of the probes attached at multiple sites of a protein. The calculation was performed using in-house MATLAB codes.

Positional fluctuations of residues and the displacement correlations among the residues, referred to as dynamic cross correlation (DCC) (Kasahara et al., 2014), can reveal how different parts of the protein are structurally coupled to each other (McCammon, 1984). In this study, such analysis was performed to reveal how other parts of the protein are coupled to the fluctuation of the HBC gate-forming residues (I177 and M181). For the same advantageous reasons mentioned above, instead of instantaneous displacement of each Cα atom from its mean Cartesian coordinates, instantaneous diagonal distance displacement from the mean for each residue was employed (⁠$ΔCαj_AC = Cαi_AC – <Cα_AC>$⁠, where j is residue index and AC is the diagonal pair, and the angle bracket (< >) indicates time average).

The calculation for DCC coefficients (r_j) was performed using the built-in function, xcorr, in MATLAB for each trajectory and averaged over the last 400 ns of the simulations for with PIP₂ and without PIP₂ simulations, respectively (Fig. 5 b and Fig. S7 a).

Two diagonal pairs (subunit A/C and subunit B/D pairs) generated two sets of Cα distance fluctuation data (ΔCα_{j_AC} or ΔCα_{j_BD}) per trajectory. Fig. 5 a depicts how DCC analyses were performed. DCC coefficients were computed between each residue j (ΔCα_{j_AC}) and a reference residue (ΔCα_{177_AC}or ΔCα_{181_AC}). The calculation was repeated for every residue and performed in two different manners, either with the two entities from the same diagonal pair (“Intra” DCC) or with one from the AC pair and the other from the BD pair (“Inter” DCC). In doing so, how HBC gate motions are coupled with the rest of the protein within the same subunit pair or across the pair can be examined.

Two metrics were employed for quantitative comparison of the degree of fourfold symmetry in simulations with PIP₂ and without PIP₂. These metrics were computed using in-house MATLAB scripts after excluding flexible parts: the N-terminal loop (residues 41–60), the extracellular turrets and SF (residues 109–155), the extended loops (residues 238–242 and 334–337), and the C-terminal helix (residues 351–376).

A perfect fourfold symmetry will ensure that the two diagonal Cα distances of a residue are identical. As the protein drifts away from fourfold symmetry, the difference between the two diagonal distances will increase. The difference was averaged for each residue, and the averaged differences were summed over all residues, excluding the flexible parts as noted above. The mean and deviation of the sums were then compared between simulations with PIP₂ (n = 5) and without PIP₂ (n = 4) (Fig. 5 d).

Another geometric measure to enumerate the fourfold symmetry is the deviation of the angle (θ) formed by the two diagonal vectors from a right angle. In perfect fourfold symmetry, the two diagonal vectors bisect each other at 90°, but this angle will deviate as the fourfold symmetry is disrupted. For every residue, the deviation from 90° was computed and averaged. The averaged deviations were then summed over all residues, excluding the flexible parts as noted above. The mean and standard deviation of these sums were then again compared between simulations with PIP₂ and without PIP₂ (Fig. 5 e).

The motions were determined by measuring the angle between the two vectors as shown in Fig. S8 a. The TMD vector is a unit vector connecting the center of geometry (COG) of the rigid parts of the TMD (residues 79–100, 130–148, and 156–185) and the COG of the SF (residues 143–148). The CTD vector is a unit vector connecting the COG of the G-loop (residues 307–308) and the COG of the whole CTD (residue 187–350). The angle was computed for every snapshot, and the histograms computed for each trajectory were averaged. For simulations without PIP₂, only the last 400 ns were analyzed (Fig. S8 b).

Low-frequency protein motions were computed using the DynOmics webserver (Li et al., 2017), which allows assessment of the lowest frequency motions in the presence of environmental factors (Lezon and Bahar, 2012) by integrating two widely used ENMs—the Gaussian network model (GNM) (Yang et al., 2005) and the anisotropic network model (Eyal et al., 2006). The input files for this server were downloaded from the Orientations of Proteins in Membrane database (https://opm.phar.umich.edu/) (Lomize et al., 2012), in which all PDB files of membrane proteins contain dummy atoms designating the membrane location and enclosing the TMD. The first four lowest modes with the membrane implemented as an environmental factor were analyzed and compared. The conformational changes along each of the three unique modes (the second and the third lowest modes were redundant) were visualized individually (Videos 1, 2, 3, and 4). 29 snapshots showing the conformational changes along a mode with the maximum RMSD of 4 Å from the crystal structure were automatically generated by DynOmics server for each lowest mode. The trajectory was captured as a series of ppm files by VMD for each unique mode, which were converted into JPEG images, concatenated, and exported into an MP4 file by QuickTime Player 7 at 12 frames per second. The HBC and G-loop gate dimensions were determined for the two extreme structures by measuring the diagonal Cα distance of the gate-forming residues (I177 and M181 for the HBC gate, and A307 and M308 for the G-loop gate) (Table 2).

Statistical significance was analyzed using unpaired parametric t tests unless otherwise stated, using PRISM 10. Statistical significance of P < 0.05, P < 0.01, and P < 0.001 is indicated by single, double, and triple asterisks, respectively.

Fig. S1 shows the α-helical contents of TM2 helices in simulations without PIP2. Fig. S2 shows the α-helical contents of TM2 helices in simulations with PIP2. Fig. S3 shows the exemplary ribbon diagrams showing deformation of TM2 from α-helix in simulations without PIP2. Fig. S4 shows the time traces of K ion permeation in simulations without PIP2. Fig. S5 shows the time traces of the HBC gate dimension in simulations with PIP2 simulations. Fig. S6 shows the time traces of the HBC gate dimension in simulations without PIP2. Fig. S7 shows the DCC of all residues to the HBC gate residue (M181). Fig. S8 shows the clamshell motion in simulations with or without PIP2. Fig. S9 shows the time traces of wetting at the HBC gate in simulations with PIP2. Fig. S10 shows the time traces of wetting at the HBC gate in simulations without PIP2. Fig. S11 shows the time traces of K ion permeation in simulations with PIP2. Fig. S12 shows the time traces of individual residues, showing continuous deviation from the fourfold symmetry. Fig. S13 shows the properties of GNM mode spectrum. Fig. S14 shows the CTD twisting motions observed in Kir6.2 open channel structures. Fig. S15 shows the structural details of Kir2 channels behind the SF. Video 1 shows the first normal mode viewed from the extracellular side with every subunit colored uniquely. Video 2 shows the second and third normal modes viewed from the side with every subunit colored uniquely. Video 3 shows the fourth normal mode viewed from the bottom with every subunit colored uniquely. Video 4 shows the fourth normal mode viewed from the side with every subunit colored uniquely. Video 5 shows the cartoon diagram of Kir2.2 color coded by ∆∆Dist of every residue. Video 6 shows the cartoon diagram of Kir2.2 color coded by r_Intra for subunit A and C or by r_Inter for subunit B and D for simulations without PIP2. Video 7 shows the cartoon diagram of Kir2.2 color coded by r_Intra for subunit A and C or by r_Inter for subunit B and D for simulations with PIP2.

Full atomistic MD simulations were performed to probe the specific conformational changes that underlie channel gating within the PIP₂-bound compact conformation. We previously generated open and conducting Kir2.2 channels by MD simulation of the crystal structure of a gain-of-function mutant, K62W/G178D (KW/GD), with 4 PIP₂ molecules bound (Zangerl-Plessl et al., 2020) (see Fig. 1 b for the locations of the mutants and bound PIP₂). In silico, we reverted the introduced Asp (G178D) to the native Gly (KW/GD(G)) to obtain open, conducting conformations of the wild-type pore. This reversion resulted in the HBC gate being narrower and frequently dewetted, with reduced overall single channel conductance compared with KW/GD, paralleling experimental findings, but the channels stayed open (Fig. 3, left) (Zangerl-Plessl et al., 2020). With these open and conductive systems, we have now approached the question of whether we can obtain Kir2 channel closure following the removal of all PIP₂ molecules and, if so, identify PIP₂-dependent conformational dynamics that are directly coupled to Kir2 channel gating.

During the simulations with PIP₂ bound (“wPIP₂”), channels stayed open and conductive throughout the 1-μs simulations, evidenced by the mean Cα diagonal distances of the HBC gate-forming residues (I177 and M181) staying greater than those in the crystal structure (PDB accession no. 5KUM), i.e., ∆Dist > 0 (Fig. 3 a, left) and the HBC gate region staying predominantly wetted (albeit with intermittent brief dewetting, Fig. 3 b, left), and K⁺ ions continuing to permeate throughout the simulation (Fig. 3 c, left). In marked contrast, KW/GD(G) channels gradually closed at the HBC following removal of the four bound PIP₂ molecules (“woPIP₂”): ∆Dist at both I177 and M181 decreased from ∼2 Å to <1 Å (Fig. 3 a, right) on average and to less than that in the crystal structure (∆Dist ≤ 0) (Fig. 3 a, right). Moreover, the HBC gate became dewetted at the early stage of the simulation and stayed dewetted (Fig. 3 b, right), and K⁺ ion permeation stopped (Fig. 3 c, right). Parallel observations were made for multiple replica simulations regarding the mean Cα diagonal distances (∆Dist) of the HBC residues (Figs. S5 and S6), wetting and dewetting at the HBC gate (Figs. S9 and S10), and potassium ion permeation (Figs. S11 and S4). These findings suggest that PIP₂ is functionally effective at keeping the channel open in these simulations and that channel closure indeed occurs upon the loss of PIP₂ molecules. These sets of simulations, with or without PIP₂, could provide ample structural information to elucidate PIP₂-dependent conformational changes in the CTD that control Kir2 channel gating.

All PIP₂-bound Kir2 crystal structures to date have been sterically closed to permeation at the HBC (Lee et al., 2016), and the full conformational changes that connect PIP₂ interaction and channel gating remain unclear. In the available crystal structures, the three phosphate groups of PIP₂ make extensive salt bridge interactions with six basic residues within each subunit (Fig. 2 a, left). We measured and plotted the distance changes (ΔDist) between each phosphate group and the directly interacting residues in the PDB 5KUM crystal structure (Lee et al., 2016) during MD simulations (Fig. 2 b, left). While most of the salt bridge interactions remain similar to those in the crystal structure, a drastic increase in distance was observed for K189 in relation to both phosphate groups, P4 and P5. This indicates that K189 loses interaction with PIP₂ when the channel is open. Conversely, the distance between K188 and P4 significantly decreased, indicating closer apposition between the two moieties in the open Kir2 channels. One conformation sampled from a PIP₂-bound conductive simulation clearly demonstrates the structural basis of these distance changes (Fig. 2 a, right). In this conformation, K188 is pulled closer to PIP₂, while K189 swings away.

An unforeseen PIP₂ interaction with the K220 residue of the neighboring subunit is also apparent in the open Kir2 channels (Fig. 2 b, right). In the crystal structure, the K220 residue is located far from any PIP₂ molecules with >10 Å distances between the charged moieties (Fig. 2 c), indicating very weak electrostatic interactions between them. However, during MD simulations in the open state with PIP₂, the distances decreased by 4–5 Å on average (Fig. 2 b, right), with the actual distances around 8 Å, still indicating weak electrostatic interactions. Occasionally, the distance reductions reached even 8–10 Å (Fig. 2 b, right), making the actual distances only 3 and 4 Å, indicating strong salt bridge formation between them. Such a drastic distance narrowing cannot be made solely by side chain reorientation, and a representative structure observed in the MD simulations (Fig. 2 c) shows that the whole βC–βD loop, in which K220 resides (Fig. 1 c), is pulled closer to the PIP₂ molecule through rigid body motion of the whole CTD subunit.

To assess global conformational differences associated with PIP₂-dependent gating, we determined the differential diagonal distance (∆∆Dist) for all residues (Fig. 4 a) and mapped these onto the Kir2 ribbon diagrams (Fig. 4, b and c; and Video 5). Although the closed-state simulations originated from those with PIP₂ after its removal, the comparison between the two sets of simulations will proceed from without PIP₂ (closed) to with PIP₂ (open) for the remainder of the report to better demonstrate PIP₂-dependent motions leading to channel opening.

In the TMD, a dominant pattern observed indicated obvious widening near the bottom of the TMD with slight narrowing at the top of the TMD, while the SF was unchanged (Fig. 4 b). A global conformation change consistent with this pattern of diagonal distance changes could involve the TMD subunits being held tightly at the SF region while the bottom parts move away from each other potentially via clockwise rotational motions in each subunit (when viewed from the top of the channel), as previously suggested (Perozo et al., 1999; Shen et al., 2002; Shrivastava and Bahar, 2006). Such conformational changes by PIP₂ should lead to opening of the HBC gate at the bottom of the TMD.

In the CTD, a dominant pattern was observed, especially when viewed from the extracellular side (Fig. 4 c). The residues on the left side of the dotted line in each subunit got closer (blue colored), while those on the right side of the dotted line moved farther away (red colored). Global conformation changes to satisfy this pattern could involve each CTD subunit making intra-subunit twisting motions in the clockwise direction in the presence of PIP₂.

Of note is that, while PIP₂ led to an expansion of the critical HBC gate, the stabilization of the CTD twisting motions resulted in the G-loop gate being slightly more constricted in the simulations with PIP₂ (Fig. 4, a and b). This apparently opposite dependence of the two narrow regions on PIP₂ binding has not previously been considered but may be important in the coupling of Kir channel gating to distinct ligands, as discussed further below.

Covariance analysis of the positional fluctuation of atoms is a widely used tool to assess how motions in one part are correlated to other parts (McCammon, 1984; Kasahara et al., 2014). We assessed DCC) of all residues with the HBC gate-forming residues (I177 and M181) to determine the correlation of residues within the same diagonal dimer pairs (r_Intra) or in opposite pairs (r_Inter) relative to the gate-forming reference residues (Fig. 5 a, see Materials and methods). For the Intra pairs (Fig. 5 b, green), the residues (170–185) adjacent to and including the HBC gate showed the strongest positive DCC, with a maximum coefficient of 1 observed for the reference gating residue I177. Additionally, residues (60–100), which constitute the slide helix and the TM1 helix, also showed strong positive DCC to the fluctuation of the reference residue. Interestingly, the G-loop regions (306–310) were negatively correlated with HBC gate fluctuations (Fig. 5 b). The r_Intra profiles were quite similar between the two simulation conditions with or without PIP₂ (Fig. 5 b, green, top and bottom).

The DCC coefficients of the Inter pairs (Fig. 5 b, pink) were generally small and opposite in sign to those of the Intra pairs, indicating dynamic motions with a twofold anticorrelation between the two opposing pairs (AC versus BD pairs). This anticorrelation was pronounced in simulations without PIP₂; the extent of inversion of correlation coefficients (green and pink lines) was smaller in simulations with PIP₂ (Fig. 5 b). This was also evidenced by the slope of the linear regression of the scatter plots of r_Intra versus r_Inter shifting from negative to positive with PIP₂ (Fig. 5 c), indicating augmented positive correlation between the two opposing pairs by PIP₂. Similar results were obtained for both HBC gate-forming residues I177 (Fig. 5, b and c) and M181 (Fig. S7, a and b).

The reduced anticorrelations may imply a PIP₂-driven increase in fourfold symmetry. To test this, we employed two different metrics (see Materials and methods): the sum of the diagonal distance differences (Fig. 5 d) and the sum of angle deviations from 90° (Fig. 5 e) between the two diagonal pairs. The more symmetric the structure, the lower the values will be. Both measurements were higher in simulations without PIP₂, although statistical significance was achieved only for the distance deviations. To test whether these results originate from mere fluctuations, the difference values over the course of simulations were considered (Fig. S12). Either negative or positive values were observed for some residues for long period of time, indicating that deviation from fourfold symmetry was maintained for such extended times and that the observed augmented sums in simulations without PIP₂ were not due to increased fluctuations. Collectively, these results support the conclusion that PIP₂ decreases the anticorrelated motions while increasing concerted motions among the four subunits.

To interpret the correlation coefficients in terms of structural dynamics and to probe how PIP₂ increases concerted motions, the cross-correlation coefficients for the Intra and Inter pairs were mapped onto ribbon diagrams of A/C and B/D subunits, respectively (Fig. 6 a), for simulations without (Video 5) or with (Video 6) PIP₂.

In the absence of PIP₂, the bottom of the TMD of the Inter pairs (B/D), as well as the top of the CTD of the Intra pair (A/C, with subunit C omitted for clarity), are negatively correlated with the dynamics of the HBC gate of the Intra pair (A/C) (Fig. 6 a, left). This suggests that expansion of the HBC gate between A/C subunits leads to the collapse of the nearby B/D subunits (Fig. 6 a, left, arrows), as well as the top of the CTD of A/C subunit pairs (Fig. 6 a, left, filled arrowhead). Conversely, the top of the CTD of the Inter (B/D) pairs was positively correlated with the dynamics of the HBC gate of the Intra pair (A/C) (Fig. 6 a, left, open arrowhead). These observations indicate that, while the dynamic coupling between the TMD and CTD within a subunit is weak, the top of one CTD subunit is dynamically coupled to the TMD of its neighboring subunit. This staggered coupling between the TMDs and CTDs in the absence of PIP₂ is illustrated in the schematic diagram (Fig. 6 b, left).

PIP₂ clearly changed how different parts of the protein were dynamically coupled. PIP₂ attenuated the extent of negative correlations at the bottom of the two TM helices of the Inter pairs (B/D) (Fig. 6 a, right, arrows) and at the top of the CTD of the Intra pairs (A/C) (Fig. 6 a, right, filled arrowheads), while decreasing the positive coupling of the CTD apex of the Inter pair (B/D) with the dynamics of the HBC gate of the intra pair (A/C) (Fig. 6 a, right, open arrowhead). This suggests that PIP₂ increases positive coupling between the four subunits as well as between the HBC and the G-loop gates within each subunit, as illustrated in the schematic diagram (Fig. 6 b, right).

DCC coefficients in the CTD also showed a vertically inverted pattern between the two pairs of subunits, which was more pronounced in the absence of PIP₂ (Fig. 6 b; and Videos 6 and 7). This suggests a scenario in which the top of each of the A and C subunits moves inward while the bottom moves outward, and concurrently the top of the B and D subunits moves outward while the bottom moves inward. Such anticorrelated tilting motions could directly expand the G-loop region, which might then be positively coupled to the HBC gate motions in the presence—but not the absence—of PIP₂ (Fig. 8 d), and thus could be critical for HBC gate opening.

Coarse-grained NMA (Brooks and Karplus, 1983; Go et al., 1983) or ENM analysis (Tirion, 1996; Bahar et al., 1997; Frappier and Najmanovich, 2014) have been widely used to elucidate functionally important large-scale conformational dynamics encoded in protein structures. To compare our observations from atomistic simulations, ENM analysis was performed for the PIP₂-bound Kir2.2 PDB 5KUM crystal structure, using the web server, DynOmics (Li et al., 2017), and the first four lowest frequency modes with high collectivities (Fig. S13, a and b) were examined. The conformational changes along each mode are shown in Fig. 7 and in Videos 1, 2, 3, and 4. The first (referred to as mode I) involves rotation of the CTD relative to the TMD, i.e., a global back-and-forth torsional motion around an interfacial plane between the TMD and the CTD (Fig. 7, top, and Video 1). The second and third modes (collectively referred to as mode II) are degenerate, involving clamshell motions between the TMD and CTD (Fig. 7, middle, and Video 2). The fourth mode (referred to as mode III) shows anticorrelated twofold motions within the lower part of the CTDs of the two pairs of diagonally apposed subunits. When one pair (A and C subunits) moves toward each other, the other pair (B and D subunits) moves away from each other (Fig. 7, bottom, and Videos 3 and 4).

The global conformational changes observed in the MD simulations are qualitatively similar to those in modes I and III in the ENM analysis (Fig. 7). Mode I may reflect clockwise twisting of individual CTD subunits with a pivot in each subunit in the opposite directions from anticlockwise rotational TMD motions (Fig. 4 c). Similarly, the motions in mode III may reflect the twofold tilting motions in the CTD seen in the MD simulation (Fig. 6 b). What then of the observed motions in the identified mode II, i.e., “clamshell” bending, where TMD and CTD move toward each other? We assessed the angle between the molecular axes of the TMD and the CTD in the MD simulations (Fig. S8 a, see Materials and methods). There were bending motions, but the mean bending angle distributions were unaffected by PIP₂ (and 3.81 ± 0.89° and 3.73 ± 1.38° in the absence and presence of bound PIP₂, respectively) (Fig. S8 b), suggesting that they may not be functionally important in the gating motions of Kir2 channels.

Kir channels are gated by multiple distinct ligands that bind to the CTD and control the opening and closing of the critical HBC gate, allowing pore-wetting and ion permeation. PIP₂ binding to a well-defined site located at the CTD and TMD interface is typically required for opening of all Kir channels, but conformational changes coupling PIP₂ binding and the gating transitions themselves remain unclear. Here, we analyzed Kir2.2 structural dynamics over the course of channel closure by running MD simulations that started with a Kir2 channel in the PIP₂-bound, open, and conductive state and underwent closure following removal of all four PIP₂ molecules. This approach provided conformational ensembles in both closed and open states, allowing for the elucidation of PIP₂-dependent conformational dynamics relevant to channel opening.

Our study reveals conformational differences between PIP₂-bound, open, Kir2 channels and the crystal structures of PIP₂-bound, but closed, channels (Hansen et al., 2011; Lee et al., 2016). Of the six basic residues that contact PIP₂ in the crystal structures, major deviations were observed with K188 and K189. Additionally, K220, from the neighboring subunit, which is positioned far from PIP₂ in the crystal structures, becomes much closer to PIP₂ in the open channels (Fig. 2). The increased distance between K189 and PIP₂ suggests that the interaction of K189 with PIP₂ may not be essential for maintaining the channel in an open state. This observation aligns with findings from studies on Kir3.2, where a mutation of the equivalent Lys residue to Tyr results in channel activity that is comparable with or even higher than the wild-type channel (Lacin et al., 2017). Aromatic amino acids, such as Tyr, favor membrane-water interfaces and can help to stably anchor sites, such as the C-linker of Kir channels, to the membrane surface (Hong et al., 2007). The K189 residue may be important for normally pulling the CTD toward the membrane, and once the CTD adopts the compact “up-conformation,” this Lys may interact with other anionic lipids to keep the up conformation.

In contrast, the neighboring K188 residue interacted more strongly with PIP₂ in the open channel than in crystal structures, suggesting this interaction may be critical for mediating the PIP₂-driven changes leading to the open channel. Indeed, mutation of this Lys residue resulted in loss-of-function in both Kir2.1 (D’Avanzo et al., 2013) and Kir3.2 (Lacin et al., 2017) channels.

The closeness of the K220 side chain to PIP₂ in the open channel was unexpected, but the drastic movement required could be explained by the movement of the whole CD loop on which K220 resides. Our previous study showed that the equivalent Lys residue mutation to Cys in Kir2.1 (K219C) resulted in a partial loss-of-function that could not be rescued by decyl modification (Lee et al., 2013), which pulls the protein toward the membrane and generates the up conformation. This implicates K220 in structural changes beyond generation of the up conformation and suggests it may be involved in the stabilization of PIP₂ binding in the open state.

Torsional motions in opposing directions, between the top and bottom halves, are frequently identified in proteins with cylindrical shapes (Taly et al., 2005; Shrivastava and Bahar, 2006; Chennubhotla et al., 2008; Bahar et al., 2010a). Multiple crystal structures and single-particle EM structures of Kir channels also imply marked rotations of the CTD relative to the TMD around the central pore axis (Hite et al., 2015; Martin et al., 2017; Niu et al., 2020; Zhao and MacKinnon, 2021). However, correlations between such whole-domain rotations and channel gating have remained obscure, with some studies proposing no correlation of CTD rotational motions with gating events (Bernsteiner et al., 2019; Fagnen et al., 2020). Our atomistic MD simulations indicate that rotation also occurs within each CTD in the presence PIP₂. Twisting of each CTD in the clockwise direction (viewed from the extracellular side, with a pivot residing near the C-linker of each subunit) is favored when PIP₂ is bound and is positively correlated with a wider HBC gate. PIP₂-dependent twisting within individual CTDs was first identified in the Kir3.2 (GIRK2) R201A mutant crystal structure (PDB accession no. 3SYQ) but in an anticlockwise direction (Whorton and MacKinnon, 2011), opposite to what is observed in our simulations. This structure was determined in complex with only two PIP₂ molecules, and the C-linker and slide helix interactions were missing, potentially due to these anticlockwise twisting motions, which would drive the C-linker moving away from the slide helix of the neighboring subunit. The R201A mutant channel also exhibited very low activity compared with the wild-type protein (Whorton and MacKinnon, 2011), questioning the relevance of such anticlockwise twisting motions for wild-type protein-gating transitions.

Recent wild-type cryo-EM structures of the Kir3.2 channel, each with four PIP₂ bound, showed no such anticlockwise twisting motions and maintained interaction between neighboring subunits through the C-linker and the slide helix (Niu et al., 2020; Mathiharan et al., 2021). On the other hand, recent studies on Kir6.2 channels reveal open-state structures, induced either by gain-of-function mutations (Zhao and MacKinnon, 2021) or by PIP₂ binding (Driggers et al., 2024), revealed significant clockwise twisting motions that were potentiated by PIP₂ binding (Fig. S14). In contrast to eukaryotic Kir channels, bacterial KirBac1.1 channels are closed by the binding of PIP₂ (D’Avanzo et al., 2010). We previously employed intramolecular FRET to experimentally identify conformational changes associated with PIP₂-driven KirBac1.1 channel closure (Wang et al., 2012). PIP₂ resulted in differential FRET efficiency changes among the CTD residues in a pattern that also implies clockwise twisting motions in each CTD. This is consistent with the PIP₂-dependent twisting motion observed in our current study, suggesting that PIP₂ binding drives similar conformational changes in both prokaryotic and eukaryotic Kir channel CTDs, even though the final functional consequences of the binding are opposite (D’Avanzo et al., 2010).

Notably, while we observed that PIP₂ potentiated the clockwise twisting motions in the CTD that drive opening of the HBC gate, they constricted the G-loop region (Fig. 8 c). Similar effects were observed in previous MD simulations of Kir3.2 channels, where PIP₂ widened the HBC gate but narrowed the pore in the G-loop region, unless two other cofactors, Na⁺ and Gβγ were present (Bernsteiner et al., 2019; Li et al., 2019). We recently showed that, while PIP₂ activates normally reduced Kir3.2 channels, it reverses gating in oxidized Kir3.2 channels, which now become spontaneously highly active but are inhibited by PIP₂ (Lee et al., 2023). This reversal of the PIP₂ effect could be explained by G-loop constriction in response to PIP₂ binding, leading to channel closure in this oxidized state (Lee et al., 2023).

Taken together, both experimental and modeling studies suggest that PIP₂ drives similar twisting motions in all Kir channels, but that the coupling of these major conformational changes to functional states of the channels may differ based on distinct local structures and the detailed coupling between the CTD and the TMD, analogous to the opposite voltage-sensing domain coupling to pore gating between HCN channels and other voltage-dependent channels (Männikkö et al., 2002; Vemana et al., 2004).

Analysis of diagonal distance fluctuations revealed consistently anticorrelated dynamics between the two diagonal subunit pairs in each channel. These anticorrelations were attenuated by PIP₂, particularly near the HBC and G-loop regions. This attenuation would tend to stabilize open channel states by promoting movement of the four subunits in the same radial direction (Fig. 8 e). Within each CTD, correlation coefficients were largely inverted vertically, from being positive to negative in one pair or negative to positive in the other pair (Fig. 6 b). This pattern can be explained by tilting motions in which the upper parts of one diagonal pair move inward, and the lower parts of the same diagonal pair move outward, concurrent with opposite motions in the other diagonal pair. Such tilting motions in the CTD were also suggested by the intramolecular FRET assays of prokaryotic KirBac1.1 proteins (Wang et al., 2012), again implying similar PIP₂-driven conformational changes in both eukaryotic and prokaryotic Kir channels. Such twofold anticorrelated motions may promote channel closure due to constriction along one diagonal axis, even as the other expands (Fig. 8 e).

Coarse-grained NMA has been widely used to elucidate structure-encoded low-frequency/highly cooperative motions of proteins that are often related to function (Brooks and Karplus, 1983; Brooks, 1995). Unlike previous studies done with the isolated TMD (Shrivastava and Bahar, 2006), CTD (Haider et al., 2005), or whole Kir2 structures with the CTD disengaged (Fernandes et al., 2022), we here performed the analysis for a Kir2 structure with the CTD in the compact up conformation. Dimension changes in both the HBC and G-loop region were negligible (Table 2) throughout the motions along the first three unique lowest frequency modes, even though modes I and III matched to the two lowest modes from previous studies that widened the gates (Haider et al., 2005; Shrivastava and Bahar, 2006). This result implies that tight engagement with another domain restricts large deformation in general and may allow only small and hence reversible conformational changes underlying the actual gating transitions. Indeed, a very restricted opening of HBC gates due to tight engagement of regulatory domains has been identified in other ion channels; the opening at the HBC gate was 10 Å narrower for the full-length KcsA protein than for a CTD-truncated mutant (Uysal et al., 2011), and a very wide opening of MthK channel was only enabled by release of the TMD from the regulatory RCK domain upon Ca²⁺ binding (Fan et al., 2020). Clear cyclic nucleotide–triggered reorganizations around the cyclic nucleotide–binding region were seen in isolated CNG channel CTDs, but these were also absent following cyclic nucleotide binding in the whole protein (Lee and MacKinnon, 2017; James and Zagotta, 2018). Also, a recent apo Kir3.2 EM structure (PDB accession no. 6XIS) showed striking widening of the G-loop region with the CTD fully disengaged from the TMD (Niu et al., 2020).

Conformational changes consistent with the lowest frequency motions were found to be important for gate opening in our all-atom MD simulations; mode I was largely consistent with the observed twisting motions, and mode III with the observed intra-domain tilting motions. This suggests that the low-frequency motions implied by coarse-grained analysis are indeed valid and may be useful in future to guide examination of large-scale motions in all-atom simulations, in which too many details can screen the cardinal information.

Cooperative subunit motions driving channel gating is a common feature of many ion channels, potentially important in allowing exquisitely sensitive responses to changes in the levels of different regulatory stimuli (Tytgat and Hess, 1992; Niu and Magleby, 2002; Qian et al., 2006; Xie et al., 2008; Benndorf et al., 2012; Guan et al., 2022). PIP₂-dependent gating of Kir2 channels is also cooperative (Xie et al., 2008), and here we identify the primary effect of PIP₂ binding as being to potentiate intra-subunit twisting about a pivot residing in each CTD subunit, with the four subunits making the movement simultaneously (Fig. 8, b and c). PIP₂ binding facilitates channel opening by decreasing the anticorrelated motions between pairs of diagonally opposed subunits (Fig. 8, d and e), thereby increasing the channel fourfold symmetry, and by enhancing the physical coupling between the TMD and the CTD within each subunit.

In our simulations, we observed two unexpected and irreversible conformational deformations. First, the SF became deformed if Y146 side chain atoms were not restrained, and second, in some simulations, the TM2 helices lost α-helical conformation during channel closure after PIP₂ removal.

SF deformation was due to rotation of the Y146 residue, which resulted in the loss of hydrogen bonding to T140 of the neighboring subunit (Fig. S15 a). This was observed in our previous simulations (Zangerl-Plessl et al., 2020) and has also been observed by another group using different starting structures, force fields, and simulation packages (Jogini et al., 2023). In contrast, simulated Kir3 channels (Bernsteiner et al., 2019; Chen et al., 2020) maintain the conductive SF conformation without restraints, suggesting that the SF deformation may be a unique feature of Kir2 channels. Comparison of the Kir3 and Kir2 structures reveal a void space behind the Kir2, but not Kir3 SF (Fig. S15 b), with smaller hydrophobic residues in the vicinity (Fig. S15c). This SF deformation mimics C-type inactivation, which, however, has not been observed in Kir2 channels when examined by electrophysiological approaches.

Y146 rotation might occur in Kir2 channels but be reversible and short-lived. While simulated proteins can become trapped in local minima during the achievable timescales (i.e., a few μs) of conventional all-atom MD simulations, reversibility on a slightly longer timescale would be hidden in the open-state noise or could be undetectably fast in single channel experiments with 1–10 kHz resolution. Consistent with this suggestion, recent time-resolved x-ray crystallography showed transient deformations of the K channel SF by strong electric fields on the timescale of a few hundreds of nanoseconds (Lee et al., 2025).

M2 helix deformation, as we observed in our simulations, was not reported in a previous study examining PIP₂ removal and subsequent channel closure (Jogini et al., 2023). We suggest the possibility that the protonation state of the rectifying controller residue (D173) (Fig. 1 b) in the inner cavity might explain this discrepancy. In the simulations by Jogini et al. (2023), two of the Asp residues were set to be negatively charged, while the other two were protonated and neutral. In contrast, we set all four Asp residues to be negatively charged based on reasons described in our previous report (Zangerl-Plessl et al., 2020). While all four Asp residues may remain negatively charged in the open conformation, protonation might occur as the channel closes and the distances between these Asp residues decrease (Cymes and Grosman, 2011; Maksaev et al., 2023). Strong repulsive forces between four negatively charged Asp residues could have caused α helical deformation as the channel closed following PIP₂ removal. Occasionally, charge screening by mobile K⁺ ions among the four Asp residues in the inner cavity might have led to 4 (out of 15) without PIP₂ simulations that maintained the M2 helical state throughout. Carrying out simulations at constant pH could be an intriguing approach to determine whether the rectifying controllers of Kir2 channels indeed switch charge states during channel activation cycles.

The initial and end-stage PDB files and the preprocessed .TPR file of each simulation are openly available in a public repository, Zenodo (https://doi.org/10.5281/zenodo.17092388) for download, and the individual MD trajectories and in-house MATLAB codes will be available upon requests.

Crina M. Nimigean served as editor.

This work was supported by CIMED Pilot and Feasibility grant CIMED-21-02 and the National Institutes of Health (NIH) R03 grant TR003670 (to Sun-Joo Lee), NIH grant R35 grant HL171542 (to Colin G. Nichols), and Austrian Science Fund grants W1232 (to Anna Stary-Weinzinger) and ZK-81B (to Eva-Maria Zangerl-Plessl). Simulations were carried out using the Vienna Scientific Cluster and XSEDE (TG-MCB190046).

Author contributions: Eva-Maria Zangerl-Plessl: conceptualization, data curation, formal analysis, funding acquisition, investigation, software, visualization, and writing—review and editing. Anna Stary-Weinzinger: resources, supervision, validation, and writing—review and editing. Colin G. Nichols: conceptualization, funding acquisition, supervision, and writing—review and editing. Sun-Joo Lee: conceptualization, data curation, formal analysis, funding acquisition, investigation, methodology, project administration, resources, software, supervision, validation, visualization, and writing—original draft, review, and editing.