The Intrinsic Flexibility of the Kv Voltage Sensor and Its Implications for Channel Gating

Protein model

The initial voltage-sensing domain model (residue 19–147) was taken from the 1.9 Å isolated VS crystal structure Protein Data Bank ((PDB) code 1ORS). Side-chain ionization states were determined based on pK_A calculations performed using the program WHATIF (12) combined with locally written code. All of the ionizable residues were in their default states (i.e., their ionization states were those expected at pH 7 based on the standard pK_A values of the residues), except for D¹⁴⁶, which was judged to have approximately equal probabilities of adopting either an ionized or neutral state (see below for details).

Setup of micelle systems

Each protein/detergent micelle was constructed using 94 DM (decyl maltoside) detergent molecules. The DM molecules were randomly rotated about their tail axes and arranged into semicircular planes, forming an expanded micelle-like torus around the isolated VS. Each detergent headgroup was approximately equidistant from its nearest neighbor. The terminal methyl groups of the DM tails were placed so that a minimum distance of 5 Å was maintained between each detergent molecule and the exterior of the isolated VS, so as to avoid unfavorable van der Waals overlaps. Similar initial protein/detergent micelle constructs have been generated in a number of other simulation studies (13,14) and have been shown to yield equilibrium structures similar to those generated by more extensive self-assembly simulations (15).

Each protein-micelle system was energy minimized and then solvated by superimposition of a box of preequilibrated SPC (16,17) water molecules, followed by removal of any water molecules too close to either protein or detergent molecules. The resultant systems were then subjected to further energy minimization. Potassium and chloride ions were added to each system by random replacement of water molecules, until a final concentration of ∼0.1 M was achieved. The resultant systems were subjected to a further stage of energy minimization. This was followed by 1 ns of restrained MD during which all of the heavy (i.e., not H) protein atoms were positionally restrained. Restrained MD runs were performed at 300 K and 368 K for each protein-micelle system (Table 1). Finally, all positional restraints were removed and production run simulations were performed for each system.

Simulation details

Simulations were performed using GROMACS 3.1 (18) (www.gromacs.org). The GROMOS96 force field (19) was used for all simulations. The parameters for DM were based on GROMOS96, supplemented with additional bond, angle, and dihedral terms to rigidify the sugar groups. Initial DM coordinates were derived from the structure of maltose when bound to a maltodextrin transport/chemosensory receptor (20) (PDB code 1ANF) combined with results from semiempirical calculations performed on the lipid chain using the CFF module in Insight II (www.accelrys.com/insight). All energy minimization procedures used <1000 steeps of steepest descent method to relax any steric conflicts generated during system setup. During restrained MD runs, harmonic restraints with force constants of 10 kJ mol⁻¹ Å⁻² were applied to all heavy atoms. Long-range electrostatic interactions were calculated using the particle mesh Ewald method (21), with a 10 Å cutoff for the real space calculation. A cutoff of 10 Å was used for the van der Waals interactions. The simulations were performed at constant temperature, pressure, and number of particles. The temperature of the protein, detergent, and solvent (waters and ions) were separately coupled using the Nose-Hoover thermostat (22,23) at 300 K or 368 K, with a coupling constant τ_T = 0.1 ps. System pressures were isotropically coupled using the Parrinello-Rahman barostat (24,25) at 1 bar with a coupling constant τ_P = 1 ps and compressibility = 4.5 × 10⁻⁵ bar⁻¹. For all simulations the time step for integration was 2 fs, and the coordinates and velocities were saved every 5 ps. The LINCS algorithm (26) was used throughout to constrain bond lengths.

Simulations were performed on a 68 node PC cluster with dual Xeon4 processors. All analyses used GROMACS and/or locally written code. Simulations were analyzed by means of principal component analysis (27): trajectories were projected along selected eigenvectors to filter the relevant motions (28). The program SWINK has been used to measure the swivel and kink angles of the helices (29). Secondary structural alignments were performed using the secondary structure matching (SSM) server (30). Molecular graphics images were generated using VMD (31) or Pymol (www.pymol.org).

Simulations in detergent micelles

