INTRODUCTION

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION CHALLENGES REFERENCES

Potassium channels play a central role in the membrane physiology of both excitable and nonexcitable cells. They are a large family of integral membrane proteins (for example, there are thought to be ~60 K channel genes in Caenorhabditis elegans; Choe et al., 1999). The sequence diversity within the K channel family reflects its functional diversity, particularly in terms of gating mechanisms and kinetics. However, all K channels share a common functional feature in that they are selective for K⁺ over Na⁺ ions. This is mirrored by a common sequence motif, the presence of a TVGYG or related fingerprint sequence in the pore-lining P-region. The x-ray structure of a "minimalist" bacterial K channel (KcsA from Streptomyces lividans; Doyle et al., 1998) revealed how the TVGYG motif forms the selectivity filter of K channels. The pore-forming domain has the topology M1-P-M2, where M1 and M2 are transmembrane (TM) -helices and where the re-entrant P-hairpin is made up of an -helix (the P-helix) plus a TVGYG-containing loop. This loop adopts a noncanonical extended conformation made possible by the two glycines. In the channel tetramer the four TVGYG motifs are brought together to form a narrow, oxygen-atom lined pore, i.e., the selectivity filter. This sits inside a bundle of eight TM -helices (two helices, M1 and M2, from each subunit). In more complex K channels, such as the voltage-gated (Kv) channels, there are four additional TM helices per subunit, presumed to be packed around the central eight-helix bundle. In the two P-domain K channels, such as the background conductance channel TWIK, the functional channel is made up of two subunits, rather than the more usual four, with each subunit containing two copies of the M1-P-M2 fold.

The P-loops of the various K channels have strong sequence homologies, indicating that they share a similar architecture. This suggests that it may be possible to generate accurate homology models of mammalian K channels by using the bacterial KcsA structure as a template. Such models may provide insights into functional (e.g., permeation) properties of mammalian K channels. For example, Repunte et al. (1999) have used a model of Kir6.2 (see below) to help explain how extracellular residues (in the M1-P loop) contribute to the single-channel conductance difference between Kir6.1 and Kir6.2. To investigate this further we decided to look at a family of mammalian K channels that are thought to share the same topology as KcsA, namely the inward rectifier K channels (the Kirs; Doupnik et al., 1995; Reimann and Ashcroft, 1999). Kir channels have a number of physiological roles, including regulation of cardiac activity and control of insulin release by the  cells of the pancreas. They also provide a test of the validity of the assumption underlying homology modeling, i.e., that there exists a shared architecture between KcsA and other K channels, as this has recently been questioned by Minor et al. (1999) with respect to Kir2.1.

In this paper we are concerned with one Kir in particular, namely Kir6.2. Kir6.2 is the channel-forming subunit of the K_ATP channel of the  pancreatic cells. In vivo, K_ATP channels are formed by a complex of four Kir6.2 subunits plus four sulfonylurea receptor (SUR1) subunits. However, Tucker et al. (1997) have shown that a C-terminally truncated form of Kir6.2 is able to homotetramerize to form functional K channels in the absence of SUR1, at least in vitro. Thus it is reasonable to model a Kir6.2 tetramer, using KcsA as a template. Other Kir channels may form heterotetramers (Reimann and Ashcroft, 1999), which would make modeling a little more difficult.

There have been a number of site-directed mutagenesis (SDM) studies of Kir channels that aid the evaluation of homology models (reviewed by, e.g., Reimann and Ashcroft, 1999). We reserve detailed discussion of these until later. However, one key residue identified by SDM is N160 (see Fig. 1), which is in the central section of the M2 helix of Kir6.2. This residue is replaced by an aspartate (D) in Kir2.1. It has been shown that mutation of this D to N (or Q) in Kir2.1 changes channel permeation properties, particularly block by Mg²⁺ ions, such that they become closer to those of Kir1.1 and Kir6.2, both of which have an N at this position (Wible et al., 1994; Stanfield et al., 1994). Mutation of N to D in Kir1.1 confers Mg²⁺ block (Lu and MacKinnon, 1994). These results suggest that N160 should form part of the pore-lining region in Kir6.2. Of course, there does remain an element of uncertainty in exactly how to interpret this SDM result in structural terms. In particular, it seems likely that the KcsA x-ray structure may represent the closed form of the channel (Perozo et al., 1998, 1999), whereas the Mg²⁺ block mentioned above occurs when the Kir is in its open conformation.

View larger version (39K): [in this window] [in a new window]   FIGURE 1   Sequence alignment of Kir6.2 vs. KcsA. Only the sequences around the modeled TM domain are given. The numbers above refer to the Kir6.2 (human) sequence. The vertical arrows indicate the beginning and end of the Kir model (175) and of the KcsA x-ray structure (173). The gray boxes indicate the positions of the helices in KcsA. The black-outlined boxes around the Kir sequence indicate the positions of the TM helices identified in an earlier study (Shrivastava et al., 2000). The gray-outlined box indicates the selectivity filter sequences. The star indicates the N160 residue of Kir6.2 discussed in the text.

In this paper we present a homology model of the pore domain of Kir6.2. This model is used as the starting point for nanosecond-duration molecular dynamics (MD) simulations in a phospholipid bilayer. Preliminary analysis of the interactions of the Kir6.2 channel model with K⁺ ions and water molecules during these simulations suggests that single-file motion of K⁺ ions and water through the selectivity filter takes place, in a fashion similar to that observed in simulations of KcsA (Shrivastava and Sansom, 2000).

	METHODS

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION CHALLENGES REFERENCES

Sequence alignment

