Molecular Dynamics Simulations of KirBac1.1 Mutants Reveal Global Gating Changes of Kir Channels

Prokaryotic inwardly rectifying (KirBac) potassium channels are homologous to mammalian Kir channels. Their activity is controlled by dynamical conformational changes that regulate ion flow through a central pore. Understanding the dynamical rearrangements of Kir channels during gating requires high-resolution structure information from channels crystallized in different conformations and insight into the transition steps, which are difficult to access experimentally. In this study, we use MD simulations on wild type KirBac1.1 and an activatory mutant to investigate activation gating of KirBac channels. Full atomistic MD simulations revealed that introducing glutamate in position 143 causes significant widening at the helix bundle crossing gate, enabling water flux into the cavity. Further, global rearrangements including a twisting motion as well as local rearrangements at the subunit interface in the cytoplasmic domain were observed. These structural rearrangements are similar to recently reported KirBac3.1 crystal structures in closed and open conformation, suggesting that our simulations capture major conformational changes during KirBac1.1 opening. In addition, an important role of protein–lipid interactions during gating was observed. Slide-helix and C-linker interactions with lipids were strengthened during activation gating.


■ INTRODUCTION
Inwardly rectifying potassium (Kir) channels are intrinsic membrane proteins that control the selective permeation of potassium ions across otherwise ion impermeable cell membranes. The primary role of Kir channels is the regulation of outward directed K + current. Under physiological conditions, Kir channels generate a large inward K + conductance at potentials negative to the equilibrium potential of potassium (E K ) but permit less outward current flow at potentials positive to E K . These channels are regulated by many different cellular factors such as ATP, intracellular pH, phosphatidylinositol-4,5bisphospate (PIP 2 ), and nonspecific secondary anionic phospholipids. 1−3 In particular, lipid modulators such as PIP 2 , as well as cholesterol, have been shown to regulate bacterial and eukaryotic Kir channels. Remarkably, while cholesterol has been suggested to inhibit both pro-and eukaryotic Kir channels, the effect of PIP 2 is opposite. Namely, while PIP 2 is essential for activation of eukaryotic Kir channels, bacterial channels such as Kirbac1.1 have been shown to be inhibited by this phospholipid. 4−7 Over the last 10 years, several Kir crystal structures of the cytoplasmic domain as well as several full length structures of prokaryotic and eukaryotic channels have been published. 8−18 Interestingly, prokaryotic homologues, share similar architecture and have functional activity with eukaryotic channels, 19 despite only moderate sequence conservation on the amino acid level. All Kir channels undergo dynamical changes to regulate ion flow. This process, referred to as "gating", involves structural rearrangements of the transmembrane (TM) as well as the cytoplasmic domains (CTD). In the closed conformation, ion flux is prevented by a narrowing of the inner TM2 helices, which form a constriction site at the helix bundle crossing (HBC) gate close to the intracellular side (see Figure  1). Computational studies on KirBac1.1 channel models provided insights into ion selectivity and gating dynamics. 20−24 A limitation of all these studies was the lack of open state X-ray structures. In 2012, the first X-ray structure of a bacterial Kir channel in a presumably open conformation was crystallized, using a known activatory "gain-of-function" mutant. 16 Comparison of this open structure with various Kir channels in closed conformation 13 provides insights into gating induced changes of these channels. In the open structure, global conformational changes are observed, including a rotational movement of the CTD relative to the plane of the membrane; in addition, a bending of the TM2 at a highly conserved glycine opens the HBC gate. The sequential process of gating events remains a major open question. Especially, the cross-talk between TM and CT domains and how this leads to channel opening is still unknown. To investigate these events, we performed MD simulations on the KirBac1.1 channel, for which only closed state X-ray structures are available. 13 ■ METHODS Molecular Dynamics Simulations. The closed (pdb identifier: 2WLL) crystal structure, 13 comprising residues 38 to 308, was used as starting point for MD simulations. The G143E mutant in protonated and deprotonated conformations and the R153A mutant were generated with the software Swiss-PdbViewer. 25 pK a values for all titratable amino acid side chains were calculated with PROPKA. 69 The structures were embedded in an equilibrated membrane consisting of 256 palmitoyloleoylphosphatidylcholine (POPC) lipids using the g_membed tool, 26 which is part of the gromacs package. K + ions were placed in the SF at K + sites S0, S2, and S4 with waters placed at S1 and S3. 27 Cl − ions were added randomly within the solvent to neutralize the system. All simulations were carried out using the gromacs simulation software v.4.5.4. 28 The amber99sb force field 29 and the TIP3P 30 water model were employed for the protein and water, respectively. Lipid parameter for the POPC membrane were taken from Berger et al. 31,32 The corrected monovalent ion Lennard-Jones parameters for the amber force field were used. 33 Electrostatic interactions were calculated at a distance smaller than 1.0 nm, long-range electrostatic interactions were treated by the particle-mesh Ewald method at every step. 34 Lennard-Jones interactions were calculated with a cut off of 1.0 nm. The LINCS algorithm 35 was used to constrain bonds. Modeling hydrogens as virtual sites 36 allowed for an integration step of 4 fs. 28 The Nose-Hoover thermostat 37,38 was used to keep simulation temperature constant by coupling (tau = 0.2 ps) the protein, lipids and solvent (water and ions) separately to a temperature bath of 310 K. Likewise, the pressure was kept constant at 1 bar by using the Parrinello−Rahman barostat algorithm 39 with a coupling constant of 1 ps. Prior to simulation, 1000 conjugate gradient energy-minimization steps were performed, followed by 5 ns of equilibrium simulation in which the protein atoms were restrained by a force constant of 1000 kJ mol −1 nm −2 to their initial position. Lipids, ions, and water were allowed to move freely during equilibration. Four times 200 ns MD simulations were performed for the full length WT channel as well as the G143E d , G143E p , and G143E d -R153A mutant channels.
Salt Bridge Analysis. Electrostatic interactions were analyzed by measuring the center of mass distances between positively and negatively charged functional groups of amino acids. A distance cut off of 6 Å was set which represents three different types of ion pair interactions, namely salt bridge, N−O bridge, and long-range ion pair 40 which are all referred to as "salt bridges" in this study. The occurrence of interaction is normalized to the most prominent electrostatic interaction in the protein (R193 and E187 of adjacent SU in WT simulations). Interaction partners that contributed more than 1% to the total electrostatic interactions in the protein are plotted in the star graphs.
Calculation of the Rotational Angle of the CTD. Measuring the changes in the angle between TM and CTD was done by calculating the torsion angle between two planes. The first plane was assigned to three points: the center of mass of the backbone atoms of resides 65−68 (sequence ASLA) in the TM region of one subunit, the center of mass of the same backbone atoms residues of all subunits and the center of mass of the backbone atoms of residues 225−227 (sequence: GWN) of the cytoplasmic region of all subunits. The second plane was assigned as follows: the center of mass of the backbone atoms of residues 225−227 (sequence: GWN) of one subunit and the same two points that include all subunits as described before. At the beginning of the simulation, the angle between these two planes was defined as zero in order to calculate the changes during the simulation.
Energy Profile Calculations. Potential of mean force (PMF) calculations were performed as described previously. 41 Briefly, the main conformational changes in the most prominent G143E d opening simulation were obtained by principal component analysis. The first eigenvector obtained by this PCA was used as a reaction coordinate. Along this reaction coordinate, 45 windows were chosen for umbrella sampling and simulated for 50 ns. Umbrella sampling simulations were performed by applying a harmonic restraint force along the transition pathway with force constants between 1 and 100 kJ mol −1 nm −2 . The first 30 ns of each window were discarded for equilibration. The potential of mean force and the statistical errors of the activation gating energy profile were estimated by making use of the g_wham tool of gromacs and the integrated bootstrap analysis method. 42 The number of bootstraps was set to 100.