The pore-forming and voltage-sensing domains of Kv channels seem to be, to some extent, independent functional modules. Thus, all K⁺ channels appear to share the same pore domain (4), whereas only Kv channels contain the VS domain. From a structural viewpoint, the isolated VS domain of KvAP forms high resolution crystals (2) and also forms a stable folded domain when in detergent micelles (2,11) and when in lipid bilayers (32). Thus, simulations of the isolated VS may provide insights into the intrinsic conformational flexibility of this domain, which in turn will inform attempts to formulate atomic resolution models of the mechanism of voltage gating of Kv channels. Furthermore, the crystal structure of the isolated VS seems to be broadly consistent with molecular physiology and spectroscopic data and thus is likely to represent a major physiological conformation of the VS when within Kv channels (33,34).

Detergents mimic a biological membrane environment and as such play an important role in a number of biochemical techniques that have been applied to the study of membrane proteins. As DM was used to solubilize both the full length KvAP molecule and its isolated VS for crystallization, we have simulated the KvAP VS in a DM micelle. The exact size of the VS-detergent micellar aggregate is unknown. Using the surface of the VS as a guide, 94 detergent molecules were needed to form a complete toroidal micelle around the protein. Relative to pure detergent micelles, increases in numbers of detergents of 10–25 have been seen for peptide-containing micelles (35). The experimentally determined size of pure DM micelles is ∼70 molecules (www.anatrace.com). Thus, 94 DM molecules in a VS/DM micelle would seem to be a reasonable estimate. The geometry of protein/detergent micelles are not known in detail. Simulation studies, both of preformed (13,14) and of self-assembled (15) protein/detergent micelles, suggest an approximately toroidal geometry for the detergent, wrapped around the protein. Therefore we inserted the VS into a preformed DM micelle, using the membrane protein aromatic residues as a guide to the location of micelle around the α-helical VS domain.

Charged residues in the VS play a crucial role in voltage sensing (36,37). Therefore pK_A calculations (see Methods) were performed on the isolated VS to ascertain the most probable side-chain ionization states when the VS was placed in a membrane mimetic environment. All of the ionizable residues were predicted to exist in their default charge states, with the exception of an aspartic acid residue (D¹⁴⁶) which was predicted to have an approximately equal likelihood for both its ionized and protonated states (its predicted pK_A was 7.9 and 5.2 in vacuo and in a membrane mimetic, respectively). We therefore have performed simulations with residue D¹⁴⁶ in each of these states (Table 1).

KvAP is from the thermophilic bacteria Aeropyrum pernix, which thrives at temperatures ∼368 K (95°C) (38). Therefore, in addition to simulations at 300 K, we have run simulations at 368 K to explore the conformational dynamics of KvAP VS at a temperature closer to its physiological operating temperature and also to explore enhanced protein flexibility at the higher temperature. In total we have performed 0.14 μs of MD simulations of the KvAP VS in a detergent micelle, a system containing ∼44,000 atoms.

Visual inspection of each VS/DM micelle system revealed that the initially ordered torus of detergent molecules in the micelle (Fig. 1 B) adjusted its shape during the course of each simulation so as to increase the solvent exposure of the ends of the helices (Fig. 1 C). This is comparable to the behavior observed in a number of other membrane protein/micelle simulations, e.g., of the β-barrel outer membrane protein OmpA (13) and of the α-helical membrane protein GlpF (14), both in micelles with octyl glucoside. The shape of the detergent component of the micelles was quantified by calculation of the eccentricity, defined as η = 1 − I I_MIN/I_AVG where I_MIN is the smallest of the three principal moments of inertia (I₁, I₂, I₃), and I_AVG is the average of all three. For a perfectly spherical object, η = 0. For the four simulations the initial eccentricity was η ≈ 0.07. At the end of the simulations, the eccentricity was η ≈ 0.19. Thus, the micelles become more ellipsoidal with time, as has also been observed with OmpA/DPC micelles (for which η ≈ 0.2—see Bond and Sansom (13)) and with GlpF/OG micelles (14).

Conformational drift

To evaluate the conformational stability of the isolated VS in these micelle simulations, we have measured the degree of drift from the initial (crystal) structure. This may be evaluated as the root mean-square deviation (RMSD) of the Cα atoms from their coordinates at the start of the simulation (Fig. 2). The general pattern is the same in all simulations—an initial jump followed by a slower rise to a plateau. This has been observed in many MD simulations of membrane proteins, both in lipid bilayers and detergent micelles. At 300 K the Cα RMSD over the final 5 ns (for all residues, i.e., including the loops between helices) is ∼2 Å regardless of the ionization state of D¹⁴⁶. This is comparable to that seen for other membrane proteins in micelles (e.g., ∼2 Å for GlpF (14) and a similar value for OmpA is obtained if one excludes the long flexible extracellular loops from the RMSD calculation for the latter protein (13)). Thus, even though the current simulations of the KvAP VS are 3–4 times longer than those of GlpF and OmpA, the conformational drift of the VS is about the same. We may therefore conclude that the VS is conformationally stable in DM micelles on a ∼40 ns timescale at 300 K.

