INTRODUCTION

TOP ABSTRACT INTRODUCTION MODEL SYSTEM SIMULATIONS RESULTS AND DISCUSSION CONCLUSIONS REFERENCES

The recent determination of the structure of the potassium channel from Streptomyces lividans (Doyle et al., 1998) has inspired this molecular dynamics (MD) investigation, examining the ion selectivity, structure, and dynamics within this complex biological ion channel. The reason for the considerable interest in this channel is its role in neuronal signaling and its existence in all biological organisms. From the available structural information, the potassium channel appears to provide a unique mechanism for ion permeation, far different from the usual ion channels studied via MD simulation. The potassium channel is classed as a long pore channel where ions queue to pass the pore in a single-file fashion. This multiple ion mechanism, combined with the high throughput of 10 pA (Schrempf et al., 1995) and selectivity margin of 10⁴ for potassium over sodium, make this an interesting channel to study by computer simulation.

Molecular dynamics simulations exist for a variety of biological channels from the study of water properties in simple structureless cylinders (Sansom et al., 1996), and small bacterial pores such as gramicidin A (Roux and Karplus, 1994) to complex channel proteins like the sodium and nAChR channels (Singh et al., 1996; Sankararamakrishnan et al., 1996) and the OmpF porin (Tieleman and Berendsen, 1998). The intensive nature of this computational technique leads to an inevitable tradeoff between system size and simulation time scales. In this study we use a combination of experimentally determined protein structure, for the segment of the protein where microscopic properties are deemed vital to the operation of the channel, and simplified models for the less sensitive regions. In this way it is possible to simultaneously treat important microscopic details and generate the statistics required for the physical properties sought.

Within the short time frames permitted by MD simulation, it is unlikely that conduction can be observed (each event requiring 10 ns on average). However, MD trajectories may reveal the existence of multiple ion configurations that can be confirmed by experimental charge densities. These simulations also unveil the water structure (not seen in x-ray diffraction at 3.2 Å resolution), as well as channel water self-diffusion and rotational correlation, ion coordination, and total and free energy profiles from which ion selectivity can be inferred. Finally, we determine the ion diffusion throughout the potassium channel. By combining these findings with the longer time-scale Brownian dynamics simulations in the accompanying paper (Chung et al., 1999), we hope to further our understanding of both the microscopic properties and the conduction mechanisms underlying this important biological ion channel.

	MODEL SYSTEM

TOP ABSTRACT INTRODUCTION MODEL SYSTEM SIMULATIONS RESULTS AND DISCUSSION CONCLUSIONS REFERENCES

Our model system, shown in Fig. 1 A, is based on the protein structure detailed in Doyle et al. (1998). The model includes a narrow selectivity filter, formed by four peptide chains, augmented by a long hydrophobic channel wall consisting of a wide chamber and interior pore. Much of the protein forming the potassium channel has not been treated explicitly in these simulations. The focus of this study is the key narrow segment that is expected to determine ion selectivity, and is certain to have considerable influence on the channel conductance.

View larger version (24K): [in this window] [in a new window]   FIGURE 1   (A) Cross-section of the model channel. Shown are the hydrophobic pore and selectivity filter protein. The solid curves outline the effective radius of the Lennard-Jones 5-3 hydrophobic pore. Double arrows indicate effective pore and reservoir diameters (taking into account van der Waals radii). Reservoirs extend outside the range of this figure to ±35 Å where periodic boundaries are applied. Helix and mouth dipoles are represented by vector arrows. Solid lines of length 5 Å oriented perpendicular to the channel axis at z = 22 Å represent the harmonic blockers that prevent water from moving behind the protein. (B) Projections of two subunits forming the selectivity filter. Carbon, nitrogen, and oxygen atoms are shown in white, gray, and black, respectively, while hydrogen atoms are not shown. Initial (r, z) coordinates of the carbonyl oxygen centers on the G1, Y, G2, V, and T residues, and of the T hydroxyl oxygen are (5.25, 22.97), (3.20, 20.38), (2.22, 17.26), (2.34, 14.68), (2.76, 10.23), and (3.11, 9.21) Å, respectively. Constrained groups of atoms are indicated with arrows. Gray and black arrows refer to endpoint constraints and interactions with the pore helix, respectively.

The selectivity filter (Fig. 1 B) is reasonably flexible, molding itself according to the positions of water molecules and ions within the pore. It typically has a mean radius of ~1.5 Å and length 12 Å (extending from ~z = 10-22 Å) and exhibits fourfold symmetry about the channel axis (z). Each subunit of the tetramer of peptide chains consists of 158 KcsA signature sequence amino acids. In this study we implement only those five residues lining the selectivity filter, namely glycine (G₁ and G₂), tyrosine (Y), valine (V), and threonine (T). The polar C==O and C-O-H groups are directed toward the channel axis. Initial coordinates for this segment of the protein have been extracted from the Protein Data Base set of Doyle et al. (1998). The average radial positions of the carbonyl-oxygen atoms that line the pore in this coordinate set are as small as 2.2 Å near the center of the selectivity filter, corresponding to an effective pore radius of ~0.8 Å. Given that the Pauling radius of a potassium ion is 1.33 Å and that of an oxygen is ~1.4 Å, it is clear that significant bending of the subunits is required for conduction to occur. The low resolution of the diffraction experiments from which the protein coordinates have been derived means that atomic positions must remain relatively free to adjust and settle into optimal values that ensure functionality of the channel. We do this by using weak harmonic constraints on the selectivity filter that permit large deformations of the narrow central region of the selectivity filter while maintaining its overall structural integrity.

To prevent displacement of the selectivity filter subunits, harmonic constraints are applied to the endpoints, which are bonded to the remainder of the protein. G₁-C and T-N atoms are held with 40 kT/Ų force constants and neighbors G₁-O and T-C held with 20 kT/Ų throughout all simulations (1 kT = 4.14 × 10²¹ J = 0.595 kcal/mol for T = 300 K). These constraints mainly hold the endpoints of the filter near their observed positions, but otherwise have no effect on ion selectivity. The Y side chain is expected to experience van der Waals and possible hydrogen bonding interactions with tryptophan residues on the pore helix. We mimic this by applying constraints to the Y-OH oxygen and four ring carbon atoms. Polar N-H bonds, directed away from the pore axis, are expected to interact with side-chains on the pore helix (a combination of electrostatic interaction and hydrogen bonding). All other atoms, including the important carbonyl groups lining the pore, are free to move. We show in the Results section that ion selectivity is an intrinsic property of the filter structure independent of the applied harmonic constraints. This provides justification for our truncation of the channel protein to the peptide chains in the filter. We emphasize that our purpose in this article is not to measure ion selectivity accurately, but simply to demonstrate that selectivity arises from the experimental protein geometry.