Sequence alignments were performed using Amps (Barton and Sternberg, 1990), aided by various other web-based alignment tools. Alscript (Barton, 1993) was used to display final alignment. Prediction of the location of the TM helices in Kir6.2 employed a number of methods, namely, TopPred2 (http://www.biokemi.su.se/~server/toppred2; von Heijne, 1992); TMAP (http://www.embl-heidelberg.de/tmap/tmap  info.html; Persson and Argos, 1994; Persson and Argos, 1997); DAS (http://www.biokemi.su.se/~server/DAS; Cserzo et al., 1997); PHDhtm (http://www.embl-heidelberg.de/predictprotein; Rost et al., 1995, 1996); memsat (http://globin.bio.warwick.ac.uk/~jones/memsat.html; Jones et al., 1994); TMHMM (http://www.cbs.dtu.dk/services/TMHMM-1.0/; Sonnhammer et al., 1998); and TMPred (http://www.isrec.isb-sib.ch/software/TMPRED_form.html; Hofmann and Stoffel, 1993). These sequence-based methods were supplemented by MD simulations of isolated M1 and M2 helices embedded in a lipid bilayer. In this manner we optimized our predictions of the extent of the TM helices of Kir6.2 (Shrivastava et al., 2000). Analysis of sequence periodicities for the predicted TM helices used Perscan (Donnelly et al., 1993).

Homology modeling

Homology models were generated using Modeller v4 (Sali and Blundell, 1993). Fourfold rotational symmetry was imposed. An ensemble of 25 model structures was generated. These were ranked by analysis of their stereochemistry, using Procheck (Morris et al., 1992), and by their agreement with the KcsA coordinates. From the ensemble, a single structure was selected for further analysis and as a starting structure for the MD simulations.

pK_A calculations

The pK_As of ionizable residues within the Kir6.2 homology model were estimated as described previously (Adcock et al., 1998; Ranatunga et al., 1998). For these calculations we used UHBD, version 5.1 (Davis et al., 1991) (with some local modifications), and partial atomic charges from the Quanta/Charmm22 parameter set.

Molecular dynamics simulations

The Kir6.2 model was embedded in a palmitoyloleoylphosphatidylcholine bilayer and used as the starting structure for nanosecond MD simulations, employing the same protocols used previously for KcsA (Shrivastava and Sansom, 2000). MD simulations were run using GROMACS v1.6 (http://rugmd0.chem.rug.nl/~gmx/gmx.html). A twin-range cutoff was used for long-range interactions: 1.0 nm for van der Waals interactions and 1.7 nm for electrostatic interactions. The time step was 2 fs. The LINCS algorithm was used to constrain bond lengths. Normal pressure and temperature were used in the simulation. A constant pressure of 1 bar independently in all three directions was used, with a coupling constant of _p = 1.0 ps (Berendsen et al., 1984). Water, lipid, and protein were coupled separately to a temperature bath at 300 K, using a coupling constant _T = 0.1 ps. MD simulations were performed on an SGI Origin 2000, typically taking ~21 days of cpu time on a single 195-MHz R10000 processor.

The lipid parameters were based on those used in MD studies of dipalmitoylphosphatidylcholine bilayers (Berger et al., 1997; Marrink et al., 1998), with the addition of some GROMOS parameters for the double bond in the acyl tail. The lipid-protein interactions used GROMOS parameters. The water model used was SPC (Berendsen et al., 1981), which has been shown to be a reasonable choice for lipid bilayer simulations (Tieleman and Berendsen, 1996). The K⁺ parameters (kindly supplied by Dr. Peter Tieleman) were as in Straatsma and Berendsen (1988). These parameters have been used in MD simulations of KcsA (Shrivastava and Sansom, 2000) and of the channel-forming peptide alamethicin (Tieleman et al., 1999a,b).

Other computational details

Structural diagrams were prepared using Molscript (Kraulis, 1991) and Raster3D (Merritt and Bacon, 1997). The projections of the inner lining of the pore were generated using PORELINING (Smith and Sansom, unpublished; also Adcock et al., 1998; Mongan et al., 1998). Pore radius profiles were calculated using HOLE (Smart et al., 1993). Secondary structure analysis used DSSP (Kabsch and Sander, 1983). Other analyses used GROMACS and/or locally written code.

	RESULTS

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION CHALLENGES REFERENCES

Sequence alignment

An optimal sequence alignment is essential to the success of homology modeling. The sequence identity between KcsA and Kir6.2 is ~30% in the P-hairpin region, making this element of the alignment relatively straightforward. However, for the M1 and M2 helices the sequence identity is rather low. To optimize the alignment of the M1 and M2 regions, the initial alignment was adjusted so as to maximize overlap between the predicted locations of the TM helices in Kir6.2 and their locations in the x-ray structure of KcsA. This provided a first-pass refinement of the alignment of the helices. We then used periodicity analysis to further refine the alignment. For M1 a clear -helical periodicity of substitution patterns for surface and buried residues of membrane proteins (Donnelly et al., 1993) was observed when the sequences of the seven classes of Kir channels were aligned. The M1 alignment was adjusted so that the residues of the Kirs predicted to be lipid exposed were aligned with the exposed side chains on the external surface of KcsA. However, such periodicity analysis yielded ambiguous results for M2, presumably reflecting the more complex environment of M2 in the channel protein. Instead, the initial alignment of M2 was adjusted to optimize the overlap of predicted (Kir) and observed (KcsA) M2 helices, and such that N160 of Kir6.2 corresponded to a pore-lining location in KcsA (see above).

The alignment thus obtained is shown in Fig. 1. Note that the extent of the alignment used in subsequent homology modeling did not extend beyond the corresponding termini in KcsA x-ray structure. The first 57 residues and the last 215 residues of Kir6.2 are missing from the homology model, which thus corresponds to the central pore-forming domain of the channel. The loop between the end of M1 and the start of the P-helix is longer in Kir6.2 than in KcsA. Within the P-helix the LWW motif of KcsA is replaced by FLF, and the TVGYG selectivity filter sequence of KcsA is replaced by TIGFG. The V-to-I substitution in the latter is conservative. However, the Y-to-F substitution in the filter may have interesting functional implications (see below). Note that the alignment of M2 conserves G165 of Kir6.2. It is possible that this conserved glycine may confer some degree of flexibility upon M2, which in turn may be related to channel gating (Shrivastava et al., 2000).

Comparison of folds

As expected, the KcsA and Kir6.2 folds are similar (C RMSD = 0.17 nm for the 90% of the C atoms falling within an RMSD upper cutoff = 0.35 nm; Fig. 2). The P-helix and selectivity filter exhibit very similar (main chain) conformations. The most striking difference between the two structures is the longer M1/P loop in Kir6.2 (19 residues in Kir6.2 versus nine residues in KcsA). This loop has probably adopted a somewhat arbitrary conformation in the homology model. This aspect of Kir homology models in general will need further attention, especially if these studies are to be extended to other Kirs and their interactions with channel-blocking toxins (which generally bind in the vestibule to which the M1/P loops contribute), such as tertiapin (Jin and Lu, 1998; Jin et al., 1999) and d-dendrotoxin (Imredy et al., 1998). The other noticeable difference between KcsA and the Kir6.2 model is the rather longer pre-M1 "tail" in the latter.

View larger version (71K): [in this window] [in a new window]   FIGURE 2   Schematic diagrams of the structures of KcsA (A and C) and Kir (B and D). In A and B the view is down the pore axis, from the extracellular end of the protein. In C and D the view is perpendicular to the pore axis, with the extracellular mouth at the top of the diagram and only two subunits shown (for clarity). In D the M1 and M2 transmembrane helices, the P helix, the selectivity filter loop (F), and the N and C termini of the peptide chain are indicated. This and subsequent structure diagrams were drawn with Molscript (Kraulis, 1991) and Raster3D (Merritt and Bacon, 1997).

Comparison of outer surfaces

As both KcsA and Kir6.2 are integral membrane proteins, then, at least in vitro, the outer surfaces of both structures are exposed to a lipid bilayer environment. As this is an anisotropic environment one would expect a corresponding anisotropy in the outer surfaces of the two channel proteins, as is the case in other integral membrane proteins. Such considerations were not (explicitly) included in the alignment and homology modeling procedure. So, visualization of the outer surface of the Kir6.2 model and comparison with that of KcsA provides a simple, if somewhat crude, evaluation of the plausibility of the homology model. Examination of the distribution of polar versus apolar side chains on the surfaces of the two structures (Fig. 3) reveals that both KcsA and Kir6.2 exhibit a segregation of side chains common to most integral membrane proteins. Thus the polar residues are clustered at either end of the molecule (corresponding to the lipid headgroup/water interfacial region of the bilayer), and in each case there is a broad central band of predominantly hydrophobic side chains.

View larger version (100K): [in this window] [in a new window]   FIGURE 3   Molecular surfaces of KcsA (A and C) and Kir (B and D) compared. In each case the molecule is viewed perpendicular to the pore axis, with the extracellular mouth at the top. In A and B the hydrophobic residues are in pale gray and the polar residues are in dark gray. In C and D all of the residues are in pale gray, except for Trp and Tyr, which are in dark gray. Both polar and Trp/Tyr side chains can be seen to form bands on the surface of the molecules that correspond to the presumed location of the lipid headgroups in a bilayer.

We also examined the distribution of amphipathic aromatic residues (i.e., Trp and Tyr). Several studies (Schiffer et al., 1992; Yau et al., 1998) have indicated that these two residues are preferentially located in two rings on the surface of membrane proteins, corresponding to the position of the lipid headgroup/water interface when the protein is embedded in a bilayer. This location may reflect the ability of such residues to form favorable H-bonding interactions with lipid headgroups (as suggested by some simulation studies, e.g., Forrest et al., 1999) or may be a more general consequence of the shape and aromaticity of such side chains (Yau et al., 1998). Comparison of the distribution of Trp and Tyr residues on the outer surfaces of KcsA and Kir6.2 reveals rings of such side chains at either end of the molecule. The approximate distances between the two rings are 3.0 nm for KcsA and 3.4 nm for Kir6.2. These values are consistent with a bilayer thickness of ~3 nm. Thus the outer surface of the Kir6.2 molecule looks like that of an integral membrane protein. Of course, this is not strong proof that the homology model is correct, but it does suggest that the global distribution of side-chain types is compatible with what we expect for a membrane protein.

Comparison of pore-lining surfaces

We have also examined the inner, pore-lining surfaces of KcsA and Kir6.2 (Fig. 4). Both structures have five rings of O atoms that form the selectivity filter, namely the O of T, and the carbonyl oxygens of the T(V/I)G(Y/F) backbone. These rings of oxygens occupy very similar locations in the two structures, indicating that the homology model has preserved the essentials of the selectivity filter structure. Turning to the central cavity, in KcsA this is lined by predominantly apolar side chains, with the exception of T107 (which may form a transient interaction site for K⁺ ions inside the cavity; Shrivastava and Sansom, 2000). Interestingly, N160, which has been implicated in interactions with ions within the Kir6.2 channel (see above), occupies a position within the Kir6.2 pore close to (i.e., one turn of the M2 helix away from) that of T107 in KcsA. Thus in both the bacterial channel and the homology model, the properties of the lining of the cavity appear to be conserved.

View larger version (30K): [in this window] [in a new window]   FIGURE 4   Pore-lining diagrams of KcsA (A) and Kir (B). In these diagrams the four subunits are spread-eagled onto a plane, with their pore-lining surfaces uppermost. The extracellular mouth corresponds to the center of the diagram. The selectivity filter oxygens are shown in dark gray, key cavity-lining residues (T107 in KcsA and N160 in Kir) are shown in mid-gray, and all other residues in pale gray.

In particular, the equivalents of the following residues in Kir 6.2 have been identified as pore-lining by a combination of Cys substitution and Cd²⁺ block (Lu et al., 1999; Loussouarn et al., 1999): N153, I154, L157, A161, I162, L164, C166, F168, M169, T171, A172, and Q173. All of these are part of the pore lining in our Kir6.2 model. Furthermore, Yang et al. (1997) have suggested that a salt bridge is required for selectivity filter stabilization in Kir2.1, the equivalent residues in Kir6.2 being E126 and R136. In the homology model the mean distance between these pairs of residues is 0.65 (± 0.07) nm. This is rather long for a salt bridge interaction, and so it will be important to monitor the interactions of these residues during the MD simulations (see below).

One may also examine the nature of the side chains at either mouth of the pore. In KcsA there is a ring of acidic side chains at either mouth of the pore. The extracellular mouth of KcsA is guarded by a ring of D80 side chains; the intracellular mouth is guarded by a ring of E118 side chains. In Kir6.2 the general feature of acidic residues at either mouth is conserved, although the detailed disposition of these side chains differs from that in the bacterial channel. Thus the extracellular mouth has rings of E104, E108, and E141 side chains, and the intracellular mouth has a ring of D58 side chains.

Pore dimensions

It is informative to compare the dimensions of the pore in KcsA and the Kir6.2 model. This may be done by calculating the pore radius as a function of position along the channel axis (Fig. 5). Such a diagram makes clear the three regions of the channels: 1) the intracellular "gate" (ICG), 2) the central cavity (C), and 3) the extracellular selectivity filter (F). In both regions ICG and F the radius of the pore falls below the ionic radius of a K⁺ ion (0.13 nm). The radius profiles in the region of F are very similar for KcsA and Kir6.2, again confirming that the filter structure is maintained in the homology model. In contrast, the ICG is somewhat narrower for Kir6.2 than for KcsA (minimum radii ~0.05 versus 0.1 nm, respectively). This appears to be due to the tight packing of the F167 side chains that form the ICG in Kir6.2, in contrast to the V115 side chains that occupy the same position in KcsA.