Article
Mutations and Growth Assay. KirBac1.1 WT coding DNA was inserted between NcoI and HindIII of pQE60 vector. 43 All mutations were introduced into KirBac1.1 by sitedirect mutagenesis kit (Agilent Inc.) and confirmed by DNA sequencing.
For growth assay, 20 ng of KirBac1.1 plasmids were transformed into E. coli BL21-gold (DE3) host strain following the protocol provided by the manufacture and 10 μL of transformants were spotted on LB agar plates containing 100 μg L −1 of ampicillin and 15 μg L −1 of tetracycline. The plates were incubated at 37°C overnight, and then, pictures were taken by digital camera.

■ RESULTS
To probe the mechanism of KirBac1.1 gating, we made use of the known activatory ("gain-of function") mutant G143E. 44 This mutant was selected due to its equivalent position to activatory mutant S129R in KirBac3.1, which was used to obtain open state crystals of this channel. 16 G143E is located in transmembrane helix 2 (TM2) at a hydrophobic interface between two adjacent TM2 helices (see Figure 1A and B and supplemental Figure 1). The activatory effect of this mutant was investigated using four times 200 ns unbiased full atomistic MD simulations of the full length KirBac1.1 WT crystal structure and mutant G143E in deprotonated (denoted as G143E d ) and protonated (G143E p ) form.
Mutant G143E d Induces Opening of the HBC Gate. MD simulations show that mutant G143E d induces global conformational rearrangements of the protein. Bending at a highly conserved glycine hinge (G134) in TM2, leading to opening at the HBC, was observed in all four simulations. A HOLE plot showing the pore diameter after 200 ns is shown in supplemental Figure 2A and B. To monitor the changes at the gate over time, we measured the Cα−Cα distance between opposing F146 residues, lining the narrowest point of this gate. As shown in Figure 1C, the distance rapidly increased to 13.8 ± 0.9 Å, compared to WT simulations, where the gated remained fully closed (Cα−Cα distance at F146: 11.0 ± 0.9 Å). The end state of the G143E d mutant was compared with the KirBac3.1 open X-ray structure. Figure 1D shows a structural superposition of the TM2 helices of the two structures, revealing good overlay between the structures. Next, the χ1 angle distribution of the F146 side chain over time was analyzed. As shown in Figure 2A−C, the χ1 angle switched from predominantly ∼160°in WT (cavity facing) to predominantly ∼270°(cavity lining) in the G143E d mutant channel.
To further investigate the consequence of these structural changes on the HBC gate, we monitored the water flux through the gate in WT and G143E d simulations. While water flux was not observed in the WT simulations (see Figure 3A and C), considerable water migration through the gate occurred in the G143E d mutant ( Figure 3B and C).
Global Conformational Changes in the Cytoplasmic Domain. In addition to the rearrangements at the HBC, our simulations revealed large conformational changes at the CTD. A rotational movement of the CTD relative to the plane of the membrane was seen in all four G143E d mutant simulations. The degree of this twisting motion amounted to 15°on