The hydrophobic pore (Fig. 1 A) extends from z  20 to z  10 Å with a minimum radius of 3 Å at the intracellular entrance (z = 20 Å), to a maximum of 5 Å at z = 5 Å. The lack of polar groups and the increased dimensions of this region of the pore permit a one-dimensional potential model to be used. Our hydrophobic wall potential has been derived by Boltzmann averaging the potential due to a regular array of Lennard-Jones (LJ) 12-6 centers on finite length cylindrical channels and fitting an LJ potential. The resulting wall potential that is applied to water oxygen atoms and ions is an LJ 5-3 function with potential minimum of 1.1 kT at 3.45 Å (Allen et al., 1999). This function provides good fits to the structured hydrophobic wall over the range of pore dimensions experienced within the potassium channel. Reservoirs are also formed by an LJ 5-3 potential wall, have maximum diameter 16 Å, and are joined by a periodic boundary at z = ±35 Å. These regions are large enough to provide bulk-like ion environments and thus suitable reference points for energy profiles. A small harmonic blocker (see Fig. 1 A) applied at z = 22 Å and r > 5 Å with force constant 8.5 kT/Ų prohibits reservoir waters from entering the protein from behind the selectivity filter.

The pore helices have a net dipole field that is represented by a set of four positive charges at z = 22.35 and radial position r = 17.35 Å, and a set of four negative charges placed on the T-N atoms (the "connection points" of the selectivity filter and pore helices) ~8 Å from the center of the wide hydrophobic chamber. The magnitude of the charge at each end of an isolated -helix is 0.8 × 10¹⁹ C (Hol et al., 1978). During initial simulations we use a slightly reduced charge of 0.6 × 10¹⁹ C and subsequently investigate the effect of the magnitude of this charge on energy profiles. The pore helix macro-dipole field is not subject to a long-range cutoff. Due to the close proximity of the negative terminals and the atoms of the selectivity filter, pore helix fields are not applied to protein atoms. This will be justified later when we examine the influence of external fields on protein structure.

Acidic amino acids at the channel entrances result in the exclusion of anions (Doyle et al., 1998). The effect of mouth dipoles near the extracellular channel entrance on energy profiles is investigated. An additional charge of 0.6 × 10¹⁹ C is added to the G₁-N, with a matching positive charge at a distance of 5 Å in the radial direction. An equivalent set of dipoles with negative terminals at z = 20 Å and r = 4 Å are used during permeation studies.

Uniform electric fields are applied during permeation tests to drive ions across the channel. A voltage difference of ±500 mV is put across the planes z = ±50 Å resulting in a 5 mV/Å field, which is much larger than the field due to membrane potential. This strong electric field helps speed up ion permeation during the short simulation periods.

	SIMULATIONS

TOP ABSTRACT INTRODUCTION MODEL SYSTEM SIMULATIONS RESULTS AND DISCUSSION CONCLUSIONS REFERENCES

MD simulations are performed using the CHARMM v25b2 code (Brooks et al., 1983) with the CHARMM19 united atom force field. A time step of 1 fs is used with the Verlet MD algorithm (Verlet, 1967). Temperatures are maintained by velocity rescaling at 1-ps intervals, leading to temperatures of 300 ± 1 K for all simulations. Coordinates are saved each 0.1 ps for trajectory analysis. Periodic boundaries are enforced with a single image at each end of the channel.

Potential models

In deciding which water model is most suitable for this study we consider two major factors. The first is the ability of the model to reproduce bulk water and aqueous solution properties, and the second is the ability of the water model to enable solvation of the hydrophilic interior of the selectivity filter. We consider three popular models for this survey: TIP3P (Jorgensen et al., 1983), ST2 (Stillinger and Rahman, 1974) and SPC/E (Berendsen et al., 1987).

The TIP3P model (Jorgensen et al., 1983) is the least successful among the three models in reproducing the structural properties of bulk water (Daggett and Levitt, 1993). It has a relatively high diffusion constant (0.40 Å²/ps; Daggett and Levitt, 1993) compared to the experimental value (0.23 Å²/ps; Lide, 1994). The modified-flexible version of this model has a reduced level of self-diffusion (0.32 Å²/ps; Sansom et al., 1996), but this is still a poor estimate in comparison with other water models. In our MD simulations with the TIP3P and the modified-flexible TIP3P models, we find that they are unable to solvate the lining of the selectivity filter. After 0.5 ns of equilibration with either model, water molecules could not reach the center of the selectivity filter. Similarly, water molecules have difficulty approaching a potassium or a sodium ion located in this region. The inability of water molecules to enter the hydrophilic selectivity filter region and the resulting energy barriers for both ion species eliminates any possibility of using the TIP3P model in this study. Since TIP3P has been used successfully in the study of other hydrophilic pores such as gramicidin (Roux and Karplus, 1991), we attribute this shortcoming to the smaller dimensions of the selectivity filter (it is narrower and shorter than gramicidin). Modifications to the water-protein interaction would be required for the TIP3P model to be used in studies of the potassium channel.

Both the ST2 and SPC/E water models adequately solvate the protein pore lining. The ST2 model has the ability to accurately represent water in aqueous solutions across the range of static and dynamic ion and water properties. For example, a bulk water self-diffusion value of 0.285 Å²/ps has been found for ST2 water (Heinzinger, 1985), which is in reasonable agreement with experiment. Hydration numbers, radii, and ion-water geometry are well-reproduced with the ST2 model for a variety of ion species (Heinzinger, 1985), despite the apparent exaggeration of charge distribution with the five-site model. A drawback of this model is the increased trajectory storage and simulation times associated with the explicit treatment of lone pairs. However, CHARMM offers special fast routines for this water model, overcoming this problem to some degree.

The SPC/E model also provides good agreement with experiment for bulk self-diffusion with the value 0.25 Å²/ps (Berendsen et al., 1987). The model accurately reproduces bulk sodium ion hydration geometry and ion diffusion (see, for example, Guàrdia and Padró, 1990), as well as adequately representing water structure and dielectric response. Nevertheless, we prefer the ST2 model in our main simulations because of the availability of tested, specially designed ion-water interaction potentials for this model (Heinzinger, 1985), and the apparent advantages of using a water model implicit in our MD package. The standard rigid five-site ST2 water model has been used, with bond lengths, angles, charge distribution, and LJ 12-6 parameters taken from Stillinger and Rahman (1974). Ion-water LJ 12-6 potential functions are taken from Heinzinger (1985). Ion-protein oxygen, carbon, and nitrogen LJ 12-6 interactions are determined from ion-ion and protein-protein LJ 12-6 parameters by CHARMM with the use of standard combination rules. Ion-ion LJ 12-6 parameters are taken from Heinzinger (1985), while CHARMM19 protein-protein parameters are used. Residues in the peptide chains are given CHARMM19 bonding, atomic partial charges, and potential parameters. A switched force cutoff for electrostatic interactions and a switched potential cutoff for LJ 12-6 interactions are applied to groups at 12 Å.