View larger version (18K): [in this window] [in a new window]   FIGURE 5   (A) Pore radius profiles, calculated using HOLE (Smart et al., 1993), of KcsA () and Kir (·····). The intracellular mouth of the pore is at z  30 Å and the extracellular mouth at z  +15 Å. The approximate locations of the intracellular gate (IGC), cavity (C), and filter (F) regions are indicated. (B) Cumulative number of water molecules (see text) from the intracellular mouth to the extracellular mouth. Note the step in each curve corresponding to the central cavity (z  5 Å).

A further difference between KcsA and Kir6.2 is the size of the central cavity. This is somewhat smaller in Kir6.2. This difference has a consequence in terms of the number of water molecules that may be expected to occupy the cavity, which may influence the energetics of ion permeation. A lower bound on this number may be obtained by integrating the pore radius profile (to obtain a cumulative volume as a function of distance moved along the pore axis) and dividing this volume by that of a water molecule (i.e., 0.03 nm³). If we perform this calculation for KcsA, there is a step in the number of water molecules at the cavity corresponding to ~12 waters (note that 14-18 waters were fitted in the cavity by the solvation procedure used in the MD simulations; Shrivastava and Sansom, 2000). In contrast, for Kir6.2 the step in cavity volume corresponds to ~8-10 water molecules. The solvation procedure added 12 waters. Thus the Kir6.2 cavity is somewhat smaller than that of KcsA. (Note that this gives a lower bound on the number of waters because of the way the pore radius profile was evaluated (Smart et al., 1993, 1997), by fitting a largest possible sphere at each point along the axis and then reporting the volume of that sphere. Use of an alternative procedure, via calculation of a solvent-accessible surface, yields an upper bound of ~23 waters for the cavity of KcsA and ~18 for the cavity of Kir6.2 (Smart and Sansom, unpublished results).) Of course, it must be remembered that these structures may resemble the closed conformations of the two channels more closely than the open conformations (see below), and it is conceivable that in the open conformations more waters may be accommodated.