Journal of Chemical Information and Modeling
Article average, with maximum values of 23°in one run (see Figure  4A). These values are in good agreement with data obtained from several KirBac3.1 X-ray structures, 13,16 suggesting that the rotational movements of these two channels are conserved.
Moreover, rearrangements at the subunit interface, especially salt bridge formations, were analyzed. In this study, the term "salt bridge" denotes nonbonded, N−O bridged, and longrange electrostatic interactions between acidic carboxyl groups and basic amino groups in the same subunit (sSU) or adjacent SUs (aSUs) as described by Kumar et al. 40 In the WT closed structure, R271, located in the β I strand (see Figure 4B, D), forms a salt bridge with E262 (G-loop of the adjacent subunit). Further, E187 interacts with K191 and R193 from the β D strand of the aSU. Moreover, hydrogen bonds between R193 and E218 from the adjacent CD-loop were observed.
In the G143E d mutant, global conformational rearrangements of the CTD led to an additional salt bridge between R271 and E187 of the aSU ( Figure 4B, E). This salt bridge formation occurred within the first 80 ns in all four simulations between all four interfaces. Due to the R271−E187 salt bridge formation, interactions between K191 and E187 were weakened.
Salt Bridge Formation between G143E d and R153. Further, the structural changes at the TM−CTD interface were examined. Our analysis revealed that the G143E d side chains form a stable salt bridge (see Figure 5A and C) with residue R153, located in the C-linker of the neighboring subunit, within the first half of the simulations. To investigate the importance of this salt bridge for the cross-talk between TM and CTD, R153A was introduced in the background of the G143E d mutant in all four subunits. In these simulations, the HBC gate opens on average to 13 Å (not shown), but the observed rotation of the CTD was rather small with ∼5°( Figure 6). This suggests that the strong electrostatic interactions between G143E d and R153 are important for the twisting motion observed in the G143E d mutant.
Influence of Protonation State on Channel Conformation. The above-described observations indicate that opening involves a two-step process. First, strong repulsion between G143E d and the surrounding hydrophobic residues triggers opening at the bundle crossing region. In a second step, electrostatic interactions between G143E d and R153 of the adjacent C-linker induce rotation of the CTD. It was previously reported that the activity of KirBac1.1, 43 as well as of a close homologue 18 are both pH-dependent, thus we investigated the influence of the protonation state of G143E on channel gating. In repeated simulations with G143E p (protonated), neither opening at the HBC nor twisting at the CTD were observed (see Figure 7A, B, supplemental Figure 2C). Moreover, no water flux was observed within 200 ns (not shown).
Energetics of the G143E d Mutant Channel Opening. To investigate the coupling between the HBC gate and the CTD twisting motion in more detail, we calculated the free energy landscape of activation gating (Figure 8). The main conformational changes in the most prominent G143E d opening simulation (CTD rotation of 23°) were obtained by principal component analysis and used as reaction coordinate for umbrella sampling calculations. At the beginning of the simulation a steep energy decrease of ∼7 kcal/mol was observed. During this phase, the HBC gate opened and the rotameric state of the F146 side chain changed from a cavity facing to a cavity lining conformation. Further, a first rotational movement of the CTD of ∼12°occurred. In addition, monitoring of the G143E d −R153 salt bridge along the reaction coordinate revealed that in 3 of the 4 subunits a salt bridge between TM2 and the linker of the adjacent subunit was formed during this phase ( Figure 8D). From 6 to 7 nm, a plateau phase ( Figure 8A) was observed, where no rotational movement of the CTD occurred ( Figure 8C). Subsequently, a second rotational movement of the CTD led to a total rotation of 23°compared to the starting structure and a further decrease in energy of ∼3 kcal/mol.
Protein−Lipid Interactions during Gating. The importance of KirBac1.1 phospholipid interactions at the TM−