In simulations with the SPC/E water model, we use the geometry, charge distribution, and LJ 12-6 parameters provided by Berendsen et al. (1987). Ion-water interactions of Heinzinger (1985) are found to be unsuitable for use with this water model. The parameters of Pettit and Rossky (1986), as tested with SPC/E by Guàrdia and Padró (1990), are used to describe ion-water interactions. Nonbonded interactions are atom-based for this water model using the same cutoff procedure as for the ST2 calculations.

Simulation protocol

The water in the channel and reservoirs is sampled from pre-equilibrated bulk coordinates by keeping only those waters that have oxygen atoms within the effective radius as shown in Fig. 1 A, and shifting the channel so as to obtain 1 g/cm³ density. Before the introduction of ions to the channel, equilibration of the water-filled system consists of 1000 steepest descent minimization steps, 10 ps of heating to 300 K, and a further 110 ps of dynamics. An ion is included by placing it at a particular point on the z axis and moving any waters closer than 2 Å into the reservoirs. The number of waters in the system remains unchanged throughout this suite of calculations, which is essential when using total system energies to determine energy profiles. An ion is held at a particular z-value with a harmonic constraint but is free to move in the xy plane to permit maximum solvation. A large value of 170 kT/Ų is used for the energy profile reported here. Energy profiles are seen to be independent of the strength of the z-restoring force. After a further 100 steepest descent minimization steps and 20 ps of equilibration and heating, between 0.5 ns and 1 ns of dynamics is carried out for each system from which mean energies are obtained. Ion energies are compared to a reference value at z = 29.0 Å in the extracellular reservoir, calculated from 1 ns of MD simulation. The choice of ion positions for energy profiles is such that the ion is placed at a point adjacent to a ring of oxygens lining the pore, or between neighboring rings, based on initial coordinates. This selection ensures that the most distinct energy states are sampled, and that intermediate ion positions have energies that are likely to be well-approximated by interpolation.

Bulk reference simulations with 528 water molecules in a periodic box of side 25.09 Å at density 1 g/cm³ and temperature 300 K are performed. After 110 ps of heating and equilibration, 100 ps of MD is used to determine bulk values of self-diffusion, rotational correlation, hydrogen-bonding, and water coordination. A single ion is added to the system, removing water molecules within 2 Å of the ion center, leaving 525 molecules. After 1000 steps of steepest descent minimization and 20 ps of heating and equilibration, 1 ns of trajectory data is used to produce bulk ion diffusion estimates for both potassium and sodium ions.

Analysis

The water-water and ion-water radial distribution functions (RDFs) are determined via the standard procedure (Allen and Tildesley, 1987). We refer to two different coordination numbers: a hydration number and a solvation number. The hydration number (number of nearest neighbors) is defined as the number of oxygen atoms (from the integrated RDF) falling within the first minimum of the bulk RDF. The "solvation" number of an ion is the total number of oxygen atoms (of any molecule) within this radius. A representative 100-ps sample from each simulation is used to calculate distribution functions.

Hydrogen bonding is measured by the time-averaged number of water-water pair interactions with total electrostatic plus LJ 12-6 energy exceeding a minimum attractive interaction. Cutoffs of 3.4, 5.1, or 6.8 kT have been used near hydrophobic and polar interfaces (Jönsson, 1981; Lee et al., 1984; Kjellander and Marelja, 1985). Here we choose an intermediate cutoff value of 5.1 kT applied in both bulk and channel water simulations. An estimate is also made for the number of hydrogen bond interactions between water molecules and protein carbonyl and hydroxyl groups. The neutral groups C==O and C-O-H are considered to behave as water molecules within the channel, the latter having negligible contribution.

Water self-diffusion and ion diffusion coefficients are calculated by averaging the mean square fluctuation of oxygen and ion position vectors over the range in time intervals 0.1  t  10 ps. Statistical uncertainty is minimized by using an overlapped data procedure (Rapaport, 1995). We consider the range within 0.1-10 ps after the initial shoulder, corresponding to inertial motion within the hydration "cage," and before the diffusion plateau, which may be a result of finite system size or simply due to statistical noise (fewer samples for large t). Before calculating diffusion coefficients, the center of mass of the system is subtracted to remove any net momentum associated with the periodic boundary. Simulations are broken into 10-12 samples and uncertainties are one standard error of means.

Water rotational correlation is determined by analyzing the decay of the dipole autocorrelation. The dipole autocorrelation function is defined as the time and system average of the cosine of the angle a water dipole at time t₀ + t makes with its dipole at time t₀. The sampling procedure is identical to that used in the diffusion calculation. A monoexponential decay is assumed and appears justified by analysis of the autocorrelation functions.

A more complete description of the channel's ability to discriminate between potassium and sodium ions comes from the calculation of free energy differences:

(1)

where U is the potential energy difference and S is the entropy difference. By comparing profiles of A and potential energy profiles, we may determine how significant entropy contributions are in reducing energy barriers that restrict ionic motion through the channel. Free energy perturbation (FEP) calculations are carried out with the PERT facility of CHARMM. The technique is based on a statistical mechanics relation between the free energy difference and the ensemble average of a function of the potential energy difference (see, for example, Mezei and Beveridge, 1986). The quantity we calculate with the FEP technique is A_KNa(z)  A_KNa(z₀). This is the free energy change associated with the creation/annihilation transformation of a potassium ion into a sodium ion with respect to the reservoir reference value at z₀. Potential parameters for the potassium ion are transformed into those for the sodium ion by means of a coupling parameter , also referred to as a molecule creation-annihilation coordinate, where  = 0 corresponds to the initial potassium ion state and  = 1 corresponds to the final sodium ion state. The potential energy function is given by the following linear combination:

(2)

It can be shown that integration of the ensemble average of the quantity dU/d from  = 0 to 1 reveals the total free energy difference (Mezei and Beveridge, 1986).

