Proteins, including ion channels, often are described in terms of some average structure and pictured as rigid entities immersed in a featureless solvent continuum. This simplified view, which provides for a convenient representation of the protein's overall structure, incurs the risk of deemphasizing important features underlying protein function, such as thermal fluctuations in the atom positions and the discreteness of the solvent molecules. These factors become particularly important in the case of ion movement through narrow pores, where the magnitude of the thermal fluctuations may be comparable to the ion pore atom separations, such that the strength of the ion channel interactions may vary dramatically as a function of the instantaneous configuration of the ion and the surrounding protein and pore water. Descriptions of ion permeation through narrow pores, which employ static protein structures and a macroscopic continuum dielectric solvent, thus face fundamental difficulties. We illustrate this using simple model calculations based on the gramicidin A and KcsA potassium channels, which show that thermal atomic fluctuations lead to energy profiles that vary by tens of kcal/mol. Consequently, within the framework of a rigid pore model, ion-channel energetics is extremely sensitive to the choice of experimental structure and how the space-dependent dielectric constant is assigned. Given these observations, the significance of any description based on a rigid structure appears limited. Creating a conducting channel model from one single structure requires substantial and arbitrary engineering of the model parameters, making it difficult for such approaches to contribute to our understanding of ion permeation at a microscopic level.
For the sake of simplicity, proteins often are pictured as rigid entities, corresponding to some average structure, immersed in a featureless solvent continuum. This simplified representation, which may be convenient and useful in some circumstances, should be used with extreme caution when considering ion permeation through narrow molecular pores. It long has been known that proteins have a rather fluid, dynamic structure with rapid conformational fluctuations (Cooper, 1976). In fact, proteins are composed of discrete atoms, which constantly are undergoing thermal fluctuations, from the rapid (picosecond) vibrations, through slower (multi-nanosecond) global reorientations and side chain isomerizations, to long timescale (microsecond to second) conformational changes (Karplus and McCammon, 1981). The reality of these fluctuations is evident in the Protein Data Bank (PDB) structure database, which reports not only a set of fixed coordinates but also the temperature B factors that denote (among other factors) the protein's thermal fluctuations. Understanding how protein dynamics, of all magnitudes and timescales, can influence protein function remains a significant challenge in modern biophysics (Karplus and Petsko, 1990).
In addition to thermal fluctuations, the finite size of solvent water molecules gives rise to structural effects on length scales smaller than 8–10 Å (Israelachvili and Wennerström, 1996), which means that errors may be introduced by approximating water as a continuum solvent. Therefore, even though the simplified views of rigid proteins and continuum solvent can be exceedingly useful in some cases, it is important to know when (and how) these notions break down and become inappropriate. These issues become particularly acute when trying to describe ion permeation through narrow single-file pores, where even small displacements in the atomic coordinates would be expected to have large energetic consequences.
Atomic motions have been shown experimentally to be important for the function of many proteins. Among the first examples of this is the observation by nuclear magnetic resonance (NMR) and computer simulation of side-chain flipping in the bovine pancreatic trypsin inhibitor (BPTI; see Northrup et al., 1982, and references therein). Hemoglobin and myoglobin also offer striking examples, as the crystal structures (Perutz and Mathews, 1966; Takano, 1977) provide no obvious path for O2 moving to the heme pocket, which led to the suggestion that the protein undergoes small fluctuations that allow for transient paths for O2 access (Perutz and Mathews, 1966). This suggestion was supported by the experimental observation that myoglobin undergoes thermal fluctuations (Frauenfelder et al., 1979) and by molecular dynamics simulations that identified a path for O2 movement to the heme pocket (Case and Karplus, 1979). Likewise, numerous simulations have reported significant atomic fluctuations in ion channels such as gramicidin A (gA) (Mackay et al., 1984; Roux and Karplus, 1991) and the KcsA potassium channel (Guidoni et al., 1999; Allen et al., 2000; Åqvist and Luzhkov, 2000; Bernèche and Roux, 2000; Shrivastava and Sansom, 2000; Noskov et al., 2004). Thermal fluctuations long have been acknowledged to be important for channel gating, whereas their quantitative effects on ion–protein interactions and their importance for ion permeation, though recognized, have received less attention (Hille, 1992).
In fact, though everybody “knows” that ion channels, like all other proteins, are flexible macromolecules, and that water molecules are not infinitesimally small, selective ion permeation through narrow pores tends to be discussed, at least at the conceptual level, with the underlying assumption that selective ion solvation implies that channels are quite “rigid” (Hille et al., 1999). Because of its conceptual simplicity, this picture of ion permeability and selectivity has influenced the interpretation of experimental results and the design of many computational models. The exact range of validity of this assumption/approximation and the implications of its breakdown remain poorly understood, however, in cases where the dimensions of the pore are comparable to the ion sizes, as is the case for the cation-selective, voltage-dependent channels (Hille, 1971, 1972, 1973). While some computational approaches, such as molecular dynamics (MD) simulations (Mackay et al., 1984), naturally incorporate the flexibility of macromolecular channel structures, extensive efforts have been spent developing and engineering computational models of ion permeation through narrow pores starting from the approximation of a unique fixed rigid structure with continuum solvent (Kurnikova et al., 1999; Allen and Chung, 2001; Hollerbach et al., 2001; Mashl et al., 2001; Burykin et al., 2002; Edwards et al., 2002; Nadler et al., 2003). The limitations inherent in this approximation are increasingly being recognized and efforts are made to address the worst problems (Mamanov et al., 2003). It also has been stated explicitly that the structure of channels are truly rigid and fixed, i.e., fluctuating by <0.1 Å (Eisenberg, 2003). It thus becomes important to examine critically how narrow pores behave at the microscopic level and on this basis adopt a perspective of ion permeability and selectivity that is in better accord with the physical and chemical characteristics of proteins and solvent.
The purpose of this communication, therefore, is to highlight the significance of protein flexibility and solvent representation for ion permeation through narrow channels. We will conclude that the expected atomic fluctuations (of sub-Angstrom order) in the positions of the pore-lining atoms have major consequences for the energetics of ion permeation. Thus, any description of ion permeation based on a single, rigid structure should be viewed with caution. More generally, the root-mean-square (RMS) thermal fluctuations set a minimal length scale for mechanistically significant differences in atomic sizes and positions.
RELATION BETWEEN FLEXIBILITY AND THERMAL FLUCTUATIONS
What is the magnitude of thermal fluctuations in proteins and what is their importance? To put this question in perspective, let us compare and contrast atomic fluctuations in proteins and various solid materials. Thermal motions typically are reported in terms of the RMS fluctuations of the three-dimensional atomic coordinates
The intrinsic flexibility of a protein, and how it responds to external perturbations, is fundamentally related to the magnitude of atomic thermal fluctuations. To illustrate this important concept, we consider a single atom fluctuating in one dimension around a mean position x with mean square amplitude 〈Δx2〉. At the simplest level, the atom in the protein can be pictured as being attached to its average position by a harmonic spring with effective force constant Keff, which is determined by the average fluctuations (Parak and Knapp, 1984; Zaccai, 2000):
where 〈Δx2〉 = 〈Δr2〉/3 in the case of isotropic motion in 3D. The concept of an effective thermal restoring constant, and its implications for protein flexibility, has been previously used to analyze experimental data (Parak and Knapp, 1984; Zaccai, 2000; Gabel et al., 2002). If the fluctuations are small, the atom behaves as if it were held near its mean position by a very stiff spring. As the RMS of the fluctuations increases, this effective spring becomes weaker.
Because the atom position fluctuates, the action of an external force Fext, such as that due to the electric field arising from a nearby ion, will perturb the atom and shift its position. From Hooke's law, the atom initially at x0 will reach a new equilibrium position x0′ when the restoring force from the effective spring Keff will be equal and opposite to the external perturbative force.
That is, the average displacement of the atom in response to an external perturbative force is fundamentally related to the magnitude of the equilibrium thermal fluctuations. While the approximation expressed in Eq. 2 eventually will break down, it remains valid for small (sub-Angstrom) displacements from x. Furthermore, although Eq. 2 was derived without explicit consideration of thermal fluctuations, it is valid also when all thermal fluctuations are taken into account (see PP1appendix).
Because the effective thermal spring constant Keff is an equilibrium property defined via the atomic RMS fluctuations, Eq. 2 is independent of kinetic factors (see PP1appendix). Therefore, even though the atom in the unperturbed protein may oscillate around its mean position at some high frequency, ω, the magnitude of the average displacement in response to some external perturbing force depends only on the magnitude of the fluctuations determined from an equilibrium averaging as described by Eq. 2.
To illustrate the importance of protein flexibility in the context of (the strong) ion–protein interactions, we consider a one-dimensional system comprising one K+ interacting with a single carbonyl (C=O) group. We assume that the oxygen atom initially is at a distance x from the fixed ion. The protein is represented by an effective harmonic restoring force, which limits the ability of the carbonyl group to move in response to the field around the cation. The displacement of the carbonyl group in response to its interaction with the ion is plotted in Fig. 1 as a function of the initial separation between the ion and the oxygen.
Materials with different various RMS fluctuations are considered: 0.1 Å, similar to a covalently bonded solid, e.g., silicon and diamond (Lu et al., 1993); 0.2 Å, a fairly rigid crystalline solid, e.g., LiF (Butt et al., 2003); 0.3 Å, a more flexible crystalline solid, e.g., KCl (Butt et al., 2003); and 0.6 Å, corresponding to the RMS observed for many proteins (Ringe and Petsko, 1986; Gabel et al., 2002) including the value obtained in the gA simulations (described below). For comparison, perfectly rigid (RMS = 0 Å) and flexible (RMS = ∞) materials are shown as the limiting cases. Rigid solids do not respond significantly to the perturbation. In striking contrast, a protein-like “soft” material (0.6–1.0 Å RMS) behaves almost like an ideally flexible material for sub-Angstrom deformations. For example, if the ion and carbonyl group initially are separated by 3.0 Å, the carbonyl will deflect by 0.39 Å, returning almost to the contact distance (of 2.58 Å). These considerations show that the (atomic) thermal fluctuations in proteins easily allow for sub-Angstrom displacements in response to strong imposed forces, such as those due to the electrostatic field around an ion. Because the elastic energy depends on the square of the strain, the response would be only prohibited for materials having RMS fluctuations close to 0.1 Å, much less than the magnitude of observed protein fluctuations.
Whereas the conformational stability of many proteins implies that proteins can be approximated as rather rigid bodies on the multi-Angstrom to nanometer length scales, the ease with which a protein structure can distort on the sub-Angstrom level has serious implications for understanding the mechanisms underlying ion selectivity in cation-selective channels, such as potassium channels. According to the simple arguments laid out above, assuming RMS fluctuations in the range of 0.75–1.0 Å (corresponding to the results of Zhou et al., 2001, and calculations of Bernèche and Roux, 2000, and Noskov et al., 2004), the displacement of a single carbonyl group by 0.4 Å (corresponding to the size difference between K+ and Na+) will incur an elastic strain energy of only a fraction of a kcal/mol. In contrast, the energy of interaction between a K+ and a single N-methylacetamide (NMA) molecule differs from the energy of interactions between a Na+ and a single NMA by 5–10 kcal/mol (Bernèche and Roux, 2002, and references therein). That is, the flexibility of proteins allows for small (local) displacements in response to the perturbing influence of an ion. In order for the strain energy to become comparable to the electrostatic energy, the channel would have to be about as rigid as diamond or silicon!
These considerations indicate that explanations of ionic selectivity in terms of the precise geometry of a (rigid) pore perfectly adapted to snugly fit K+, but not Na+ (Doyle et al., 1998; Hille et al., 1999; Yellen, 2002), are incomplete and unlikely to be correct. This conclusion raises fundamental questions about the physical mechanism enabling potassium channels to maintain a high selectivity and compels us to seek alternatives to the snug-fit hypothesis. Clearly, the spatial organization of the pore-lining groups is important (Doyle et al., 1998; Zhou et al., 2001), as it provides for an environment that both energetically and kinetically is well adapted to solvate (and desolvate) the permeating ions. It has long been recognized that channels use a combination of geometric, chemical, and energetic factors to select among different ions, and that the size of the ion and a narrow ion channel protein pore must be similar in order to allow for selective ion permeation (Mullins, 1960; Hille, 1973); but descriptions of ion selectivity must be consistent with the dynamic and atomistic nature of proteins. This entails considering alternative perspectives that also include the properties of the microscopic interactions and their influence on solvation and selectivity (Eisenman, 1961; Eisenman and Horn, 1983; Yamashita et al., 1990). Developing further the original ideas of G. Eisenman, an attractive possibility is that selectivity might arise locally, from the intrinsic dynamical and electrostatic properties of the microscopic interactions between the cation and the carbonyl groups (Noskov et al., 2004).
IMPLICATIONS FOR CALCULATING ENERGY PROFILES
In the following, we consider two ion channels (gA and KcsA) to further illustrate quantitatively the importance of protein flexibility for understanding ion permeability. To focus the presentation, we will highlight the limitations of rigid pore models in studies of ion permeation. We further will show that it is necessary to treat the pore water molecules as molecular entities, as opposed to a continuum solvent.
Gramicidin A Channels
The amount of experimental information together with the channel's relative simplicity (Andersen and Koeppe, 1992; Busath, 1993; Koeppe and Andersen, 1996) have motivated numerous computational studies of ion permeation through gA channels (Roux, 2002), including several studies that employ rigid protein models and approximate electrostatic representations (Kurnikova et al., 1999; Edwards et al., 2002; Nadler et al., 2003).
The most fundamental limitation in setting up a model in which the solvent is represented as a continuum dielectric, with the goal of studying ion permeation through a narrow pore, is the obvious problem of trying to define macroscopic dielectric constants for microscopic length scales (on the order of Å). This problem is exemplified well in the case of gA channels, which clearly demonstrate the fundamental and well-documented difficulties associated with attributing a uniform macroscopic dielectric constant (of any value) to water molecules organized in single file along a narrow pore (Partenskii and Jordan, 1992; Dorman and Jordan, 2004). The simplest (and most commonly used) choice is to assume a bulk-like value for the dielectric constant of the solvent region, such that εw = 80. At the microscopic level, however, there is no justification for this or any other fixed value (Partenskii and Jordan, 1992; Dorman and Jordan, 2004), a situation that already was perceived by K. Wilson and collaborators (Mackay et al., 1984). Similarly, there is uncertainty concerning the protein dielectric constant εp. Experimental determinations on dry polypeptides range from 4 to 10 (Tredgold and Hole, 1976), and a value of 4 was deduced from the effects of side chain dipoles on the conductance of gA channels (Koeppe et al., 1990). A priori calculations suggest values ranging from 2 to 5 depending on the situation (Simonson and Brooks, 1996; Pitera et al., 2001; Schutz and Warshel, 2001), though a value of 2 has been used in many theoretical studies. One may note that the value of εp has a small effect for large and medium-sized pores (Bastug and Kuyucak, 2003), where one can approximate the pore water to resemble bulk water; but a large effect in the case of narrow pores (Corry et al., 2003), where the application of continuum models is questionable anyway.
Nevertheless, assuming that one is willing to use continuum dielectric approximations, as has been done in numerous studies of ion permeation, a second difficulty arises from the need to define the dielectric boundary separating the protein and solvent regions, i.e., to decide which region of space that will be assigned the different dielectric constants (εw or εp). Conceptually, the dielectric boundary corresponds more or less to the van der Waals envelope of the protein, which defines the volume from which the water molecules are physically excluded. For example, if all the atoms of the protein are represented as hard spheres, the solvent-accessible region corresponds to the region of space that is located outside any of those hard spheres. The advantage of this procedure is that it is operationally well defined and very easy to implement. The disadvantage is that it ignores the basic fact that the water molecule itself has a physical dimension and that the solvent is not a continuous medium of infinitesimal points (able to fill the tiniest cavities and cracks in the interior of the protein). Therefore, a more sophisticated method consists of also representing the water molecule as a hard sphere. The dielectric boundary between the solvent and the protein then is defined by “rolling” the chosen water probe onto the protein, mapping all the regions of space that are physically accessible to the solvent. The procedure is illustrated schematically in Fig. 2; it is the standard procedure used in applications of Poisson-Boltzmann theory to globular proteins (Davis et al., 1991; Sitkoff et al., 1994; Nina et al., 1997). This procedure is intrinsically more physical because it avoids the spurious effects of assigning a high dielectric value of εw to small spatial regions in the interior of the protein that are not even able to contain a single water molecule.
We performed Poisson equation solutions with the PBEQ module (Nina et al., 1997) of CHARMM (Brooks et al., 1983) for a single gA ion channel structure, PDB:1GRM (Arseniev et al., 1986), embedded in continuous dielectric membrane and water, to demonstrate the difficulties in assigning dielectrics. A schematic, resembling the model gA system, is illustrated in Fig. 3.
Fig. 4 A shows the effect of probe size on the assignment of dielectric constants throughout the pore of the PDB:1GRM gA structure. Though water can permeate gA channels (Rosenberg and Finkelstein, 1978), a water probe of 1.4 Å radius will not fit throughout most of the gA pore, an artifact that results from choosing a single rigid set of atomic coordinates (see Fig. 3). This contradiction illustrates well the significance of thermal fluctuations of the protein and reveals (again) the problems associated with the concept of a rigid structure. As a consequence, most of the lumen will not be assigned the high dielectric constant of water. Potential energy profiles (also referred to as the static field or fixed charge field) were obtained from single solutions to Poisson's equation, in the absence of an ion, and reported as the electric potential along the z axis, multiplied by the unit charge e. Fig. 4 B shows that the effect of the probe size on the electrostatic potential energy may be as much as 8 kcal/mol, emphasizing the difficulties/uncertainties associated with using a continuum solvent model in a narrow pore. The significance of defining dielectric boundaries where one has difficulties fitting even a single water molecule is unclear.
In the procedure described above, channel atoms, ions, and water molecules were represented as simple hard spheres. Additional difficulties therefore arise when choosing the particular value of the radius that should be assigned to the various atom types. Those atomic radii are required for the construction of the dielectric boundary, and their values have a very significant impact on the results of any model based on macroscopic electrostatics. In the present calculations, we employ the atomic radii of Nina et al. (1997), which were optimized from MD free energy calculations with explicit water molecules to accurately reproduce the solvation free energy of all amino acids. In this library, the backbone atomic radii are 1.52 Å for O, 2.04 Å for C, 2.23 Å for N, and 2.38–2.86 Å for Cα. Furthermore, a water probe radius of 1.4 Å was used to determine the dielectric boundary. In contrast, Nadler et al. (2003) arbitrarily used a unique value of 1.5 Å for all nonhydrogen atoms, providing no justification to explain why all atoms can be treated as having the same size. Edwards et al. (2002) used atomic radii derived from atomic contact distances (Li and Nussinov, 1998). Those radii are somewhat smaller than the radii of Nina et al. (1997), and their ability to accurately reproduce solvation free energies is not known. In addition, those previous studies constructed the protein–solvent dielectric boundary without using a water probe (i.e., assuming that water molecules are infinitesimal “points”). The specific choice can lead to a factor of two difference in the calculation of the dielectric self-energy (also referred to as the reaction-field energy, electrostatic barrier, image charge repulsion, and dielectric boundary force) contributions to the energy profile that are reported below, again reinforcing the limitations of the continuum solvent approach.
Once a model based on a unique rigid channel structure and continuum solvent has been adopted, the results are extremely sensitive to the slightest variations in atomic coordinates. Fig. 5 A shows the potential profile along the channel z-axis for the available high resolution NMR structures of the single-stranded β6.3-helical bilayer– or micelle-incorporated gA dimer: PDB:1GRM (blue; Arseniev et al., 1986); PDB:1MAG (green; Ketchem et al., 1997); and PDB:1JNO (red; Townsley et al., 2001).
These profiles were calculated without using a water probe (i.e., probe radius 0.0 Å) to define the dielectric interface and differ by as much as 6 kcal/mol. Fig. 5 B shows that profiles obtained using a 1.4 Å probe differ by as much as 14 kcal/mol. Electrostatic self energy profiles (Fig. 5, C and D) were calculated by solving Poisson's equation for each position of a K+ ion, placed at regular points along the channel z-axis, with the protein charges set to zero. These self energy profiles reveal the sensitivity to the precise shape of the protein–water dielectric boundary when an ion is brought inside the narrow pore from the bulk. For probe radii 0.0 Å and 1.4 Å, the self energies differ by up to 2.2 kcal/mol and 2.7 kcal/mol, respectively. Moreover, the total energies (Fig. 5, E and F) show that there is no cancellation of errors between the potential and self energy components. Of particular importance, the PDB:1GRM and PDB:1JNO coordinates correspond to essentially the same structure (based on NMR observables and backbone pitch; Woolf and Roux, 1996; Allen et al., 2003), yet they yield quite different potential profiles. The problem only becomes accentuated by noting (Fig. 5 F) that the total energies, calculated using any of the rigid structures, are significantly above the one-dimensional potential of mean force (PMF) profiles obtained from MD simulations with explicit solvent (Allen et al., 2004).
But gA channels are not rigid. Based on MD simulations, the RMS fluctuations in carbonyl oxygen atomic coordinates vary between 0.4 and 0.65 Å (Fig. 6). The fluctuations are lowest for residues 6–8, in general agreement with experimental solid-state NMR studies (Lee et al., 1995), and larger near the COOH-ethanolamide termini where the cations spend most of their time (thus the value of 0.6 Å used in the earlier illustration) and the formyl-NH2 termini where the two subunits join. These fluctuations in the absence of an ion, also are evidenced by 7–13° fluctuations in backbone ϕ and ψ dihedral angles, in agreement with experimental estimates (North and Cross, 1995). The corresponding flexibility of the protein (and response to an imposed force) is revealed by MD simulations with a K+ placed near the entrance of the channel, which causes the backbone ϕ and ψ dihedral angles to deviate by 2–8°, relative to the empty channel; estimates that are within the range deduced experimentally (Tian et al., 1996).
These small, Angstrom-scale fluctuations in structure that result from thermal motions produce large variations in the energetics of ion permeation, when used in models based on any one rigid structure and continuum electrostatics. To illustrate the effect of thermal fluctuations on ion energetics, we also report Poisson solutions for multiple frames from an MD trajectory. Simulations were initiated with the PDB:1JNO gA structure, as described previously (Allen et al., 2004). Each sample was oriented such that the channel axis coincides, as closely as possible, with the fixed z-axis of the system. A total of 200 samples were taken from 4 ns of MD simulation (six of which were excluded because protein atom centers came within 1.4 Å of the fixed z axis). The variations in the potential profile are as much as 15 kcal/mol with no probe (Fig. 7 A), and 39 kcal/mol with a 1.4 Å probe (Fig. 7 B).
The self-energy profiles calculated with 0.0 and 1.4 Å probes are plotted in Fig. 6 (C and D). They vary by up to 3.6 kcal/mol when calculated with no probe, and by up to 6.0 kcal/mol with a 1.4 Å probe. Similar variations in continuum energetics (of up to 20 kcal/mol) were calculated from MD trajectories with an ion included within the pore (Mamanov et al., 2003). These variations are an order of magnitude too large to permit meaningful analysis of ion permeation. In contrast, when the atomic fluctuations, which are intrinsic features of any molecular system, are integrated into the equilibrium PMF calculated from MD simulations (Allen et al., 2004), the resulting PMF (Fig. 5, E and F) is well converged and well defined because it corresponds to the reversible work done by the mean force. The PMF is well behaved because it is the force that is properly Boltzmann-weighted averaged over a large number of configurations for each ion position; hence, the name potential of mean force. A properly constructed PMF thus becomes a critical first step in the formulation of nonequilibrium permeation models (Roux et al., 2004). In particular, the PMF is consistent among modern all-atom force fields, and integrating over 1 or 2 ns of simulation does not alter the profile by more than a fraction of a kcal/mol (unpublished data). Lastly, one should stress that a key structural feature of the equilibrium PMF such as the binding sites near the entrance of the channel (Fig. 5, E and F) cannot be reproduced by the continuum model, regardless of the chosen coordinates or probe size, as shown in Fig. 5.
KcsA Potassium Channels
As a second example of a narrow pore, we consider the KcsA channel. Since the structure of this channel was reported by Doyle et al. (1998) and at higher resolution by Zhou et al. (2001), KcsA channels have been the focus of many simulation studies, including continuum electrostatic approximations that do not consider the consequences of protein thermal fluctuations (Allen and Chung, 2001; Mashl et al., 2001; and Burykin et al., 2002). As would be expected from the preceding analysis of the importance of thermal fluctuations in gA channels, atomic fluctuations also are important in KcsA. To illustrate the inherent problems faced by models based on a rigid structure and continuum electrostatics, we compare the Poisson solutions obtained with the two available X-ray structures.
Fig. 8 shows the potential profile calculated for the 3.2 Å PDB:1BL8 (Doyle et al., 1998) and 2.0 Å PDB:1K4C (Zhou et al., 2001) resolution KcsA potassium channel structures. The positions of the carbonyl oxygen atoms in the selectivity filter differ by an RMS of 0.7 Å. This difference is similar to the RMS thermal fluctuations of the selectivity filter backbone atoms (Zhou et al., 2001), indicating that the two structures are essentially identical (Bernèche and Roux, 2000). The electrostatic potential profile along the channel axis through the selectivity filter with a probe of radius 0 Å (no probe), is plotted in Fig. 8 A. The potential profiles differ by as much as 15 kcal/mol. With a probe radius of 1.4 Å (Fig. 8 B), these differences can be as large as 47 kcal/mol.
As was the case for the gA channel, fundamental difficulties arise when attempting to define the dielectric interface in the narrow selectivity filter of the KcsA channel. In both the 3.2 Å and 2.0 Å resolution structures, the 1.4 Å water probe, and even a smaller 0.8 Å probe, could not fit at many positions within the selectivity filter. For example, the Val76 carbonyl oxygen atoms on opposite subunits are separated by <4.6 Å in 1K4C, corresponding to an approximate pore diameter of 1.8 Å (pore radius of 0.9 Å). This is too narrow for either a water molecule or a K+ to pass, implying that the selectivity filter must be flexible to permit conduction. Comparison of Fig. 7 (A and B) reveals that the effect of the probe size on the potential profile for the KcsA structures can be as large as 83 kcal/mol. Thus, even the simple task of assigning the dielectric regions cannot be resolved. As in the case for the gA channel, such extreme sensitivity to a unique structure is avoided in the PMF calculated from MD (Bernèche and Roux, 2001, 2003).
Though the simplified representation of proteins in terms of their average structures (neglecting the thermal fluctuations) provides invaluable insight, it also causes major difficulties for the analysis of ion movement through narrow pores. Considering only the average structure of an ion channel leads to serious problems because it ignores not only the thermal fluctuations but also the distortion of the average structure in response to the presence of an ion, both of which are essential ingredients in a proper description of selectivity (Noskov et al., 2004) and channel function generally (Immke et al., 1999; Boccaccio et al., 2004).
The reason why there are problems associated with an analysis based on any unique, fixed structure may be illustrated in simple terms in the case of the electrostatic interactions, which plays a dominant role in permeation and selectivity. Though Coulomb's law between two charges varies instantaneously as 1/r, then 〈1/r〉 ≠ 1/〈r〉 whenever r undergoes (thermal) fluctuations. Moreover, the inequality is independent of the chosen timescale. Using a unique fixed structure may be an acceptable approximation in the case of very wide pores, such as Escherichia coli OmpF porin (Im and Roux, 2002) where 〈1/r〉 ≈ 1/〈r〉; but the situation is very different when the dimensions of the pore are comparable to the ion sizes. This is why, despite the sub-Angstrom fluctuations being much faster (picoseconds) than the timescale for ion conduction (nanoseconds), a representation based on a single static structure, even one corresponding to the average position of the atoms, cannot provide a valid model for ion permeation through a narrow pore. Experimental results on thermal fluctuations in proteins (Frauenfelder et al., 1979; Parak and Knapp, 1984; Ringe and Petsko, 1986; Gabel et al., 2002), the ion-induced reorganization of the KcsA pore (Zhou et al., 2001) as well as simulation-based calculations (including those presented here), show that proteins are flexible on sub-Angstrom length scales. This means that the protein's response to perturbative forces, such as those arising between permeant ions and the charged or dipolar groups lining the pore, cannot be neglected in narrow ion channels.
The present analysis indicates that ionic selectivity in K+ channels cannot be a property of a rigid channel structure, in which the selected ion “fits” particularly well, such that the oxygens lining the pore coordinate the selected cation with sub-Angstrom precision. Although the conserved residues surrounding the selectivity filter are essential for the overall stability of the three-dimensional fold (within 1 Å, or so), the MD calculations show that the architectural sub-Angstrom rigidity of the protein is not a key factor in making the channel selective for K+ over Na+. Rather, selectivity is controlled by the dynamical and electrostatic properties of the carbonyl ligands surrounding the ion (Noskov et al., 2004), a perspective that shares some of the basic ideas originally proposed by Eisenman (1961).
The huge variation in electrostatic profiles, over the ensemble of structures sampled by MD, shows that calculations based on static structures and continuum solvent can lead to vastly differing energy profiles, even for surprisingly small variations in atomic coordinates. In addition to these difficulties, the implementation of strategies for defining electrostatic continuum regions encounters severe unresolved problems, such as whether a high dielectric constant should be used where a single water molecule cannot fit.
Though it is possible to “engineer” simplified models to yield (seemingly) reasonable results, this involves making numerous arbitrary decisions about dielectric constants, atomic radii, charge distribution, diffusion coefficients, and channel structure. By arbitrary decision, we mean that the values of the microscopic parameters of the model are adjusted so as to “fit” the final outcome of the model (e.g., the channel conductance) to the experimental results, as opposed to being chosen on the basis of more fundamental chemical and physical knowledge. The microscopic insights into the mechanisms of ion permeation and selectivity are severely limited in such models. It becomes necessary to treat explicitly the atomic interactions and thermal motions and deformations upon ion binding, as is done in MD simulations. The use of a rigid protein model leads to extreme sensitivity to the precise (sub-Angstrom) values of the atomic coordinates, for example, but such spurious problems disappear when one considers an equilibrium thermal average corresponding to a PMF calculated by MD to describe ion permeation.
We wish to establish a connection between atomic thermal fluctuations of a protein and its flexibility in response to a perturbing force. Consider a solvated protein in a membrane with coordinates and velocities in the system denoted by R and V, respectively. According to Maxwell-Boltzmann statistics, the joint probability P(R,V) for the velocities configuration is given by
The above argument can be extended to any multidimensional case of perturbation affecting the set of coupled coordinates of the protein atoms,
where Xγ denotes the Cartesian coordinate (x, y, z) of protein atom γ. The average for the coordinate Xα is then
Eq. A4 also can be written in terms of a set of cross-correlation matrix elements
Comparison with Eq. 2 shows that the inverse of the matrix C acts as effective spring constants maintaining the protein structure.
We thank Drs. Peter C. Jordan, Zhe Lu, and Chris Miller for stimulating discussions and suggestions.
This work was supported by the Revson Foundation (T.W. Allen) and National Institutes of Health grants GM21342 (O.S. Andersen) and GM62342 (B. Roux).
Ted Begenisich served as guest editor.
T.W. Allen's present address is Department of Chemistry, University of California, Davis, Davis, CA 95616-5295.
Abbreviations used in this paper: gA, gramicidin A; MD, molecular dynamics; NMR, nuclear magnetic resonance; PMF, potential of mean force; RMS, root-mean-square.