At the higher temperature, the RMSDs for both simulations are substantially higher (∼4 Å). (This is comparable to the figure for OmpA at 300 K if the flexible loops are included in the calculation.) However, if only the α-helical regions are taken into consideration, the final Cα RMSDs are ∼3 Å. This indicates that at 368 K the greater structural drift is largely, but not entirely, due to conformational changes occurring in the loop regions. Again, the protonation state of D¹⁴⁶ does not seem to influence the conformational drift at the higher temperature.

In summary, from consideration of conformational drift we can see that removal of the VS from its crystallographic interactions with the Fab molecule and embedding it in a detergent micelle does not lead, on a ∼40 ns timescale, to a substantial change in conformation even at an elevated temperature. By way of comparison, chymotrypsin inhibitor 2 has been observed to unfold on a ∼25 ns timescale in simulations at 373 K (39). Taken together, this suggests that the isolated VS domain is intrinsically stable in a micelle environment.

Secondary structure and flexibility

A more detailed analysis of the simulations reveals the relative flexibility of the different regions of the VS. Analysis of the secondary structural elements as defined by DSSP (40) (Fig. 3 A) showed that helices S1, S2, and S4 were stable over the course of each simulation, regardless of the temperature or protonation state of D¹⁴⁶. In contrast, the S3 helix (which in the crystal structure is distorted to form distinct S3a and S3b helices, Fig. 3 D) and the S3/S4 loop show enhanced flexibility at the higher temperature. In particular, S3a unfolds at 368 K, regardless of the protonation state of D¹⁴⁶. (There is also some loss of helicity from the N-terminus of S4 at 368 K, but only if D¹⁴⁶ is ionized.) Thus, the elevated temperature simulations suggest a degree of conformational instability for S3a. Interestingly, a consensus TM helix secondary structure prediction analysis of KvAP (41) suggested that S3a may not form a canonical TM helix.

It is possible to explore in more detail the relative flexibility of different regions of the VS by calculation of the root mean-square fluctuation (RMSF) of each Cα atom from its time averaged coordinate (Fig. 3 C). The overall pattern for each simulation is as anticipated, with higher RMSFs for the loops and lower RMSFs for the cores of the α-helices. However, there are some more significant details. First, the S2/S3a loop and the S3a region show a significantly higher RMSF at 368 K than at 300 K. This is consistent with the secondary structure analysis and suggests that S3a is unstable at the higher temperature. The second significant detail is that the S3b-S4 loop shows a relatively high RMSF in each of the simulations.

A comparison with the corresponding crystallographic temperature factor (B-value; Fig. 3 B) profile is informative. Of course, one should guard against overinterpretation of comparisons of crystallographic B-values and simulation RMSFs. However, a striking difference in the overall profiles is evident, namely that the crystallographic B-value for the S3b/S4 loop is low relative to the remainder of the protein, whereas the RMSF for this region is high relative to the protein as a whole. As shown in the inset in Fig. 3, the S3b-S4 loop corresponds to the primary interaction between the VS and the Fab molecule. Thus it would seem that either the intrinsic flexibility of the S3b/S4 loop is suppressed within the Fab/VS crystal and/or the local conformation of this region in the crystal may be altered relative to its equilibrium structure in the micellar environment. It is of possible functional significance that this region, the S3b-S4 loop, corresponds to the tip of the voltage-sensing paddle that is suggested to be capable of changing its conformation/orientation relative to the membrane in the paddle model of voltage gating (11).

The results of the analysis of secondary structure and flexibility may be summarized with reference to Fig. 3 D. Comparison of the crystallographic structure with those at the end of the VS300 and VS368 simulations for the S2–S3 region shows that there is a local change in conformation in S3a and a local distortion of the C-terminus of S3b. Thus, the simulations reveal that the S3a and S3b helices possess a degree of conformational flexibility. As we will see, this is important in allowing the S4 helix to move relative to the remainder of the VS domain, a process that is in turn of relevance to voltage gating in the channel.