After 520 ps of equilibration, potential parameters of the potassium ion are transformed into those of the sodium ion in six windows with  values 0, 0.2, 0.4, 0.6, 0.8, and 1. Such a calculation is carried out for several key ion positions with the ion held by a harmonic constraint, as done with the total energy profiles. Similar calculations have been carried out in the past for the gramicidin channel (Roux and Karplus, 1991). Free energy differences between potassium and sodium ions are reported with respect to the reservoir value at z = 29.0 Å. Error estimates come from the measure of hysteresis found by determining the free energy change when the sodium ion is transformed into a potassium ion. It is found that free energy hysteresis remains low, even for very short simulation times, and that simulations of 120 and 600 ps generally give the same results to within ~1 kT. Results reported in this study involve 120 ps simulations for each ion, position, and transformation direction (KNa and NaK). Window sizes are 20 ps with the first 5 ps of each window used as equilibration and the last 15 ps for ensemble generation. Post-processing via the CHARMM Weighted Histogram Analysis Method (WHAM; a generalization of the histogram method of Ferrenberg and Swendsen, 1989) is used to find accurate estimates of free energy change after dynamics simulation.

	RESULTS AND DISCUSSION

TOP ABSTRACT INTRODUCTION MODEL SYSTEM SIMULATIONS RESULTS AND DISCUSSION CONCLUSIONS REFERENCES

We begin with a summary of protein structure and water self-diffusion results before examining ion coordination and energy profiles. Multiple ion configuration tests and ion diffusion estimates complete our calculations. All bulk reference values have been summarized in Table 1.

Protein structure

Fig. 1 B shows the initial protein structure in the absence of water and ions. When the channel is immersed in water, the narrow segment of the pore widens. The extent of deformation depends on the magnitude of the harmonic constraints applied to the filter. It is found that for N-H constraints of <10 kT/Ų the protein cannot maintain its shape. The selectivity filter tends to bulge and form a spherical-like cavity filled with many water molecules. When the constraints are increased above 10 kT/Ų, the filter maintains its structural integrity and the equilibrium geometry becomes more or less independent of the applied constraints. Initial tests involving 220-ps simulations showed that high selectivity can be achieved throughout the selectivity filter with the set of N-H group constraints K_G1 = 10, K_V = K_G2 = K_Y = 20 kT/Ų. Using this set as a starting point, we have investigated whether changes in the constraint values influence the selectivity property in an appreciable way. For this purpose we have carried out 520-ps simulations for each constraint value with an ion placed at the reservoir reference point z = 29.0 Å and the points immediately above and below the N-H group in question, adjacent to carbonyl oxygen rings. For example, for each value of K_G2 (the force constant applied to G₂-N and H atoms), the total system energy is recorded at z = 29.0, 17.3, and 14.7 Å (reservoir, adjacent to the G₂ oxygen ring, and adjacent to the V oxygen ring). The G₁ constraint K_G1 is studied at one point (z = 20.4 Å) because it has only one neighboring carbonyl ring position. These simulations are performed for both potassium and sodium ions so as to observe the effect of constraints on the ability of the channel to discriminate between these two ions.

Fig. 2 A shows the N-H constraint dependence of the total system energy for an ion at a particular position with respect to the energy of the same ion at the reservoir reference point. Solid curves (filled circles) represent the potassium ion relative energies while dotted curves (open circles) represent sodium relative energies. Inspection of the constraints K_V, K_G2, and K_Y shows that sodium energy difference is consistently higher than that of potassium at all four test positions. Thus a clear preference for potassium over sodium is observed independent of the constraint strengths applied to the V, G₂, and Y residues. In the case of K_G1, a similar discrimination occurs for 10 kT/Ų, but for larger values, potassium and sodium experience similar energy barriers. While energy barriers near the mouth of the channel display some dependence on the G₁ constraint strength, this has little impact on the overall selectivity of the filter. This survey has shown that the strength of harmonic constraints has little influence on channel function, and that our original choice of N-H constraints is suitable for simulations to follow.

View larger version (20K): [in this window] [in a new window]   FIGURE 2   Effect of harmonic constraints on energy. Each graph shows the dependence of the relative potassium (solid curve, filled circles) and sodium (dotted curve, open circles) energies with respect to the reservoir reference point (z = 29 Å), on the strength of harmonic constraints. (A) The uppermost panel displays the effect of the G1 force constant (KG1), applied to G1-N and H atoms, on the relative energies of potassium and sodium ions at z = 20.4 Å. The second row shows the effect of KY on the relative energies at z = 20.4 and 17.3 Å, the third row shows the effect of KG2 at positions z = 17.3 and 14.7 Å, and the fourth row shows the effect of KV on relative ion energies at z = 14.7 and 10.2 Å. (B) The dependence of relative potassium and sodium total system energies at z = 17.3 Å on the constraint strength applied to tyrosine side-chain atoms (ring carbons and hydroxyl oxygen). All error bars are one standard error of means.

The mean positions of carbonyl oxygens and the relative energies for potassium and sodium ions within the selectivity filter are found to change very little when the constraints on the Y side chain are varied. Fig. 2 B displays the dependence of the potassium and sodium relative energies on these constraints at z = 17.3 Å. It is seen that a large ion discrimination exists in this region independent of the strength of this constraint. Even in the absence of a harmonic constraint (K = 0), the large energy difference is maintained. The Y-OH oxygen and four ring carbon atoms are held with 10 kT/Ų restoring forces in this work to mimic interactions with tryptophan residues on the pore helix and to add some stability to the protein. The size of these harmonic forces is quite small, and similar constraints have been used in model nAChR studies (Smith and Sansom, 1997).

Fig. 3 A shows the selectivity filter with equilibrated water after 220 ps of dynamics with these constraints. In this two-dimensional projection as many as three water molecules are embedded in the protein between the subunits. The remaining water molecules form a single-file structure throughout the narrow segment of the protein. A typical equilibrated multiple ion system is shown in Fig. 3 B.

View larger version (23K): [in this window] [in a new window]   FIGURE 3   (A) Equilibrated structure for water-filled channel. (B) Equilibrated multiple ion system structure example. Water molecules are shown as small dark gray spheres, potassium ions as large black spheres. Diagrams have been generated with the use of RasMol v.2.6.

Channel water properties

Analysis of water density profiles within the selectivity filter reveals that water molecules tend to remain very close to the channel axis with average density approximately twice that of bulk water. Despite this increase, the number of nearest-neighbor waters 2.73 ± 0.05 is considerably lower than the bulk value of 4.70 ± 0.06 (Table 1) due to geometrical restraints. The number of hydrogen bonds per water molecule is 2.18 ± 0.02, also considerably lower than the bulk value of 2.73 ± 0.05. Considering the small number of nearest-neighbor water molecules, the water in the selectivity filter experiences a high percentage of hydrogen bonding (80% of all nearest-neighbor molecules are hydrogen-bonded). Because of the geometrical restraints, the average water molecule would be unable to remain hydrogen-bonded to the nearest-neighbor molecules if it were to rotate by any significant angle. The considerable energy cost incurred by the breakage of hydrogen bonds would prohibit the rotation of water molecules in the selectivity filter. Also contributing to this lack of rotational freedom are 0.22 (with standard deviation 0.20) water-carbonyl hydrogen bonds per water molecule, as well as long-range electrostatic interactions with the protein and pore helix field.