Journal of Chemical Information and Modeling
Article CTD interface was reported previously. 13,45 Thus, we investigated the protein−lipid contacts in this region. Figure  9 shows the number of hydrogen bonds to the lipid head groups over time. While for R49 no gating dependent effect was seen, the number of hydrogen bonds increased for K57 (slidehelix) and R151 (C-linker) during channel opening. Downward movement of the slide-helix at the C-terminal end and a subtle outward movement repositions K57, leading to increased lipid exposure of this residue. Further, repositioning of the C-linker induced by the G143E d −R153 salt bridge, led to an upward movement of the R151 side chain and strengthened lipid contacts as shown in Figure 9C.
Experimental Testing of the G143E Mutant. We have constructed the KirBac1.1 G143E mutant and attempted to express and purify the protein for functional assay. Unfortunately, the protein appears to be very toxic to E. coli host strain, since we are unable to obtain any transformant that expresses mutant protein at detectable levels. Previously, our studies

Journal of Chemical Information and Modeling
Article showed that KirBac1.1 WT plasmid expresses active potassium channels, but G143C mutant loses channel function. 46 In the present work, we transformed pQE60 vector, KirBac1.1 WT, G143E, and G143C (20 ng of each), into equal amount (28 μL) of E. coli BL21-gold (DE3) competent cells. As shown in Figure 10, plenty of transformants were obtained for vector and G143C mutant. For KirBac1.1 WT, the transformation efficiency was decreased, but we obtained only few colonies for G143E. Our data are consistent with the G143E plasmid generating highly active channels, which are very toxic and kill the host stain.
Comparison to Experimental Data for KirBac Gating Motions. The transition pathways obtained by our MD simulations are in good agreement with experimental data. Recent FRET experiments on KirBac1.1 channels revealed major molecular motions in the CTD induced by PIP 2 binding. 47 Remarkably, during our 200 ns simulations with the G143E d mutant, we observed tilting motions of the β I sheet, as reported from these FRET experiments. Additionally, X-ray structure analyses on the homologue KirBac3.1 channel revealed a twisting motion of the CTD of 23°relative to the plane of the membrane. 13 The final KirBac1.1 state obtained by simulating the G143E d mutant is in excellent agreement with the twisted conformation described for the homologous KirBac3.1 channel. 13,16 The structural changes observed for the HBC gate (bending at a glycine hinge) are consistent with data on other K + channels. However, the extent of channel opening varies among crystallized structures. 48−52

■ DISCUSSION
In this study, we investigated conformational changes of KirBac1.1 gating by taking advantage of the activatory mutant G143E d in TM2. Simulations with this mutant revealed detailed mechanistic insights into the gating of Kir channels.
In all G143E d mutant simulations, HBC opening occurred prior to conformational changes at the CTD. The introduction of a negatively charged glutamic acid in a hydrophobic pocket between the TM2 helices led to strong repulsion, which enabled opening of the HBC (Figure 1 and 2). This finding is supported by a previous MD simulation on a KirBac6.1 homology model. 44 Further, this region has previously been shown to have dramatic effects on gating. For example, a KirBac3.1 open state X-ray structure was crystallized by mutating the equivalent position S129 to an arginine. 16 Interestingly, an activatory mutant (A108T/S) was also identified at this site in the bacterial K + channel KcsA 53,54 as well as in eukaryotic Kir channels. 55 Moreover, we recently showed that the equivalent position is conserved in voltage gated calcium channels 56 and mutation of this position has substantial effects on channel gating. We and others have indicated that the small size at this position seems critical for stabilizing the closed gate. 53,57−59 Comparing G143E simulations in protonated and deprotonated state (Figure 1 and 7), revealed that the effect of the mutant in KirBac1.1 on gating is primarily resulting from the negative charge of the side chain and to a lesser extent a size effect. This suggests that several factors contribute to the effects of mutants close to the HBC gate.
The limited simulation time of 200 ns makes it difficult to assess, whether the HBC gate is fully or only partially opened in our simulations. However, water flux observed in all G143E d mutant simulations indicates an open state (Figure 3). Substantial conformational variation of the open HBC gate have been reported in X-ray structures 9,16,60,61 indicating that subtle variations between channel species might exist. For example, structural differences in CTDs might influence gating. Additionally, there is accumulating evidence that several open states exist for each channel, as reported for KirBac1.1. 19 In previous X-ray structures of the KirBac3.1 channel, twisted and nontwisted CTD conformations were obtained only with a closed HBC gate. 13 This led to the conclusion that the CTD rearrangements trigger HBC opening. Contrary to this previously suggested gating model, 16,18 our simulations revealed