The crystal structure of the isolated VS domain of KvAP superimposed on the corresponding VS domain of Kv1.2 reveals that the two structures are similar (RMSD 2.9 Å). However, they diverge most notably in the region of the S3 helix (Fig. 4). Thus, the Kv1.2 crystal structure has captured the S3 helix in a different conformation from that of KvAP, which may further support the notion of a flexible S3 helix.

S4 helix movements and flexibility

The S4 helix is believed to move relative to the remainder of the channel in some way that couples voltage-sensing to channel gating. It is therefore useful to examine the simulations from the perspective of the intrinsic flexibility of the S4 helix, and in particular to see whether or not it behaves as a rigid helical rod.

Simple visual examination of snapshots of S4 suggests a degree of hinge bending may occur in this helix about the central I¹³⁰ residue (Fig. 5 A). To explore this aspect of S4 dynamics in more detail, we employed principal components analysis (PCA) to reveal the major motions of the helix. We combined PCA with analysis of the motions of S4 in terms of helix kink and swivel. Analyzing TM helix distortions in terms of kink and swivel angles has been employed to explore proline-induced distortions of TM helices in x-ray structures of membrane proteins (29) and more recently in analysis of M2 helix hinge bending as the basis of channel gating in simulations of inward rectifier K⁺ channels (42). If we take the motions of S4 represented by the first two eigenvectors (together these account for ∼50% of the motion observed in the simulations) and analyze these in terms of kink and swivel angle, we see clear evidence for both swivel motion and, to a lesser extent, kink motion about the center of S4 (Fig. 5, B and C). The central hinge point seems to be around residue I¹³⁰, i.e., in the center of the S4 helix. Thus, the S4 helix does not behave as a rigid rod but rather as two semirigid α-helical segments connected by a central molecular hinge. We note that the flexibility observed in the S3b/S4 loop may facilitate such motions of S4.

Water penetration

There have been a number of suggestions that the S4 helix sits within a water-filled crevice, which therefore focuses the TM voltage difference on a narrow central region of the S4 helix (4,10). This is evidently of some importance for possible voltage sensing/gating models of Kv and related channels (43). We have used water penetration within the protein/detergent micelles as a measure of the degree of crevice formation in the VS. As can be seen (Fig. 6 A), a number of water molecules do infiltrate the VS, forming hydrogen bonds with, e.g., R¹¹⁷, R¹²⁰, R¹²³, and R¹²⁶ of the N-terminal half of S4 (Table 2). Thus, our simulations indicate that hydrogen-bonding interactions with such waters may help to stabilize the positively charged S4 helix within a bilayer-like environment. Preliminary analysis of simulations of the VS in a lipid bilayer support this claim (Z. A. Sands and M. S. P. Sansom, unpublished results). Of course, in the intact channel the VS is packed next to the core pore-forming domain of the protein. However, it is likely that the micelle mimics the environment provided by the remainder of the protein in the intact channel.

Analysis of the radius profile along the length of the presumed crevice (Fig. 6 B) reveals that there is a central constriction which, averaged across the simulation, is significantly narrower than the radius of a water molecule (∼1.4 Å). Thus the voltage drop across the membrane must occur over a region at least ∼12 Å thick around S4. The central region corresponds to a salt bridge between a side chain in S2 (D⁶²) and a side chain in S4 (R¹³³; see Fig. 6 C), which is maintained throughout the course of each simulation. This salt bridge holds S2 and S4 in close proximity and prevents water from passing from one side to the other of the crevice. It is noteworthy that the central hinge in S4 (at residue I¹³⁰) is just one turn of the helix before R¹³³. Thus the two halves of S4 are connected by a hinge that fits within the narrow central region of the water penetrable crevice.

Detergent interactions

It is of some interest whether the gating charges on S4 are able to favorably interact with the micelle environment or whether they remain buried within the interior of the VS. In the paddle and the twisted S4 models, the charged S4 helices are significantly exposed to the surrounding lipid environment during their translocation from one side of the membrane to the other. Therefore favorable interactions between the gating charges and the lipid environment could significantly lower the energy barrier to channel activation. Analysis of the hydrogen-bonding patterns reveals that the side chains of Arg's 117, 120, 123, and 126 are able to form multiple hydrogen bonds with the uncharged headgroups of the decyl maltoside detergent molecules (Table 2).