Analysis of dipole projection distributions reveals a clear preference for water dipoles to be directed in the z-direction, aligning with the pore helix dipole field. The z-component of the net dipole moment of this water column is 2.9 D with a standard deviation of 1.6 D (1 Debye = 3.34 × 10³⁰ Cm). This is similar to the dipole moment of a single ST2 water molecule (µ = 2.353 D). In comparison, the moment of the protein making up the selectivity filter is µ_z = 11.4 ± 1.8 D and that of the pore helix is µ_z = +21.7 D. The water dipoles respond to the conflicting fields of the protein and pore helix to give only a small net moment. When a 5 mV/Å uniform field is applied to the system, negligible polarization of the protein and water column is observed (the change of net dipole moments being less than one standard deviation). The inability of the protein to polarize in the presence of a strong external field justifies to some degree our not applying the pore helix field to the selectivity filter. Future MD studies involving the entire channel protein should reveal whether protein polarization has any effect on ion energetics.

The highly compact nature of water molecules within the selectivity filter and the significant effect of hydrogen bonds within this region are expected to result in a considerable loss of self-diffusion and rotational freedom. Fig. 4 shows the self-diffusion coefficient and inverse rotational correlation time constant across the length of the channel. Both the xy- and z-components of the self-diffusion appear to be significantly reduced in the selectivity filter, where, on average, the self-diffusion coefficient is 0.019 ± 0.005 Å²/ps: ~7% of the bulk value. This greatly diminished diffusion within the selectivity filter is a feature of narrow hydrophilic pores (Allen et al., 1999), and especially of single file systems (Lynden-Bell and Rasaiah, 1996). The study of molecular diffusion in zeolites (Hahn et al., 1996), for example, where diffusion is a single-file process, reveals two orders of magnitude reduction from bulk diffusion levels. Significantly diminished self-diffusion coefficients have also been observed within a variety of other model biological pores such as porin (Tieleman and Berendsen, 1998) and several poly-alanine and amphipathic bundles including nicotinic acetylcholine receptor-like channels (Smith and Sansom, 1997; Breed et al., 1996). The first-order inverse rotational correlation time is also very low in the selectivity filter highlighting the inability of water molecules to rotate.

View larger version (14K): [in this window] [in a new window]   FIGURE 4   Self-diffusion and rotational correlation: (A) Self-diffusion Dz (solid curve, filled circles) and Dxy (dotted curve, open circles). (B) First-order inverse rotational correlation time constant 1. A 5-Å grid has been used over the potassium channel with a single value for each reservoir. The dashed lines are the bulk reference values from Table 1. Error bars are one standard error of means.

The increased self-diffusion and rotational freedom near the intracellular channel mouth is associated with a decreased water density owing to the shape and dimensions of the repulsive wall in that region. In fact, the density is as low as 60% of bulk in the zone 20  z  15 Å where self-diffusion is considerably higher than bulk. The increased rotational freedom of water molecules could result in a high effective dielectric constant in this region: a property essential to the translocation of ions across it (Chung et al., 1999). The decrease in diffusion and rotation as the pore radius increases from 3 to 4 Å (z > 20 Å) has been observed in studies of cylindrical channels with the same wall potential (Allen et al., 1999). We therefore see a correspondence between diffusion within the potassium channel and diffusion within a cylinder of radius equal to that of a particular cross-section.

Ion coordination

Analysis of the bulk coordination properties in Table 1 shows that both potassium and sodium in ST2 water require approximately seven hydration waters to complete the first coordination shell, but that there is a difference of ~0.5 Å in the radius of the hydrated ions. Ion discrimination is therefore most likely to come from the proximity rather than the number of coordinating molecules. Within the reservoirs and throughout the hydrophobic pore, ions are observed with bulk-like hydration. The ions are, however, largely dehydrated in the selectivity filter region, in places only coordinated by two water molecules. The carbonyl groups lining the selectivity filter play a key role in allowing one ionic species to traverse the channel, while rejecting another.

Despite the flexibility of the protein, with carbonyl displacements as large as 0.7 Å and the ability of ions to drift by as much as 0.5 Å toward carbonyl oxygens, sodium is unable to maintain its first solvation shell throughout the selectivity filter. To illustrate this point, typical RDFs for both water- and protein-oxygen atoms are plotted in Fig. 5 for potassium (A) and sodium (B). Shown in the insets are ion-carbonyl oxygen geometries representative of each simulation. The large size of the potassium ion means that carbonyl oxygens sit well within the first hydration shell, and this is evidenced by the overlapping bulk hydration and protein peaks in the RDF. In contrast, we see that the first maximum in the protein RDF of sodium (with Pauling radius 0.95 Å) sits far from the first hydration maximum.

View larger version (18K): [in this window] [in a new window]   FIGURE 5   RDFs for potassium (A) and sodium (B) when an ion is held at z = 20.4 Å. The RDF including water oxygen atoms is plotted with a solid line whereas that for protein oxygens is plotted with a dotted line. The dashed curve is the bulk hydration RDF. The inset shows an ion-carbonyl oxygen geometry representative of the simulation.

Fig. 6 displays the positions of the first maxima in the hydration and protein RDFs within the selectivity filter compared to bulk. Throughout the entire selectivity filter the carbonyl oxygens are seen to contribute to the solvation of potassium, compensating for hydration losses, but are unable to coordinate a sodium ion to the same extent. In fact, on average over the selectivity filter (10.2  z  21.7 Å), there is a 0.27 ± 0.09 Å gap between first protein oxygen and bulk hydration maxima in the RDFs for sodium, whereas there is only a 0.11 ± 0.08 Å average space for potassium. Calculations repeated with the SPC/E water model are found to be almost identical to these ST2 model results. There is therefore little model-dependency in this figure and gaps between water and protein RDF maxima are a result of the protein structure. We conjecture that the radial positions of the first hydration maxima are smaller than the bulk values because of the increased water density in the selectivity filter region, and due to the increased space arising from the carbonyl oxygens being held slightly away from each ion.

View larger version (16K): [in this window] [in a new window]   FIGURE 6   Positions of the first maxima in the hydration and protein oxygen radial distribution functions of potassium (A) and sodium (B). Solid curves (filled circles) are hydration RDFs, dotted curves (open circles) are protein oxygen RDFs, and the dashed lines are bulk values.