Ionizable side chains

The pK_A calculations indicated that certain side chains in the Kir6.2 homology model had their pK_As significantly perturbed from their standard values. Thus, at neutral pH, only two of the four D58 residues, three of the E126, three of the R136, and two of the K170 residues were predicted to be ionized. Furthermore, all four E140 and all four E141 were predicted to be protonated. Note that the protonation of two of the D58 residues reduces the magnitude of the ring of negative charge at the intracellular mouth of the channel. This has some similarities to the results of similar calculations for KcsA (Ranatunga and Sansom, unpublished results), which indicate suppression of ionization of the E118 residues at the intracellular mouth. The suppression of ionization of one E126 and one R136 (in the same subunit) will prevent formation of a salt bridge (see below).

Kir6.2 in a lipid bilayer

MD simulations of the Kir6.2 model embedded in a lipid bilayer (Fig. 6 A) were performed for two reasons: 1) to "refine" the homology model via simulations in an explicit bilayer plus solvent environment, by aiding identification of the less conformationally stable and hence (presumably) less well-modeled regions of the structure; and 2) to explore the nature of Kir6.2's interactions with K⁺ ions and water molecules, allowing comparison with the results of a previous simulation of KcsA (Shrivastava and Sansom, 2000). Accordingly, the Kir6.2 model was inserted into a palmitoyloleoylphosphatidylcholine lipid bilayer that contained a hole of dimensions sufficient to accommodate KcsA. Some minor adjustment of lipid molecules was needed to accommodate Kir6.2 by removing steric conflicts. The resultant system was solvated, and counterions were added by replacing water molecules at the positions of lowest Coulomb potential, to give an overall electroneutral system. Note that the solvation procedure resulted in 12 waters being added to the central cavity. To optimize the conformations of the lipids after protein insertion, a 50-ps simulation was run, during which the protein atoms were restrained to their initial positions while the lipids and water were free to move, and a lateral (xy) pressure of 500 bar was applied. This was followed by an equilibration run of 150 ps, during which the protein restraints were retained, but the pressure was returned to 1 bar along all three coordinate axes.