Journal of Chemical Information and Modeling
Article that the CTD conformational changes can occur after HBC gate opening (Figures 1 and 8). This indicates that the coupling between TM and CTD might operate bidirectionally. The mutant might influence the cross-talk between the two domains. Nevertheless, the importance of electrostatic interactions for stabilizing the twisted conformation, as predicted in our simulations ( Figure 5), is in excellent agreement with previous data on KirBac3.1. 13 In the G143E d mutant, these interactions are mainly accomplished by salt bridge formation of the mutant side chain with R153. The significance of this contact is further stressed by results from the R153A mutant simulations, where only subtle twisting motions were observed ( Figure 6). Although this interaction can only occur in the G143E d mutant, end states obtained from simulations closely resemble native twisting motions of the CTD as inferred from KirBac3.1 structures. Moreover, subunit interface rearrangements predicted by our simulations (Figure 4) are similar to KirBac3.1. 13 Another important prediction from our simulations concerns the pH-dependence of mutant G143E. Only the deprotonated glutamic acid induced global conformational changes on the nanosecond time scale, suggesting that gating of this mutant might be pH dependent (Figures 1 and 7). Similar observations were reported for a F168E mutant in the HBC gate of the mammalian Kir6.2 channel. 62 Taken together, our data provide structural details of how protonatable side chains can be used to induce channel gating by pH titration.
We attempted to validate our findings by expressing and purifying the protein in an E. coli host strain. Unfortunately, the G143E mutant appears to be very toxic and kills the host stain, which suggests that the mutant generates highly active channels in agreement with our simulations, revealing that global gating rearrangements in both the HBC gate and the CTD with the G143E d mutant are accessible via MD simulations. This indicates that the mutant might significantly decrease the energetic barrier for channel opening, since all WT X-ray structures to date were captured with a closed HBC gate. Indeed, our PMF calculations revealed an energy difference of ∼10 kcal/mol between closed and open state, with no energy barriers present (Figure 8). It is conceivable that an energetic barrier needs to be overcome in WT for channel opening as shown in a previous simulation study on KcsA. 41 There is accumulating evidence, highlighting the importance of lipid components for regulating Kir channels (for recent reviews, see refs 63−67). Analysis of protein lipid interactions in WT and G143E d mutant simulations revealed gatingdependent hydrogen bond formation. In particular, interactions of K57 located in the slide-helix and R151 from the C-linker to the lipid head groups are significantly increased upon channel opening ( Figure 9). These observations are in excellent agreement with a study by Enkvetchakul et al. 45 which reported the importance of lipid head groups in regulating KirBac1.1 gating. Interestingly, additional nonspecific anionic lipid interactions have been recently shown to be required for Kir2 channel gating. 3 This indicates that all Kir channels are strongly lipid regulated, further supported by a recent MD study on a Kir3.1 chimera. 68 In conclusion, the presented simulations unravel the progression of conformational changes during gate opening. Contrary to previous hypotheses based on static crystal structures, opening of the HBC gate can trigger twisting of the CTD. This process is mediated by electrostatic interactions between TM and CT domains. Additionally, lipid contacts with the slide-helix facilitate channel opening and presumably stabilize this conformation. One has to keep in mind however that our simulations are based on the G143E "activatory" mutant. It cannot be excluded that the gating transitions of wild type Kir channels differ. Even though, the gating transitions observed in our simulations are in good agreement with recent FRET experiments, 47 wild type open state X-ray structures of Kir channels in combination with MD simulations will be needed to validate the proposed gating mechanism.