Toward the ends of the selectivity filter, away from the G₂ and V residues, it is the width of the pore that results in poor sodium coordination. For example, a potassium ion held between Y oxygens at z = 20.4 Å leads to a mean Y oxygen radial position of 3.0 ± 0.3 Å such that the first maximum in the ion-carbonyl oxygen RDF is at 2.90 ± 0.04 Å. When a sodium ion is placed at the same position, carbonyl oxygens only approach the channel axis to within 2.8 ± 0.4 Å resulting in a carbonyl RDF maximum at 2.60 ± 0.08 Å. The Y carbonyl oxygens, originally at radial position 3.2 Å, are permitted to move in toward the channel axis by no more than 0.4 Å. This is sufficient to solvate the potassium ion, but not the sodium ion, which would require ~0.8 Å motion.

The central region of the filter (G₂ and V residues) would be expected to be narrow enough to solvate sodium. However, the carbonyls are not able to approach the ion because of the presence of hydration waters. Hydration waters, either directly or indirectly, force carbonyl oxygens attempting to solvate the ion away from the channel axis. Because the carbonyl oxygens are not free to diffuse around hydration waters, the ion cannot attain a bulk-like coordination. For example, when a sodium ion sits at z = 16 Å, the mean G₂ and V oxygen radial positions are 2.67 ± 0.22 and 2.51 ± 0.26 Å, respectively, resulting in a first RDF maximum at 2.74 ± 0.04 Å. Since the oxygens in the first hydration shell of sodium should reside near 2.34 ± 0.04 Å, it is clear that sodium should experience an energy barrier. Because the potassium ion has approximately the same radius as a water molecule, it does not encounter this difficulty. The compact nature of carbonyl groups near the center of the selectivity filter eliminates the possibility of the protein efficiently solvating an ion with a radius smaller than that of water. Thus, the widening of the narrow central region of the selectivity filter due to hydration waters appears to provide a second exclusion process for the sodium ion, which could not have been deduced from structural information alone.

The above selectivity mechanisms support the "close fit" hypothesis originally suggested by Bezanilla and Armstrong (1972). The idea is that the potassium ion fits nicely within the selectivity filter, but that the sodium ion is too small and would suffer a large dehydration barrier. Doyle et al. (1998) apply a similar interpretation to their molecular basis of the potassium channel. In particular they emphasize that the interaction of tyrosine side chains with tryptophan residues on the pore helix is responsible for holding the protein in place, resulting in this close fit. In this study we have seen little correlation between the strength of restoring force applied to the tyrosine side chains and the mean positions of carbonyl oxygens. Also, these constraints are not sufficient to maintain the integrity of the protein structure, and additional electrostatic interactions with polar N-H bonds on the backbones of the peptide chains are needed to prevent large motions of the carbonyl oxygens. Recent experimental inactivation studies of potassium channels (Kiss et al., 1999) have added weight to the close fit hypothesis. During inactivation the pore undergoes a constriction such that a shift from the nearly exclusive preference for potassium ions occurs, and the channel becomes less permeable to potassium ions and more permeable to the smaller sodium ions, eventually ceasing conduction altogether. These inactivation studies also reveal that the selectivity filter is not entirely rigid, as observed in our MD simulations.

To conclude this section on ion coordination, in Fig. 7 we show the first hydration (water oxygens) and solvation (water and protein oxygens) numbers for potassium (A) and sodium (B). Despite the almost total dehydration of the potassium ion, the solvation number remains above the bulk hydration value throughout the selectivity filter owing to the ability of carbonyls to emulate water molecules. In contrast, there are regions where sodium is unable to meet this quota. Calculations repeated with the SPC/E water model show an almost identical plot of coordination numbers for the potassium ion, but reduced solvation losses for the sodium ion due to the smaller bulk hydration value of ~5.6 (see Table 1). We therefore anticipate smaller relative energy barriers for the sodium ion with the SPC/E water model in comparison to the ST2 model.

View larger version (15K): [in this window] [in a new window]   FIGURE 7   First shell coordination (hydration and protein solvation) numbers for potassium (A) and sodium (B). Dotted curves (open circles) are hydration numbers, solid curves (filled circles) are total water plus protein solvation numbers, and the dashed lines are the bulk hydration numbers.

Energy profiles

We anticipate that the favorable interaction between potassium and carbonyl oxygens, together with the helix dipoles, will provide an attractive potential well for that ion. Sodium, however, is poorly solvated in regions where there is little stabilizing potential due to the helix dipoles and as a consequence will likely experience a significant energy barrier. The energy profile for the potassium ion over the selectivity filter is compared with that for the sodium ion in Fig. 8 A. It is clear that the channel will be impermeable to sodium, as it experiences a steep energy barrier that reaches a maximum of 30 ± 4 kT. The sodium ion experiences large energy barriers over the majority of the selectivity filter for (z > 12.4 Å) despite the influence of the strong pore helix dipole fields. The energy difference between sodium and potassium is between 6 and 36 kT in this region. The low energy of the potassium ion in the selectivity filter demonstrates the fine balance of energies associated with the lost hydration waters and the compensation by protein carbonyl groups. The smooth profile of sodium in comparison with potassium indicates the obvious involvement of the carbonyls of each protein residue in the coordination of the potassium ion as it traverses the selectivity filter.

View larger version (15K): [in this window] [in a new window]   FIGURE 8   (A) The total energy of the system with respect to the extracellular reservoir is plotted against ion position z for potassium (filled circles and dashed curve) and sodium (open circles and dotted curve). A profile for potassium corrected by placing mouth dipoles is drawn as a solid curve. All values are derived from simulations using the ST2 water model. Error bars are one standard error of means. (B) Potassium energy well depth (average of relative energies at points z = 5 and 7.5 Å) is plotted against pore helix charge. The fitted solid curve is the quadratic function E(kT) = 37 + 150q  295q2 with helix charge q in units 1019 C.

The global minimum of the potassium energy profile is 59 ± 4 kT at z = 7.5 Å, and that of sodium is 47 ± 3 kT at z = 5 Å. These deep energy minima are a result of the strong pore helix field. The depth of the energy well for potassium drops steadily as the pore helix charge is reduced, as shown in Fig. 8 B. Brownian dynamics simulations obtain physiological conductance with energy wells of ~30 kT (see Chung et al., 1999) which corresponds to a pore helix charge of ~0.45 × 10¹⁹ C in this study. The energy profile recalculated with this reduced pore helix strength shares all of the same localized features seen in Fig. 8 A, but with a reduced energy well in the hydrophobic chamber.