View larger version (62K): [in this window] [in a new window]   FIGURE 6   Simulation system MDK2 at the start of the MD run. (A) Kir (dark gray ribbons), lipid (pale gray bonds and dark gray spheres for P atoms of headgroups), and five of the 12 K+ ions (pale gray spheres). Water molecules are omitted (for clarity). (B) The selectivity filter showing two of the four subunits, two K+ ions (K1 and K2), and the intervening "x-ray" water molecule (Wx). The carbonyl oxygens of the residues that make up the selectivity filter and the four binding sites (S1 to S4) are indicated by arrows. S1 is formed by two rings of carbonyl oxygens (from F133 and G132); S2 is formed by two rings of carbonyl oxygens (from G132 and I131); S3 is formed by two rings of carbonyl oxygens (from I131 and T130); and S4 is formed by a ring of carbonyl oxygens (from the T130 main chain) and a ring of hydroxyl groups (from the T130 side chains).

Two initial simulations were performed (MDK0a and MDK0b; see Table 1). K⁺ ions were not added to the protein in these simulations. Instead, only Na⁺ counterions were present, located within the "bulk" water regions on either side of the bilayer. These two simulations differed in their treatment of the ionizable side chains within the protein. In MDK0a all residues were assumed to be in their default ionization state, i.e., acidic residues (Glu and Asp) were negatively charged and basic residues (Lys and Arg) were positively charged. Histidine residues were neutral. This gave a net charge of 20 on the Kir6.2 model. In simulation MDK0b, the residues were assigned ionization states based on their calculated pK_A values, assuming neutral pH (see above). This reduced the net charge on the protein to 12. Each of these simulations was run for 1 ns. The "drift" of the simulated structure from the initial model was evaluated by estimating the C RMSD as a function of time. The drift was significantly greater for MDK0a (final RMSD ~0.3 nm) than for MDK0b (final RMSD ~0.2 nm). Thus it seems that adjustment of the residue ionization states to match their calculated pK_A values results in a more "stable" simulation. On the basis of this preliminary comparison we decided to focus on the simulation with the pK_A-adjusted ionization states. However, examination of both the MDK0a and MDK0b models at the end of the simulation revealed that their selectivity filter regions had become distorted away from the conformation seen in KcsA. Previous simulation results for KcsA suggested that its selectivity filter became distorted over the duration of a prolonged MD simulation if no K⁺ ions were present in the filter region, but not if K⁺ ions were included in the filter. Therefore, we decided to run Kir6.2 simulations with pK_A-adjusted ionization states and with K⁺ ions included within the filter of the starting structure of the channel.

In the KcsA x-ray structure, K⁺ ions are present simultaneously at two sites of three within the selectivity filter, and a third K⁺ ion is suggested to be present in the central cavity. We therefore set up two simulations of Kir6.2 that included K⁺ ions. In MDK2, two K⁺ ions were incorporated, at sites S1 and S3 of the filter (Fig. 6 B), and with a water at site S2 between the two K⁺ ions (as is also seen in the x-ray structure of KcsA). The starting coordinates of the two ions and the water were taken as being similar to those in KcsA, i.e., in the centers of the rings of oxygen atoms that define the sites forming the selectivity filter. Thus the distance K1-O(Wx) = 0.26 nm and K1-O(Wx) = 0.44 nm. The setup for simulation MDK3 was similar, with the addition of a third K⁺ ion at the center of the cavity. Simulation MDK3 was run for 1.0 ns and simulation MDK2 for 1.6 ns. In both cases an initial 50-ps equilibration simulation was run, during which the coordinates of the protein atoms and the two (or three) K⁺ ions were restrained, while the lipid and solvent atoms were free to move. Monitoring of the dimensions of the simulation box as a function of time suggested that these parameters had fully relaxed during the equilibration period.

The progress of simulations MDK2 and MDK3 was also monitored in terms of the drift from the initial protein structure (Fig. 7). This revealed no significant difference between the two simulations. In both cases the C RMSD after 1 ns was ~0.25 nm. (For MDK2 no further increase in C RMSD took place between 1.0 and 1.6 ns.) This is similar to the degree of drift seen in simulations of the OmpF porin trimer (Tieleman and Berendsen, 1998) and of KcsA (Shrivastava and Sansom, 2000). It is encouraging that the drift of the Kir6.2 model is similar to that of two x-ray structures of membrane proteins (i.e., OmpF and KcsA). The initial rise in the C RMSD over the first 200-300 ps is common in such simulations and may be attributed to relaxation of the protein upon its transfer to a bilayer environment and/or inaccuracies in the potential function. The absence of a significant difference in overall RMSD between the two simulations suggests that the presence/absence of the third (i.e., cavity) K⁺ ion does not influence the overall conformational stability of the protein.

View larger version (20K): [in this window] [in a new window]   FIGURE 7   Drift of protein structure from the initial model. The root mean square deviation (RMSD) of all C atoms from the starting structure is shown as a function of time for simulations MDK2 (black line) and MDK3 (gray line).

Structural fluctuations

The C RMSD provides a picture of the global drift of the Kir6.2 model structure during the MDK2 and MDK3 simulations. However, it is perhaps more informative to examine the fluctuations in the structure on a residue-by-residue basis, by examining the root mean square fluctuation from its time-averaged position of each C atom (Fig. 8). One may note that the overall pattern is the same for each of the four subunits within a single simulation and that there are no major differences between the two simulations. For the -helices (M1, P, and M2), the RMSFs are low (~0.06 nm at the center of each helix). Thus these simulations do not suggest any tendency of the helices to attempt to repack, at least on a nanosecond time scale. The RMSFs are largest at the N and C termini of each chain (as one would expect) and for the M1/P loops, where they rise to ~0.3 nm. This may reflect errors in the conformation of these loops (see above). Reassuringly, the RMSFs in the filter region are low (~0.1 nm). However, this is a little higher than is seen in the immediately preceding P-helix, presumably reflecting the absence of a network of H-bonds to hold the filter rigidly in place (see below).

View larger version (26K): [in this window] [in a new window]   FIGURE 8   Fluctuations of protein coordinates. The root mean square fluctuations (RMSFs) of the C coordinates from their time-averaged values are shown as a function of residue number for the MDK2 (A) and MDK3 (B) simulations. In each graph the four lines correspond to the four subunits. Note that the residue numbers are relative to the start of the model, such that residue 1 here corresponds to residue D58 in the sequence of Kir6.2 (see Fig. 1). Thus the M1 helix corresponds to residues 11-37, the P-helix to residues 60-71, the selectivity filter to residues 73-77, and the M2 helix to residues 87-110.

An alternative view of structural flexibility during the simulations is provided by the time-dependent secondary structure of the protein (Fig. 9). This also reveals considerable fluctuations in the secondary structure (i.e., backbone torsion angles and H-bonding) of the M1/P loop. The filter (F) region shows relatively minor fluctuations in secondary structure. The M1, P, and M2 helices are maintained throughout both simulations. There are occasional localized losses of -helicity in M2 (e.g., subunits 1 and 2 of MDK2; see Fig. 9). These may be related to the flexibility seen in simulations of this helix in isolation in lipid bilayers (Shrivastava et al., 2000). It is possible that such flexibility may play a role in channel gating.

View larger version (109K): [in this window] [in a new window]   FIGURE 9   Secondary structure (analyzed using DSSP; Kabsch and Sander, 1983) as a function of time for simulations MDK2 and MDK3. The gray scale used to indicate secondary structure is given at the bottom of the figure. Black indicates -helical residues. The residue numbering is as in Fig. 8, such that the four subunits correspond to residues 1-118, 119-236, 237-354, and 355-472. The locations of the M1, P, and M2 helices and of the filter (F) are indicated to the right of each diagram.

Ion/water/pore trajectories

Visualization revealed that the K⁺ ions did not remain in their initial locations throughout the simulations. Rather, they moved within the filter and, in MDK3, within the cavity in a fashion reminiscent of that seen in simulations of KcsA (Shrivastava and Sansom, 2000). This can be seen most clearly by examination of z-coordinate trajectories (i.e., trajectories along the pore axis) of the K⁺ ions (Fig. 10). The results for MDK2 provide clear evidence for concerted, single-file motion of K⁺ ions and water molecules within the selectivity filter of the Kir6.2 model. We will focus on the two K⁺ ions in the filter (K1 initially at site S1 and K2 at site S3), on the "crystallographic" water molecule placed between the two ions at site S2 (Wx), and on two water molecules, W10391 and W10096, that transiently interact with either end of the K1-Wx-K2 column. Note that W100391 is one of the "bulk" water molecules, at the extracellular mouth of the pore, and W10096 is one of the water molecules within the central cavity.

View larger version (32K): [in this window] [in a new window]   FIGURE 10   Ion and water z trajectories for simulations MDK2 and MDK3. In each case, the coordinate system is such that the center of the selectivity filter corresponds to z = 0. Thus the selectivity filter extends from z  +0.5 nm at its extracellular mouth to z  0.5 nm. The cavity extends from z  0.5 nm to z  2 nm, and the intracellular gate is centered at z  2.5 nm. The locations of the four binding sites (S1 to S4) are shown by the thin black horizontal lines. Thick lines show the trajectories of K+ ions (black) or water molecule oxygen atoms (gray). For MDK2 (A) the trajectories of K1, Wx, and K2 and of two water molecules that transiently enter either mouth of the filter (W10391 and W10096) are shown. For MDK3 (B) the trajectories of K1, Wx, K2, and K3 are shown.

During the first ~300 ps, K1, Wx, and K2 remain at sites S1, S2, and S3, respectively. W10391 jumps in and out of a "site" just above S1, thus transiently serving as one of the two waters that solvate K1. Similarly, from ~100 to 300 ps, W10096 sits at site S4, serving as one of the two waters solvating K2. For both K1 and K2 the remaining solvation interactions are provided by Wx and by the rings of O atoms that define the sites within the filter. At ~300 ps there is a concerted transition in the locations of ions and waters. This appears to start with K2 and W10096 and is then propagated to Wx, K1, and W10391. Thus W10096 moves from S4 to the mouth of the central cavity; K2 from S3 to S4; Wx from S2 to (close to) S3; K1 from S1 to S2; and W10391 moves transiently to S1, at ~400 ps, before it is displaced by another water and diffuses into the bulk solvent. Thus, at the heart of this motion is a concerted movement of K1-Wx-K2 from S1-S2-S3 to S2-S3-S4. The latter configuration persists until ~1.5 ns, at which time K2 enters the cavity. Note that W10096 continues solvating K2 until ~500 ps, at which time it is replaced by another water from the cavity.