The small energy barrier for potassium at the channel entrance (20.4 < z < 21.7 Å) of 9 ± 3 kT can easily be annulled by charge residues or dipoles. The effect of mouth dipoles, with charge 0.6 × 10¹⁹ C, is to eliminate any potassium barrier as shown as a solid curve in Fig. 8 A. There appear to exist sharp fluctuations in the energy throughout the selectivity filter, arising from the inability of carbonyl groups to minimize the energy of an ion at every point on the channel axis. Since energy profiles are found to be fairly independent of the strength of the harmonic forces applied to the protein, the fluctuations in the energy must be a feature of the channel geometry and not an artifact of the constraints. Furthermore, it shall be seen in the following section that these local energy minima and maxima do not significantly alter ion positions or hinder the passage of ions through the pore.

An improved measure of channel selectivity can be obtained from the free energy difference between potassium and sodium ions. Since CHARMM does not permit free energy calculations with the ST2 water model, we have used the SPC/E water model for this purpose. This also provides an opportunity to investigate the water model-dependence of our findings. Before examining free energies, we show in Fig. 9 A the total energy profiles for potassium and sodium ions obtained from simulations with the SPC/E water model. The profiles are qualitatively similar to those found with the ST2 model (Fig. 8 A) in that the sodium ion experiences significant energy barriers and that there is a clear preference for the larger potassium ion throughout (note that neither of the profiles in Fig. 9 A includes the effect of mouth dipoles). As anticipated, the discrimination between these two ions is not as pronounced with the SPC/E water model; nevertheless, the function of the channel has not changed.

View larger version (16K): [in this window] [in a new window]   FIGURE 9   Energy profiles with the SPC/E water model. (A) The total energy profiles, with respect to the reservoir reference point (z = 29 Å) for potassium (filled circles and solid curve) and sodium (open circles and dotted curve). (B) Comparison of free and total system energy difference profiles. The difference in the free energy change, corresponding to the transformation KNa, between each point in this profile and the reservoir reference point (z = 29 Å) is drawn as a solid curve (filled circles). The total system energy differences relative to the reservoir are drawn as a dotted curve (open circles). Error bars are one standard error of means.

The difference between the sodium and potassium total energy profiles with respect to the difference at z = 29.0 Å (the reservoir reference point) is plotted as a dotted curve (open circles) in Fig. 9 B. The solid curve (filled circles) in Fig. 9 B represents the free energy difference between sodium and potassium ions with respect to the free energy difference at z = 29.0 Å. This figure reveals that entropy contributions alter sodium barriers with respect to potassium only very slightly throughout the majority of the selectivity filter, with only the narrow central region between the V and G₂ residues (z = 16-17.3 Å) exhibiting sizable changes. The free energy perturbation analysis has therefore revealed the same ion discrimination as the total energy analysis. The large barriers of as much as 18 kT for the sodium ion with respect to the potassium ion near the intracellular and extracellular entrances of the selectivity filter are easily sufficient to exclude sodium ions from this channel.

Multiple ion configurations

By placing unconstrained potassium ions at various positions within the channel we can observe the ability of the channel to hold one, two, or three ions at a time. Experimentally, three regions of high ion concentration are observed: 1) within the wide chamber (z  5 Å), 2) a broad region at the interior end of the selectivity filter (z  12 Å), and 3) a region near the exterior end of the filter (z  20 Å). In the absence of Coulomb repulsion of ions, the energetically favored positions for a single ion in the channel are determined by the energy profile (Fig. 8 A). It is fairly clear from this profile that the channel is likely to hold more than one potassium ion because of the depth of the energy wells. These energy wells, in particular the well within the wide chamber, appear too deep to permit a single ion to traverse the channel. We refer to the potential profiles in the presence of one or two ions in the channel reported by Chung et al. (1999) as a guide to how this channel is likely to support multiple ion configurations.

Two unconstrained ions are placed within the channel at various likely positions in the absence of an external field and in the presence of extracellular and intracellular mouth dipoles. After 200 ps of dynamics, the ions in the channel form one of two stable configurations. The first configuration is one where two ions placed in the selectivity filter remain in this region with average positions 13.1 and 21.7 Å. If the two ions are placed too close together this configuration is unstable. The second possible configuration has one ion in the chamber (at  = 2.7 Å) and one in the selectivity filter (near z = 11.9, 16.0, or 19.2 Å). A combination of these two configurations would explain the experimentally observed regions of high charge density.

Next, three ions are placed within the channel, one in the chamber at 5 Å and two at suspected positions within the selectivity filter. When the two filter ions are initially separated by >6 Å, this configuration is maintained with average positions after 200 ps of dynamics being 6.8, 11.5, and 21.6 Å. When initial positions differ from this optimal configuration, one ion is either expelled from the channel or moves down into the already occupied chamber. Applied electric fields also reduce the propensity for a three-ion system with the ion at position z = 11.5 Å falling into the chamber region without provocation. We propose that the three-ion configuration is merely an intermediate state occurring during conduction with the channel occupied by just two ions the majority of the time. This is also confirmed by Brownian dynamics studies (Chung et al., 1999).

With an average of 10 ns per conduction event, the ions would have to be perfectly primed for conduction if an ion were to pass through the channel within the short time frame considered here. Although the potassium channel shows greater conductivity in the intracellular-extracellular direction, it is known to conduct ions in both directions. To speed up permeation we consider the extracellular-intracellular direction, which allows us to utilize the repulsive potential of additional ions outside the channel (one in each reservoir, z = ± 25 Å). A uniform field of 5 mV/Å is applied in the z-direction and 400 ps of trajectory data are generated with a 2-fs time step. Fig. 10 shows a typical trajectory for an ion passing almost the entire length of the selectivity filter. The ion comes to rest at numerous positions along its journey corresponding to energy minima in the energy profile of Fig. 8 A. We also observe an ion moving from z = 23.0 (reservoir) to 18.8 Å (within the channel) at approximately the same time (to within 100 ps). Thus, the motion of one ion triggers the motion of another due to a change in the ion-ion repulsive force. This appears to be the onset of a conduction event. Also noteworthy is that when the constraints on the protein atoms are altered, no noticeable change in the ion trajectories is seen. This suggests that it is ion-ion repulsion that governs ion conduction, and not short-range interactions with the protein.

View larger version (20K): [in this window] [in a new window]   FIGURE 10   Ion trajectories. An external field of 5 mV/Å2 drives ions from extracellular to intracellular. Reservoir ions are not shown. Initial z-coordinates are 18.8, 21.7, and 23.0 Å.

Ion diffusion