A similar concerted motion of K1-Wx-K2 occurs in MDK3, although this time at the start of the simulation (Fig. 10 B). Note that in this case we have not plotted the trajectories of any other waters interacting with the K⁺ ions in the filter. The trajectory of the cavity ion, K3, reveals that for the first ~200 ps it remains in the cavity close to the "focus" of the P helix dipoles (Roux and MacKinnon, 1999). Subsequently it moves "down" the cavity, but its exit appears to be blocked by the tight ring of F168 side chains that forms the ICG (see above) in Kir6.2. Thus, in MDK3 for Kir6.2 we do not see exit of K3 from the channel within 1 ns, unlike the situation with KcsA (Shrivastava and Sansom, 2000).

Flexibility of the filter

An important feature of the concerted motion of K⁺ ions and water molecules through the selectivity filter is that the filter region seems to undergo minor changes in conformation during this process. Qualitatively, this is evident if one visualizes snapshots of the filter during the progress of the concerted ion/water motion (Fig. 11). From these snapshots it is evident that there is a degree of flexibility in the filter. In particular, the filter undergoes minor changes in conformation, so that the carbonyl oxygens appear to "track" the K⁺ ions as they move along the filter. For example, note the changes in conformation around the F134/G135 peptide bond as the K1 ion at that site is replaced by a water molecule. By plotting the backbone (,) values as functions of time for the residues of the filter one can see that changes in backbone conformation in this region are correlated with movements of K1, Wx, and K2 from site to site (data not shown). It is these, generally minor, changes in backbone conformation that appear to be responsible for optimizing the ion/protein and water/protein interactions.

View larger version (42K): [in this window] [in a new window]   FIGURE 11   The selectivity filter at 250, 300, 350, and 500 ps in simulation MDK2, showing the concerted motion of ions and water. At 250 ps the occupancy of the sites is S1 = K1, S2 = Wx, S3 = K2, and S4 = W10096, whereas at 350 ps it is S1 = W10391, S2 = K1, S3 = Wx, and S4 = K2.

	DISCUSSION

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION CHALLENGES REFERENCES

Validity of the homology model

Before discussing the implications of the MD simulation results, it is desirable to attempt to assess the validity of the Kir6.2 homology model. There is a reasonable body of SDM data for Kir channels, although some of these are for members of the family other than Kir6.2 (Yang et al., 1997; Lu et al., 1999; Loussouarn et al., 1999). Cysteine scanning data for Kir2.1 implicate residues G156, N160, and L164 of Kir6.2 as forming part of the pore lining (Lu et al., 1999). Similar data for Kir6.2 implicate residues F168 and A172 as part of the pore. If one examines the locations of these side chains in the Kir6.2 homology model, all are exposed to the pore, thus supporting the model, except for A172, which lies just below the ICG. However, it should be remembered that the homology model is probably somewhat closer to the closed state of the Kir6.2 channel, whereas the cysteine scanning mutagenesis data provide information on the open state.

The difference between the open and closed forms of the channel may account for the suggestion (Minor et al., 1999) that the architecture of the Kirs differs fundamentally from that of KcsA, with the packing of M1 and M2 helices differing markedly between the two families. We believe this to be unlikely because, despite a low sequence conservation of the helices, the sequence identity in the P-region is significantly higher (~30%), suggesting that the K⁺ channels do adopt similar architectures in the transmembrane region. Moreover, given the results of Perozo et al. (1998, 1999), one would expect the interhelix interactions to differ significantly between the open and closed states of the channels. Minor et al. (1999) state in their studies on M1/M2 contacts, "[S]ince the open probability for Kir2.1 is very near 1 at negative potentials and the selection experiments require functional open channels, this arrangement most likely reflects the open state of a Kir channel."

Minor et al. (1999) demonstrate a good structural correspondence between Kir2.1 and KcsA in the residues at the M1/M2 interface, using their alignment. However, they find a poor structural correspondence in the residues at the M2/M2 intersubunit interface and lining the pore, being unable to interpret results of their selection experiments in light of the KcsA structure, leading to their conclusion that the Kirs and KcsA adopt fundamentally different structures. Nevertheless, an attempt to interpret their results using our homology model indicates not that the Kir selection experiments cannot be interpreted in terms of the KcsA structure, but rather that a satisfactory interpretation is very much dependent on the sequence alignment used, as indicated in Table 2. Significantly, the sequence alignment of Minor et al. (1999) is based on analysis of residues participating in M1/M2 interhelix interactions (hence the good structural correspondence between the residues at that interface, despite the low sequence similarity). Our alignment benefits from analysis of the physicochemical interactions of the helices with the lipid bilayer environment via MD simulations on the single helices (Shrivastava et al., 2000), in addition to including experimental evidence and results from periodicity analysis of the Kir sequences.

Overall, as the experimental data available cannot be used to rule out the validity of our homology model, and as our model is able to satisfactorily interpret results of the selection experiments, we remain reasonably confident in our hypothesis that Kirs and KcsA have the same basic architecture. Further weight is added to this hypothesis by the results of Repunte et al. (1999), who use a homology model of Kir6.2 based on KcsA to successfully account for the large difference in unitary conductance of Kir6.2 and Kir6.1.

The stability of the model during the MD simulations also lends some support to its validity. If the M1 and M2 helices were incorrectly packed one might have expected to see a greater degree of drift over the duration of the simulations. That such drift (i.e., repacking of helices) can be observed in a nanosecond simulation is suggested by the large fluctuations in conformation of the M1/P loops and by other simulations on helix bundle models of ion channels (Forrest et al., 2000).

The other main weakness of our Kir6.2 homology model that should be remembered is that the 52 N-terminal and 215 C-terminal residues are missing. It has been suggested that the C-terminal region of Kir6.2 may form an additional membrane-attached helix or helices (Doupnik et al., 1995). However, in the absence of experimental topological data to support these hypotheses, we have not attempted to include them in our model.

There is a further similarity between our Kir6.2 model and KcsA. It has been suggested (Guidoni et al., 1999) that a salt bridge between D80 and R89 of neighboring subunits of KcsA may enhance the stability of the tetramer. In Kir6.2 the model structure suggests a similar role for residues E140 and E141 of one subunit interacting with R136 of an adjacent subunit. However, the pK_A calculations indicate that the two glutamates are protonated, and so salt bridge formation does not take place during the MDK2 simulation. This aspect of the model may merit further analysis. We note that multiple interaction sites, some within the M1-P-M2 region, have been implicated in the assembly of Kir1.1 (Koster et al., 1998).

Implications of the MD simulation results

What are the implications of the MD simulation results? In terms of the fluctuations of the polypeptide backbone, the substantial fluctuations in the M1/P loop indicate that the conformation of the latter may need further refinement. However, it is conceivable that this loop may exhibit some degree of genuine dynamic disorder in the Kir6.2 structure. It is of interest that the M1/P loop in KcsA is a locus of crystal contacts in KcsA and so may be more ordered in the crystal than in the protein in a bilayer (Declan Doyle, personal communication).

It is significant that single-file concerted motion of K⁺ ions and water is seen through the selectivity filter of the Kir6.2 model. It has long been suggested, on the basis of electrophysiological data, that K channels should contain a multiion single-file pore (Hille, 1992). It is of interest that we have now observed concerted single-file motion in simulations of KcsA (Shrivastava and Sansom, 2000) and of Kir6.2 (this study). Thus this result is rather robust in that it has been seen in two different sets of simulations, with different sets of protein coordinates and different selectivity filter sequences (see below). If such single-file motion does take place, then the net motion of ions and water through the channel should be coupled, which would be detectable via measurement of streaming potentials (Aidley and Stanfield, 1996).

Ion selectivity

The model and simulations for Kir6.2 have interesting implications for the mechanism by which K channels select for K⁺ ions in favor of Na⁺ ions. It has been argued (Doyle et al., 1998) that selectivity is achieved by the filter providing an exact fit to a K⁺ ion, whereas the carbonyl oxygens are too distantly spaced for optimal interactions with a smaller Na⁺ ion. This mechanism requires the filter to have a relatively rigid geometry, which is suggested to be maintained by a network of H-bonds and extensive van der Waals contacts between the tyrosine side chain of GYG of one subunit and the W68 side chain of the WW motif in the P-helix of the adjacent subunit. This region of the structure has been likened to a layer of springs holding the pore at its required diameter. However, two aspects of the results presented in the current paper suggest that the true picture regarding selectivity may be a little more complex. First, some degree of flexibility in the selectivity filter is seen as the K⁺ ions and water molecules move from site to site. This suggests that the filter might also be able to undergo minor changes in conformation to optimally solvate Na⁺ ions. Second, the network of H-bonds that in KcsA are proposed to maintain the filter geometry cannot be formed in Kir6.2 because the tyrosine is replaced by a phenylalanine, and the tryptophan is replaced by a phenylalanine residue. This might be expected to increase the flexibility of this region. From a physiological perspective, Kir6.2 is selective for K⁺ ions. Taken together, these results suggest that a "rigid filter" model of K channel selectivity may be a little too approximate, and that the reality may involve a subtle energetic balance of ion/protein, ion water, and protein deformation energies. This is supported by a number of experimental studies. For example Silverman et al. (1998) have shown that in Kir3.1/Kir3.4 heteromultimers the tyrosine of the Kir3.1 filter may be mutated to many different residues without loss of selectivity. (These authors also comment on the replacement of the WW motif in the P-helix of KcsA by a LF motif in Kir channels). Furthermore, Yang et al. (1997) show that replacement of E138, R148 in Kir2.1 (i.e., the equivalents of E126, R136 in Kir6.2; see above) retained channel activity but resulted in loss of potassium versus sodium selectivity.

	CHALLENGES

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION CHALLENGES REFERENCES

What do the current studies tell us about homology models of K channels? First, this study shows that homology modeling can yield interesting results concerning mammalian K channels and can enable us to formulate theories and questions that are less obvious at the level of sequence alignments alone. Similar conclusions have been arrived at from homology modeling studies of the pore domains of Kv channels (Ranatunga and Sansom, unpublished results) and of 2P-domain K channels (Bartlett and Sansom, unpublished results). Homology modeling has also been applied to bacterial K transporters (Durell and Guy, 1999). It is therefore timely to apply this approach to a wider range of K channels.

The second challenge that emerges from these studies is the need to develop a model of a K channel in an open conformation. This may be possible by combining homology modeling with restraints derived from spectroscopic (Perozo et al., 1998, 1999) and SDM (Liu et al., 1997) data. Such a model could be refined by nanosecond MD simulations in a bilayer environment. Thus simulations may become a tool in more accurate modeling of membrane proteins.