The bulk potassium diffusion result (Table 1) is only ~53% of the experimental value. A similar finding for the sodium ion is also observed. Our large box with dimensions greater than twice the long-range cutoff of 12 Å should lead to accurate bulk diffusion estimates. Therefore the likely cause for the lower diffusion rates is to be sought in the potential models used for ions and water molecules. The large hydration numbers seen with the ST2 model (also observed in Heinzinger, 1985) will certainly lead to a lower diffusion. It could also arise from the treatment of long-range forces (Del Buono et al., 1994; Perera et al., 1995). Nevertheless, comparisons between the potassium channel and bulk simulations with the same truncation scheme should remain valid.

By combining all trajectories from multiple ion and permeation tests (with and without applied electric fields), it is possible to calculate ion diffusion constants from 19.5 ns of dynamics simulation. Little difference in ion diffusion between the different sets of simulations is observed. This can be attributed to the small drift velocities associated with the driving force over the 10-ps time intervals considered in the mean square displacement analysis. Fig. 11 displays the xy- and z-components of the potassium ion diffusion across the channel. The low level of potassium diffusion in the reservoirs is a result of one or more ions being trapped near the mouth dipole potential wells at the channel entrances, and due to the finite dimensions of the reservoirs. The z-component of the ion diffusion remains close to bulk ion diffusion (>50%) throughout the hydrophobic pore, as predicted by the analysis of cylindrical channels (Allen et al., 1999). Ions appear to be able to diffuse within zones of low density between cylindrical water shells and require little translation or rotation of water molecules. High levels of ion diffusion ( of bulk) have also been observed for sodium ions in the wide nAChR channel (Smith and Sansom, 1998). The xy-components are, however, drastically reduced from bulk levels. According to cylindrical channel studies, allowing for the difference in sodium and potassium ion sizes, very little xy-ion diffusion is expected for pore radius <4.5 Å. A slight rise is seen in the wide hydrophobic chamber where the pore radius reaches 5 Å, but the shape of the wall and the influence of the selectivity filter protein has also greatly suppressed xy-diffusion there.

View larger version (17K): [in this window] [in a new window]   FIGURE 11   Potassium ion diffusion. Dz (solid curve) and Dxy (dotted curve) are compared to the bulk value (dashed line). Error bars are one standard error of means.

Within the selectivity filter the z-component of the ion diffusion drops considerably. At the intracellular side of the filter region, the ion diffusion in the z-direction remains close to that of bulk, but drops to just 8% of the bulk level near the extracellular entrance. The average value over the selectivity filter region (considering 10  z  21 Å) is 0.039 ± 0.015 Å²/ps, ~ of the bulk diffusion. As one would expect, ion mobility is strongly correlated to the water self-diffusion within what is basically a single-file system. Naively one would anticipate that the selectivity filter would therefore be the bottleneck to conduction within the potassium channel. However, the Coulomb repulsion due to other ions can lead to large drift velocities (of the order 10 m/s), as observed during permeation tests (Fig. 10). This suggests that an ion possessing enough energy may overcome the barrier to diffusion created by the strong water-water and water-protein interactions. Furthermore, Brownian dynamics simulations with varying diffusion constants have shown that high conductances can be achieved with low ion diffusion within the selectivity filter (Chung et al., 1999).

	CONCLUSIONS

TOP ABSTRACT INTRODUCTION MODEL SYSTEM SIMULATIONS RESULTS AND DISCUSSION CONCLUSIONS REFERENCES

We have presented a molecular dynamics study of the potassium channel based on recent x-ray diffraction findings (Doyle et al., 1998). By using a simplified model where only 5 of the 158 amino acids in the KcsA signature sequence have been treated explicitly, we have observed some important features of the potassium channel. Harmonic constraints are used to emulate the interactions of the selectivity filter with the remainder of the peptide chain. A one-dimensional potential representation is used for the wide hydrophobic segment of the channel, which is not expected to influence results based on the success of this model in studies of cylindrical channels (Allen et al., 1999). Of the three water models considered in this study, TIP3P is found to be unable to solvate the selectivity filter and, therefore, excluded from further analysis. ST2 and SPC/E models, however, are successful in this respect, and further give similarly adequate descriptions of water self-diffusion and ion-water geometry (hydration radii and water orientations).

Using a minimal set of constraints to represent the binding of the selectivity filter to the remaining protein, a reasonably flexible structure has been obtained, in agreement with recent experimental evidence (Kiss et al., 1999). Despite this flexibility, a high level of ion discrimination has been observed. The potassium ion is found to remain completely solvated throughout the entire channel. Carbonyl oxygens directed toward the channel axis in the selectivity filter effectively compensate for the dehydration of the ion in that region. In contrast, the sodium ion cannot remain completely solvated in regions where the stabilizing interactions with the protein are weak. This leads to a significant energy barrier as high as 30 kT (and free energy barrier with respect to the potassium ion as large as 18 kT), explaining the rejection of sodium ions from this channel. There appear to exist two processes leading to large sodium ion energy barriers. The first occurs at the wider ends of the selectivity filter where "molecular springs" hold the carbonyl oxygens away from the small sodium ion. The second process occurs near the narrow central region where hydration waters prevent carbonyl oxygen atoms from approaching the sodium ion. The first process leads to the large free energy barriers for sodium, reinforcing the "close fit" hypothesis of Bezanilla and Armstrong (1972).

The low resolution of the experimental data from which protein positions have been determined means that protein atoms must be allowed some freedom to settle into optimum positions. Given the weakness of the applied constraints and the resulting flexibility of the protein, optimum atomic positions are achieved during dynamics simulation. The strengths of constraints on the protein are found to have little bearing on ion discrimination. This suggests that once bonded to the neighboring protein residues through the endpoint constraints, the selectivity filter will offer a high level of discrimination merely as a result of its equilibrium conformation. Comparisons with molecular dynamics energy profiles and semi-empirical calculations of potential profiles involving the entire experimentally determined protein, are currently being performed and should provide interesting comparisons with the results obtained with the simple model presented in this study.

Multiple potassium ion configurations observed during dynamics suggest that the channel can easily accommodate two ions, but that the three-ion system is relatively unstable. With the application of an electric field and increased ion concentrations, ions have been observed moving throughout the selectivity filter and hydrophobic pore. Localized features of the energy profile do not appear to control permeation.

The self-diffusion and rotational correlation of water throughout the hydrophobic pore are generally of the same order as bulk. However, due to the single-file nature of diffusion within the selectivity filter, and the extent of water-water and water-protein interactions, self-diffusion is found to be reduced by an order of magnitude from bulk levels and long rotational correlation times are observed. The ion mobility is reduced to just of bulk in this region. Despite this low diffusion within the selectivity filter, high drift velocities due to Coulomb repulsion of ions are seen. Thus, channel conductance appears to be driven not by ion mobility, but rather by the strength of electrostatic interactions with other ions and electric dipoles within the protein.