This Article Abstract Full Text (PDF) All Versions of this Article: biophysj.105.077073v1 90/10/3447    most recent Alert me when this article is cited Alert me if a correction is posted Services Similar articles in this journal Similar articles in PubMed Alert me to new issues of the journal Download to citation manager Citing Articles Citing Articles via HighWire Citing Articles via Google Scholar Google Scholar Articles by Allen, T. W. Articles by Roux, B. Search for Related Content PubMed PubMed Citation Articles by Allen, T. W. Articles by Roux, B.

Ion Permeation through a Narrow Channel: Using Gramicidin to Ascertain All-Atom Molecular Dynamics Potential of Mean Force Methodology and Biomolecular Force Fields

Toby W. Allen ^*, Olaf S. Andersen  and Benoit Roux 

^* Department of Chemistry, University of California at Davis, Davis, California; and Department of Physiology and Biophysics, Weill Medical College of Cornell University, New York, New York

Correspondence: Address reprint requests to Toby W. Allen, Dept. of Chemistry, University of California, Davis, One Shields Ave., Davis, CA 95616. Tel.: 530-754-5968, Fax: 530-752-8995; E-mail: twallen{at}ucdavis.edu.

	ABSTRACT

TOP ABSTRACT INTRODUCTION METHODS AND RESULTS DISCUSSION AND CONCLUSION ACKNOWLEDGEMENTS REFERENCES

 
We investigate methods for extracting the potential of mean force (PMF) governing ion permeation from molecular dynamics simulations (MD) using gramicidin A as a prototypical narrow ion channel. It is possible to obtain well-converged meaningful PMFs using all-atom MD, which predict experimental observables within order-of-magnitude agreement with experimental results. This was possible by careful attention to issues of statistical convergence of the PMF, finite size effects, and lipid hydrocarbon chain polarizability. When comparing the modern all-atom force fields of CHARMM27 and AMBER94, we found that a fairly consistent picture emerges, and that both AMBER94 and CHARMM27 predict observables that are in semiquantitative agreement with both the experimental conductance and dissociation coefficient. Even small changes in the force field, however, result in significant changes in permeation energetics. Furthermore, the full two-dimensional free-energy surface describing permeation reveals the location and magnitude of the central barrier and the location of two binding sites for K⁺ ion permeation near the channel entrance—i.e., an inner site on-axis and an outer site off-axis. We conclude that the MD-PMF approach is a powerful tool for understanding and predicting the function of narrow ion channels in a manner that is consistent with the atomic and thermally fluctuating nature of proteins.

	INTRODUCTION

TOP ABSTRACT INTRODUCTION METHODS AND RESULTS DISCUSSION AND CONCLUSION ACKNOWLEDGEMENTS REFERENCES

 
As computational methods increasingly are used to interpret or predict biomolecular function (1), it becomes important to critically evaluate their suitability. In principle, the molecular dynamics (MD) potential of mean force (PMF) approach offers the best route from computer simulation to experiment (2) in a manner that is consistent with the atomic and thermally fluctuating nature of proteins (3,4). In practice, computational studies of ion permeation face significant challenges due to the widely varying timescales of protein and membrane thermal fluctuations that become relevant when constructing permeation models. These fluctuations range from rapid bond, angle, and torsion fluctuations to side-chain isomerizations and large-scale protein and lipid conformational changes (5). Because these fluctuations underlie all protein function, they need to be incorporated in computational models that aim to interpret biological structure-function relationships (6). Moreover, because these thermal protein fluctuations are associated with large variations in ion energetics (3,7), it becomes essential to explicitly incorporate their effects in order to obtain appropriate equilibrium averages for the parameters of interest.

A direct connection between structure and function cannot easily be obtained via MD simulation, however, because ionic fluxes correspond to transit times of 10–100 ns—meaning that it becomes difficult to establish contact with experimental results (single-channel conductances and ion binding constants). To circumvent this difficulty, macroscopic (e.g., (8–10) or semimicroscopic (e.g., (7,11–13)) physical models, which treat some or all of the system as uniform dielectric media, have been invoked. In this study, we instead chose to keep the fully microscopic treatment and demonstrate that permeation can be accurately described via an equilibrium free energy surface that incorporates all of the thermal fluctuations of the ions, water, protein, and phospholipids and which is free of parameter fitting. For this purpose, we employ MD to sample a statistical ensemble of configurations for a fully explicit, atomistic system which, with ion mobility calculations, is fed into a phenomenological conduction model consistent with the PMF calculation (14,15). We show that this approach can be used to obtain a rigorously defined, well-converged, and consistent free energy surface and demonstrate that present-day MD force fields—in conjunction with present-day computational methods—are, perhaps surprisingly, accurate in describing ion permeation through a narrow pore.

To test the approach, and its ability to predict experimental observables, we chose the gramicidin A (gA) channel as our test case because this channel, with its single file pore, poses an extreme challenge for computational investigations of molecular function. The gA channel structure is known, being a single-stranded, right-handed ß^6.3-helical dimer (16), which has been thoroughly characterized structurally (17–20) and functionally (21–25). It is also small enough to allow good sampling with rather modest computational resources (26). The gA channel thus is an excellent system for testing how well MD-PMF simulations can be used to predict complex molecular functions. Since the first MD simulations on this molecule in 1984 (27), several studies have improved our understanding of the microscopic mechanisms of ion permeation (for review, see (28)). The aim of this study was to critically examine, and hopefully validate, this approach by directly finding contact with experimental measurements.

A perennial problem in previous studies has been that MD simulations predicted free energies corresponding to rates of ion movement several orders-of-magnitude smaller than the measured rates (27,29,30). Perhaps the most notable example was the pioneering computations by Mackay et al. (27), in 1984, which revealed a barrier of nearly 40 kcal/mol opposing translocation of a Cs⁺ along the axis of the gA channel (these results were reported in (31)). In more recent studies (29,30), the MD-PMF profiles for K⁺ across the gA channel were constructed using no more than 80 ps simulations per umbrella sampling window (30) (see also (32)). These calculations also predicted barriers for ion permeation that were several kcal/mol too high to be compatible with experiment (possibly as much as 7 kcal/mol, based on the findings of this study), meaning that the predicted rates of ion movement were approximately five orders-of-magnitude too low. This could suggest that the atomic force fields used in MD simulations are not adequately calibrated and therefore unable to sufficiently stabilize ions within the narrow pore, as compared to bulk water. Methodological limitations, such as the construction of starting configurations, equilibration, and simulation times also might account for the poor prediction of experimental observables. Indeed, a recent study in which the sampling of the equilibrium distribution of ions was greatly extended (15), and in which destabilizing effects of periodicity and hydrocarbon polarizability were accounted for, led to a well-converged PMF that allowed for semiquantitative prediction of experimental observables. Therefore, as these computational studies have become more sophisticated (with improved methodologies and computational sampling), they provide increasing confidence that MD-PMF calculations could become a powerful tool for understanding (and eventually predicting) ion permeation (and, by implication, other molecular functions).

The difficulties encountered in MD simulations of ion permeation can be traced largely to the observation that the measured rates of ion permeation for gA channels are very high (at low permeant ion concentrations comparable to predictions based on a simple waterfilled pore immersed in bulk water (33)). This means that there cannot be a major energy barrier for ion movement, such that the energy profile for permeating ions result from almost complete cancellation of two very large opposing contributions: ion dehydration and protein/pore water solvation. MD simulations going back to 1984 (27) identified a possible major role of the single-file water to overcome the large dehydration barrier, and almost complete cancellation of the barrier by a combination of the single-file water and protein was demonstrated when the overall PMF was decomposed into water, protein, and membrane/bulk water contributions (15). This need to accurately represent ion solvation in both extremes of bulk water and almost complete dehydration in a narrow pore poses significant challenges to MD force fields.

The MD-PMF strategy adopted in this study is based on the assumption that the long-time behavior of an ion permeation event is dominated by some rate-determining step(s) and that its dynamical evolution can be described as a progress along some reaction coordinate (2). The true reaction coordinate for a given system is not known a priori, and the MD-PMF strategy consists in choosing some suitable order parameter (e.g., the position of the translocating ion along the channel axis) as a mathematical surrogate for the true reaction coordinate. The underlying assumption is that all other variables fluctuate rapidly, such that they can be integrated over in order to obtain a PMF for the chosen order parameter(s). Different choices of order parameters may provide adequate descriptions of the rate-determining step(s)—and usually yield slightly different PMFs—but the absolute transition rate is insensitive to such choices as long as dissipative factors are considered properly (34,35). In our case, relatively slowly varying degrees of freedom have been associated with the orientation/reorientation of the single file water column (15), which could be important for describing the microscopic dynamical mechanisms of permeation, in particular the kinetics of ion entry/exit. For now, however, we describe the permeation mechanism in terms of the ionic spatial coordinates alone, as a useful simplification.

Within this framework, we pursued strategies for identifying the most vulnerable aspects of MD-PMF simulations with the goal of improving their ability to predict experimental observables. To this end, we examine the dependence of ion conduction observables on the choice of MD force field and explore the sensitivity of the results to small changes in parameters. We conclude: first, that careful attention to the physical system at hand allows for significantly improved predictions of experimental observables; second, that simple changes in a given force field are unlikely to provide significant improvements; and third, that while semiquantitative agreement with experiment can be achieved with a modern all-atom fixed-charge force field, ultimately an electronically polarizable simulation will be required for further improvement.

	METHODS AND RESULTS

TOP ABSTRACT INTRODUCTION METHODS AND RESULTS DISCUSSION AND CONCLUSION ACKNOWLEDGEMENTS REFERENCES

 
Ion channel-membrane simulations
Simulations were done with the program CHARMM (36) using the PARAM27 (37) force field (referred to as CHARMM27 from here on), with standard protein (37) and TIP3P water (38) with ion parameters from Beglov and Roux (39). Simulations were also done with the all-atom AMBER PARAM94 (40), using TIP3P water, with ion Lennard-Jones parameters from Aqvist (41), and also with the united-atom GROMOS87 (42) force field, using the SPC water model (43) and ion parameters from Straatsma and Berendsen (44). We employ the AMBER94 force field previously imported into the CHARMM program (45) and imported GROMOS87 into CHARMM for comparison. In each case, the standard CHARMM27 lipid parameters were used (46) so as to isolate the effects of the protein when comparing PMFs. The use of particle-mesh Ewald (47), SHAKE (48), and constant pressure and temperature algorithms (49), have been described previously (15,20).

Systems that consist of a gA helical dimer (Protein Data Bank (PDB) No. 1JNO (18)) embedded in a DMPC bilayer (Fig. 1), were created using extensions of previous membrane-building techniques (50). The choice of starting gA structure is based on evidence (20) that dynamical trajectories starting with the PDB:1JNO structure reproduce experimental solid-state NMR measurements (51) better than the solid-state NMR PDB:1MAG structure (19). Membrane patches of approximately one and three shells of lipid molecules around the gA protein were used in the simulations. These patches consisted of 20 and 96 lipid molecules, and 1080 and 3996 water molecules, respectively. For the smaller one-shell system, hexagonal periodic boundaries of xy-translation length 32.1 Å, as determined from the area of the protein and lipids, and average height 74 Å, were imposed on the protein-membrane system. For the larger three-shell system, the hexagon xy-translation length was 61.9 Å, with average height 75 Å. Pressure coupling was employed in the z-direction (parallel to the membrane normal); the x,y dimensions of the hexagonal boundaries remained fixed during simulations. A 1 M KCl ionic solution was to used ensure good sampling of the ionic bath; this corresponds to 19 K⁺ and Cl^– pairs in the smaller system and 74 pairs in the larger system. Fig. 1 A shows the gA ion channel embedded in one shell of lipids. Fig. 1, B and C, shows the small and large systems, respectively, from the top with periodic images.

View larger version (145K): [in this window] [in a new window]   FIGURE 1  Gramicidin A in the bilayer: (A) one-shell system: gA dimer (yellow); DMPC bilayer atoms C (gray), O (red), N (blue), and P (green); K+ (green spheres) and Cl– (gray spheres); water O (red) and H (white). Within the channel, seven single-file water molecules are drawn as spheres adjacent to a single K+ ion at the channel entrance. The chosen MD frame has a channel axis with tilt angle 9° relative to the membrane normal vector z. Some lipid molecules and electrolyte from neighboring images are visible. (B) One-shell system (CPK color with green lipid C atoms) with hexagonal periodic images (CPK color with gray lipid C atoms) viewed along the membrane normal. Water molecules and ions have been removed for clarity. (C) Three-shell (large) system with hexagonal periodic images.

View larger version (47K): [in this window] [in a new window]   FIGURE 2  Observed rotameric states of the Trp-9 residues in the gA channel during 47 ns of unbiased simulation (adapted from (20)). The solid box highlights the correct rotameric state and indicates the placement of a two-dimensional flat-bottomed restraint to maintain this rotamer during PMF calculations.

View larger version (114K): [in this window] [in a new window]   FIGURE 3  Density of lipid heavy atoms around the gA dimer. Density is plotted as a two-dimensional histogram in axial and radial directions with respect to the protein. The one-shell histogram is an average over 10 ns of simulation whereas the three-shell system is an average over only the first 2 ns of a 10-ns simulation to highlight the presence of a lipid density over the channel entrance.

Next, while the water inside the channel maintains a single-file column without any manipulation, occasionally small gaps in that column may occur; especially when an ion is in close proximity to the channel entrance. To test if a system configuration was suitable for a starting point for a window, the maximum space between waters inside the channel was computed. This was done by checking for any gap greater than 1.5 Å in the range –10.5 z 10.5 Å. No constraints were used to maintain the pore water structure during simulations. (These breaks in the water column may become important when predicting, or interpreting, the diffusion coefficient of the ion-water column within the pore.)

Constraints were applied to ensure the membrane and protein are kept near the center of the periodic box. A weak harmonic planar constraint, of force constant 5 kcal/mol/Ų was applied to the z-coordinate of the center of mass of the lipid bilayer to prevent drifting. A weak center-of-mass constraint was also applied on the xy position of the center of mass of the channel by application of a cylindrical harmonic constraint of force constant 5 kcal/mol/Ų. These constraints have no effect on the z position of the channel relative to the membrane, nor the tilting of the channel. They only act to center the membrane and channel independently and have no impact on results.

Reaction coordinate for ion permeation
Following the general statistical mechanical equilibrium theory formulated in Roux (53) and reviewed in Roux et al. (2), the system can be separated into pore and bulk regions, which allows the definition of the free energy surface generated for pore occupation by n ions with coordinates r_i (i = 1, n). One can form a hierarchy of n-ion PMFs for different occupancy states of the pore (53). For low-to-moderate ionic concentration, the gA channel should be occupied by just one cation (25,54), and a one-ion PMF, will reveal much about the function of this ion channel in this regime. This one-ion PMF can be written in terms of a configurational integral (see Eq. 13 in (53)). (As noted above, to achieve satisfactory sampling, the simulations were done using a 1 M KCl solution where the channel may be occupied also by two ions, as noted below; this will not affect any of our conclusions, except that they pertain only to the one-ion case.)

To calculate a meaningful one-ion PMF, we must choose a pore region which is almost exclusively occupied by a single ion. Table 1 shows distributions of ion occupancies on either side of the channel (n_left, n_right) as a function of the size of an exclusion sphere, centered on the origin, based on analysis of 10.9 ns of unbiased simulation. A suitable choice for this radius appears to be 14 Å, above which it becomes possible to see two cations on one side of the channel within the sphere. Furthermore, with this radius, an anion is found inside the sphere during only 1% of simulation, which is evidence of valence selectivity of the gA channel (this question will be examined further in a separate study). Still, 19% of the time, one cation may be bound at both entrances of the channel. To calculate a one-ion PMF, starting configurations were chosen such that only the one ion was in the 14 Å sphere, and during biased simulation other ions (cations and anions) were excluded with a repulsive flat-bottom spherical harmonic restraint with force constant 5 kcal/mol, applied to other ions only when they enter this exclusion sphere. Thus, our simulation methodology is designed to compute the one-ion PMF from MD simulations that rule out multiple pore occupancy (as would occur with 1.0 M K⁺ in the aqueous solution (54)).

Because the ion is confined within a narrow region in the xy-plane within the channel, one may assume that the equilibrium distribution of lateral displacements is obtained quickly and thus may be integrated away, leading to the one-dimensional PMF W(z), or free-energy profile (2,53). However, as reported previously (15), the ion becomes unbounded in the xy-plane at a distance of 14–15 Å from the channel center. Thus, this one-dimensional profile has limited significance outside the channel because integration over all xy extents will cause the PMF to tend toward –. To obtain an unambiguous free energy profile one must restrict the lateral displacement of the ion. In the current computations, a flat-bottom cylindrical constraint with radius 8 Å (relative to the center of mass of the dimer) of force constant 10 kcal/mol/Ų (applied only outside 8 Å) was employed. Without this restraint, the shape of the PMF near the channel entrances is ill-defined and the bulk reference value is meaningless because it is determined by the extent of sampling, as it may have been in previous attempts (30). (The influence of the cylindrical restraint can be, and is, rigorously accounted for in the present analysis.)

The one-dimensional PMF and convergence
Initial configurations for the simulations to calculate the PMF were chosen by searching a 4-ns sample of unbiased MD trajectory for a frame in which the ion exists very close to the reaction coordinate and the center of the window (within 1 Å). When no configuration was found with a K⁺ ion near the center of a window, as was the case deep within the channel, a nearby water molecule was exchanged with the outermost ion. The trajectory was searched until a water oxygen atom was located close to that point, and that water molecule was exchanged with the K⁺ ion furthermost from the channel.

We calculated the PMF W(z) using umbrella sampling (56). This requires a set of equally spaced simulations {i}, biased by window functions that hold the ion near positions along the z-axis. We simulated 101 independent windows, defined by harmonic potential functions, positioned at 0.5 Å increments in z = (–20, +30) Å. Harmonic potentials have a force constant 10 kcal/mol/Ų, chosen to ensure overlap of neighboring windows. For each of the 101 windows, equilibration was performed for 80 ps before 1–2 ns of trajectory generation on separate CPUs. Ionic distributions were unbiased using the weighted histogram analysis method (WHAM) (57), which consists of solving coupled equations for the optimal estimate for the unbiased density (z). Strict attention was paid to the convergence of these WHAM equations. Achieving convergence in the free energy constants (to, say, 0.001 kcal/mol), does not guarantee a comparable convergence in the PMF, as further iterating the WHAM equations can lead to many kcal/mol changes in the PMF. To guard against this, each 100 iterations we check every point in the PMF for convergence to within 0.001 kcal/mol. This is a very strict criterion and usually requires >10,000 WHAM iterations. A total of 2 ns of trajectory was generated for each of the 81 windows between z = –20 to + 20 Å and 1 ns for each of the 20 windows from z = + 20.5 to + 30 Å. The one-dimensional PMF, W(z), shown in Fig. 4 (defined only between the dotted vertical lines) reveals much about the permeation process, including a high central barrier and local free energy minima throughout the channel.

View larger version (13K): [in this window] [in a new window]   FIGURE 4  One-dimensional PMFs from simulations with the CHARMM27 force field (without artifact correction). (A) PMFs for all 40 x 50 ps blocks in the total of 2 ns simulation/window. Panels B and C reveal asymmetry in the PMF for 1 and 2 ns per window, respectively, by plotting with the mirror image about z = 0 (dashed curves). The symmetrized PMFs for 1 and 2 ns are shown in panel D. Broken vertical lines at |z| = 15 Å indicate that the one-dimensional PMF is not defined approximately beyond those points.

We computed symmetrical PMFs by applying the WHAM equations to the biased ion-density distribution after creating duplicate windows on opposite sides of the channel. The second test of convergence is thus the difference in the symmetrized 1- and 2-ns PMFs (a temporal convergence). Fig. 4 D shows the PMFs following symmetrization for both the 1 and 2 ns calculations. The PMF is well converged with a maximum deviation near the center of 0.7 kcal/mol (lower in the 2-ns PMF). The average deviation between the two curves is 0.30 kcal/mol, which becomes our estimate for the convergence error in PMF calculation. We conclude that 1 ns per window is sufficient to obtain a well-converged, symmetrized PMF, and we use 1-ns-per-window simulations for all subsequent PMF calculations in this article. The symmetrized CHARMM27 one-dimensional PMF (after 2 ns) reveals a central barrier of 11 kcal/mol with respect to the binding site. This barrier height should be compared with previous estimates (29,30,32) of the order of 15 kcal/mol. There is a deep outer binding site at z = 11.3 Å, and a more shallow inner binding site at 9.7 Å.

As is evident from Fig. 4 A, attempts to estimate the height of the central barrier from a short simulation will result in considerable uncertainty. Though, as noted above, the one-dimensional PMF is not rigorously defined outside the channel, the variation in the value at the channel center relative to ±20 Å provides a measure of the PMF convergence. Thus, when comparing all 40 possible 50-ps blocks in the 2-ns total simulation per window (Fig. 4 A), the central barrier varied by 17 kcal/mol—from as little as just 1.2 kcal/mol to as much as 18.2 kcal/mol. The same analysis for the set of 100-ps nonsymmetrized PMFs leads to values ranging from 5.5 to 18.2 kcal/mol. This range drops to 9.6, 5.8, 5.8, and 4.0 for simulations lasting 200, 300, 500, and 1000 ps/window, respectively. The decreasing range (or uncertainty in the height of the central barrier) with increasing simulation is comforting, but it is important to establish an independent test. In the following section, the barrier height from the full two-dimensional PMF (a well-defined quantity) will be computed and compared to an independent free energy perturbation (FEP) calculation.

Two-dimensional free energy landscape
The range of validity of the one-dimensional PMF is within the bounds of the ion channel. Thus, a complete description of permeation requires a three-dimensional PMF, which can be approximated by a two-dimensional PMF, W(z, r), using radially symmetric configurational integration. Other coordinates may describe slowly varying degrees of freedom in the system. For example, the dipole moment, µ, of the single-file water column was shown, via W(z, µ), to provide a barrier to permeation that is not a function of ion position (15). Equilibrium distributions, biased in z, involving a secondary variable r, may be readily unbiased to produce a two-dimensional PMF W(z, r). Once the WHAM equations (57) have been iterated to convergence, may be unbiased via a straightforward extension of the one-dimensional equation.

The resulting two-dimensional PMF (Fig. 5 A), reveals the position of binding sites at the channel entrances and the scale of the free energy barrier experienced by the permeating ion relative to the entrances. It also reveals the extent of lateral ion motion of the ion. Because the two-dimensional PMF is determined in the laboratory frame, lateral movement of the ion relative to the channel, combined with channel tilting, lead to fairly broad free energy wells. This tilt has an impact on the shape of the one-dimensional PMF of Fig. 4 because the ion experiences a greater radial range, with lower free energies, near the channel entrances (raising the barrier with respect to the binding sites compared to the two-dimensional PMF). The action of the cylindrical constraint in the bulk is evident in this graph. It is also evident that the free energy surface is becoming flat away from the channel. Because this two-dimensional PMF is determined only up to a constant, and we need a bulk reference to establish the zero of this free energy surface. To find the correct bulk reference for the two-dimensional PMF, an additional 4-ns simulation (in the absence of a window biasing potential) was used to calculate the bulk ion density and thus the bulk limit, via where is the K⁺ density far from the channel (30 < |z| <35 Å). The result is shown in Fig. 5 B. This surface extends out to |z| = 35 Å and reveals a flat PMF away from the channel. The values W_bulk(z, r) and W(z, r) were then matched by linearly interpolating over the range 25 < |z| < 29 Å, for r < 8 Å with the resultant two-dimensional PMF shown in Fig. 5 C. Ion free energies are now known at all positions relative to the bulk. The outer binding site at z = 11.3 Å, can be seen to be –3.2 kcal/mol relative to the bulk. In the narrowest part of the channel, an ion experiences a barrier of 7.2 kcal/mol relative to the bulk.

View larger version (49K): [in this window] [in a new window]   FIGURE 5  Creating a two-dimensional PMF. (A) The PMF obtained from two-dimensional unbiasing of equilibrium distributions from umbrella sampling simulations. (B) The PMF obtained by analysis of bulk ion densities. (C) A blend of panels A and B using linear interpolation, similar to that in Allen et al. (15), but with extended range.

A comparison of reaction coordinates
While a coordinate parallel to the membrane normal is necessary for conduction calculation, the free energy surface governing ion permeation also can be studied along an instantaneous channel CoM-axis to reveal in more detail the ion-protein-water interplay during permeation. Such a PMF is interesting for describing the free energy of the ion relative to the ion channel ß-helix. We chose the time-varying vector passing through the center of mass of monomer 1 and monomer 2 of the gA dimer to create a different PMF, based on umbrella sampling along a coordinate that is the ion distance from the center of mass of the channel dimer, projected onto the axis connecting the two centers of mass of the monomers. As before, we apply a lateral constraint in the form of an 8 Å flat-bottom cylinder centered on the instantaneous channel CoM-axis to ensure a well-defined region of sampling outside the channel.

We expect the two PMFs to differ because the channel on average tilts 12° with respect to the bilayer normal (15), which contributes to a radial displacement of ions away from the channel center. Fig. 6 shows the average tilt of the channel as a function of ion position (averaged over 1-ns umbrella sampling for each window). The average tilt decreases steadily as the ion moves from the channel center (where it experiences a maximum of 16°) to the entrances. The decreased tilt near the entrances can be rationalized by considering the interaction between the ion and the bulk electrolyte as a function of depth. The attractive force between an ion in the pore and the bulk solution and membrane interface is stronger when the ion is nearer either interface (58). As a result, the ion-interface interactions will exert a torque on the tilted channel, which will tend to reduce the tilt when the ion is further from the channel center (closer to the interface), as observed in Fig. 6.

View larger version (14K): [in this window] [in a new window]   FIGURE 6  Average channel tilt angle (from average cosine of the angle separating the channel axis and membrane normal, z) as a function of ion position z for the small, one-shell, system. Results have been symmetrized by averaging windows on each side of z = 0. Error bars are not shown for clarity. The average standard deviation in ion position is 0.25 Å and that of tilt angle is 3.6°.

View larger version (63K): [in this window] [in a new window]   FIGURE 7  Two-dimensional PMFs in the pore region for two different reaction coordinates: the position along the instantaneous channel axis (monomer center-of-mass to the monomer center-of-mass) (A), and the projection of the distance to the center-of-mass onto the membrane normal, the z-axis (B). Each surface is based on a calculation with a symmetrized biased density. Only the first 1 ns of simulation for each window is included.

The one-dimensional PMF from this instantaneous CoM-axis umbrella sampling calculation is shown in Fig. 8, together with the 1-ns symmetrized PMF for the original z-coordinate PMF of Fig. 4 C. Not surprisingly, this one-dimensional PMF along a channel-axis vector is not that different to the PMF along the z-vector. The reason for this is that the lateral displacement corresponding to the average tilt is <0.2 Å near the binding sites. One would expect that the depth of the binding sites will change slightly: because of the greater lateral displacements in the z-vector frame (away from z = 0) relative to the instantaneous channel CoM-axis frame, the binding sites should be (a little) deeper (relative to the center of the channel) in the original fixed frame PMF along the fixed z-vector. This is the case in Fig. 8, although the difference is just a fraction of a kcal/mol. There also are small differences in the center of the channel where the original fixed-frame PMF experiences a slightly higher barrier. Given the uncertainties of the PMF and the nontrivial relationship between tilting and the PMF, the origin of the small differences in the PMFs remains unclear.

View larger version (11K): [in this window] [in a new window]   FIGURE 8  One-dimensional PMFs from simulations with the CHARMM27 force field (without artifact corrections) using two different reaction coordinates: the position along the instantaneous channel axis (solid) and the projection of the distance from the center of mass onto the membrane normal, z (dashed). Each curve has been symmetrized (via biased density) and only the first 1 ns of simulation for each window is included.

Correcting for simulation artifacts
A spurious destabilization of the ion is caused by the finite size and the periodicity of the system. In addition, in current MD force fields the hydrocarbon chains of the lipid molecules are nonpolarizable, meaning that they have an effective dielectric constant of 1 (62), which is quite different from the value deduced for the bilayer hydrophobic core (63,64) and measured for bulk hydrocarbons (65), 2. These artifacts can be approximately corrected using a continuum electrostatic approximation (15,66) utilizing trajectories to average over protein and single-file water configurations. We estimate these corrections for an ion near the channel axis and apply to the one-dimensional PMF calculations (with the two-dimensional PMFs remaining uncorrected). However, as we shall show in the following section, all comparisons with experimental measurements can be formulated using the one-dimensional PMF.

To implement the corrections, we establish a free energy cycle (Fig. 9). The corrected PMF W_corr(z; _m = 2; L = ), with the correct membrane dielectric constant, _m = 2, and for an infinite membrane with no periodic effect, may be obtained from the MD-PMF W_MD(z; _m = 1; L = L₀), with respect to some reference position far from the channel, z', where the dielectric constant of the membrane is _m = 1 and the periodic length is L = L₀ in the xy-plane. The relation may be written as

(1)

This may be rewritten as

(2)

where W_MD(z) and W_corr(z) are defined relative to z',

(3)

and

(4)

View larger version (21K): [in this window] [in a new window]   FIGURE 9  The free energy cycle illustrates the sequence of correction calculations required for a PMF calculated with finite system size and nonpolarizable membrane. The gA channel is shown as yellow; high dielectric ( = 80) bulk water as blue; membrane core with = 1 as white; membrane core with correct hydrocarbon dielectric constant ( = 2) as gray; pore water molecules as red (O) and white (H) circles; and K+ as green circles with + sign.

Periodic boundaries are orthorhombic in the PBEQ algorithm whereas they are hexagonal in the MD-PMF. We chose to match the area of the unit cell so as to equate the number of channels per unit area. To do so, the orthorhombic xy-translation vector for PBEQ images was set to 29.9 Å for the one-shell system, instead of 32.1 Å for the hexagonal boundaries. In the three-shell case the length was set to 57.6 Å instead of 61.9 Å. To test the significance of this approximation, we performed calculations with explicitly included hexagonal and cubic periodic xy images of the MD averaged protein and ion in the PBEQ module of CHARMM. The results for hexagonal and cubic images differed by just 0.19 kcal/mol, suggesting that the corrections are not significantly affected by this choice.

To make these profiles more readily applicable as corrections (requiring interpolation onto a finer grid), we fitted the MD-averaged Poisson solutions with solutions for a single structure using optimized dielectric constants of protein (and channel water). This was done with a single MD-averaged structure of the protein with implicit water inside the pore. The MD-averaged protein structure was obtained by constrained minimization of the PDB:1JNO structure to MD-averaged heavy atom positions. The dielectric constants _p = _c were optimized to obtain the best numerical fit the MD data.

Estimates, obtained from averages over MD trajectories, show that correcting for the spurious destabilization leads to a –1.6 kcal/mol correction (at the channel center, relative to the bulk) for the one-shell system (Fig. 9 A). Fits using values of _m = _p = 1 and 2 envelop the MD solutions, and the optimal choice was found to be _p = _c = 1.25 (Fig. 9 A). For the three-shell system, the correction for periodicity had a maximum amplitude of 0.03 kcal/mol. The spurious destabilization has decayed almost completely when the system size is doubled.

Correcting for the effect of the dielectric constant of the hydrocarbon chains leads to a further –2.1 kcal/mol stabilization of the ion in the center of the channel, as revealed in Fig. 9 B. The MD-averaged data were fitted with a Poisson solution using _p = _c = 1.75. To test the extent of the effect of polarizability of the phospholipid hydrocarbon tails on the stability of the ion, we previously carried out calculations using Drude oscillators (68) to describe the hydrocarbon chains (15). A similar approach using a grid of polarizable point-dipoles was used by Åqvist and Warshel in 1989 (69). In this case, the change in the stabilization of the ion was –3.6 ± 0.3 kcal/mol, which can be compared to our estimate G_diel(0) from the Poisson solutions (–2.1 kcal/mol). The continuum electrostatics approximation captures the effect of the polarizability, but may underestimate its magnitude.

It is also possible to estimate the effect of high electrolyte concentration (70) by approximating the effect of reducing the concentration from 1 M to a level of 0.1 M that better corresponds to the single-ion regime

(5)

using the Poisson-Boltzmann equation (see (70)). As expected, the effect of reducing the ionic concentration from 1 M to 0.1 M is a small additional stabilization of –0.2 kcal/mol at the channel center, as seen in Fig. 10 C.

View larger version (10K): [in this window] [in a new window]   FIGURE 10  Corrections applied to the one-dimensional PMF to correct for simulation artifacts. (A) Poisson size correction (one shell of lipids infinite bilayer); (B) Poisson membrane dielectric constant correction (m = 1 2); (C) Poisson-Boltzmann concentration correction (1 M 0.1 M). Data points represent calculations which are ensemble averages of Poisson solutions using a set of MD protein and channel-water coordinates. Short-dashed curves in panels A and B are corrections that use a single MD-averaged structure with the dielectric constant of the protein and channel water of 1. Long-dashed curves are corrections that assume the dielectric constant of the protein and channel water to be 2. The solid curves in panels A and B employ protein/channel-water dielectric constants of 1.25 and 1.75, respectively, as a fit to the calculated MD averages. In panel C, the dielectric constant was assumed to be 1.5.

View larger version (25K): [in this window] [in a new window]   FIGURE 11  (A) One-dimensional CHARMM27 PMFs (from 1 ns simulation per window) before (dashed) and after (solid) artifact corrections. The dotted curve shows the corrected PMF with scaled hydrocarbon dielectric correction (see text) to be referred to as the "Drude-corrected" CHARMM27 PMF. (B) comparison of PMFs, before and after corrections, between small (20 lipids, red) and large (96 lipid, blue) membranes. The PMFs have been matched at z = 20 Å. PMFs from 2 ns/window simulation (not shown) experience barriers that are 0.7 kcal/mol less than the plotted 1 ns/window PMFs (with or without Drude-oscillator-scaled corrections).

To test the accuracy of our size correction of Fig. 9 A, and to further evaluate the consistency and convergence of the results, we obtained a PMF with the larger three-shell membrane system. Fig. 11 B compares the one- and three-shell PMFs (using 1 ns of simulation in both cases) before corrections as dashed curves. The two PMFs differ by 1.6 kcal/mol at the channel center, a difference that is similar to the size correction for the one-shell system (Fig. 10). The solid curves in Fig. 11 B show our final results with all corrections (size, dielectric of the membrane, and concentration). The two PMFs now look very similar and differ on average by 0.28 kcal/mol. The correction obtained with Poisson solutions therefore accounts for reduced membrane size.

Channel binding and conductance
Experimental observables may be predicted from the calculated MD-PMFs to ascertain their agreement with electrophysiology. Here we outline the strategy by which we use the equilibrium free energy calculations to predict macroscopic observables and apply these methods using our corrected PMF.

Equilibrium dissociation constants
The equilibrium single-ion dissociation constant K_D can be expressed either in terms of the three-dimensional one-ion three-dimensional PMF, or in terms of the one-dimensional PMF, W(z), (see (53)) as long as the sampling of the lateral motion (i.e., with a cylindrical restraint) and corresponding offset constant are incorporated correctly (15). The dissociation constant for the channel (encompassed by –15 < z < 15 Å) from the corrected 1-ns PMF (Fig. 11) is 0.30 M. The 2 ns/window PMF for the CHARMM27 force field from Fig. 4 d, after artifact corrections (not shown) yields a dissociation coefficient of 0.21 M. This estimate is in within the range of experimental values determined from NMR and conductance studies: 0.017 M (71) (measured in dodecylphosphocholine micelles); 0.019–0.73 M (72) (in aqueous lysophosphatidylcholine dispersions); 0.035 M (54) (in glycerylmonoolein bilayers); and a recent estimate of 0.07 M (25) (in diphytanoylphosphatidylcholine/n-decane bilayers). When integrating over the regions defining the individual sites, ions will bind to the outer (10.2 < z < 12.5 Å) and inner (6.9 < z < 10.2 Å) sites with dissociation constants of 0.83 M and 3.6 M, respectively. Previous studies have found mixed results for the most stable position for a Na⁺ (29,73,74). Experimentally, the major cation binding site for both Na⁺ and K⁺ is the inner site (60,71). Significant K⁺ binding may also occur further from the channel center, as Tian and Cross (59) show that the Na⁺-induced chemical shift of Leu-10 to Trp-11 linkage exceeds that of Leu-12 to Trp-13, whereas the reverse is the case for K⁺.

Channel conductance
To ascertain the magnitude of the current that can pass through the channel, the net stationary flux of ions across the channel can be calculated using a one-dimensional one-ion pore Nernst-Planck theory of Levitt (75) (see also (2) for application to MD-PMF). The ion flow is governed by the total PMF, W_tot, which can be expressed as a sum of the equilibrium PMF W, dominated by local molecular interactions, and the interaction of atomic charges with the transmembrane scalar potential, _mp; see Eqs. 38 and 39 in Roux (53). The one-ion pore NP theory of Levitt (75) provides a convenient route to quickly obtain a rough estimate of the maximum conductance given a one-dimensional PMF. There are alternative approaches that, however, must be considered with caution. For example, utilizing the one-dimensional MD-PMF in any three-dimensional model of permeation as done by Corry and Chung (76) requires additional assumptions and approximations about the width, the length, and the shape of the pore, as well as the absolute value of the free energy in order to construct a complete three-dimensional PMF (the ambiguities in such an operation are well illustrated by considering Fig. 5).

The charge distribution changes considerably as the ion moves, with most of the change being due to the dipole moment of the pore water (15). To gauge the contributions from the coupling to the system dipole moment to the trans-membrane potential, we express the total PMF as a cumulant expansion in the dipole moment of the system (see Eq. 36 in (53),

(6)

where µ(z) is the deviation from the mean value of the dipole moment of the pore system as a function of the ion position. The value _mp is estimated by solving the modified Poison-Boltzmann equation (2,53) for an MD-averaged gA structure in a 25 Å membrane slab of dielectric constant 2, using 1 M KCl salt in the bulk regions with dielectric constant 80. Poisson-Boltzmann solutions for _mp, together with histograms of the dipole moment of the single-file water column obtained from all umbrella sampling trajectories (with symmetrization), have been used to estimate the perturbations to the equilibrium PMF of Eq. 6. The charge distribution within the pore region is represented by an expansion in the dipole moment, dominated by the single-file water column in Eq. 6, and the dielectric constant is set to 1 inside the protein and channel regions. We also evaluate the value where L_p is the length of the pore region. The value L_p is equal to the diameter of the sphere defining the single-ion region (28 Å) and is within the range where the one-dimensional PMF is meaningful (|z| < 15 Å).

Conductance calculations are complicated by the coupling of the dipole to the membrane potential, with the linear term in µ(z) influencing the flow of ions and the quadratic term creating an additional barrier as a consequence of the bimodal nature of the dipole distribution in the bulk (15). Fig. 12 A shows the membrane potential term q_ionmp of Eq. 6, increasing linearly across the channel. Fig. 12 B shows the first correction of Eq. 6 (due to the linear term in the dipole moment), and Fig. 12 C the second correction (due to the quadratic term). The sum of all three terms representing the coupling of the system charge distribution to the applied voltage is shown in Fig. 6 D (with curves from Fig. 6 A superimposed as thin curves). When 100 mV is applied across the membrane (solid curves in Fig. 12), the linear and quadratic dipole corrections have maximum amplitudes of only 0.15 and 0.04 kcal/mol, respectively. Because of the different voltage dependencies, the quadratic term becomes substantial at larger applied potentials. At 500 mV (dash-dot curves), the corrections approach 1 kcal/mol and, as the voltage applied across the membrane is increased, these corrections may have considerable effect on the conductance.

View larger version (15K): [in this window] [in a new window]   FIGURE 12  Terms in the cumulant expansion Eq. 6. (A) The membrane potential for 100 mV (solid), 250 mV (dashed), and 500 mV (dash-dot) applied potential difference from solutions to the modified Poisson-Boltzmann equation with MD-averaged structure. The linear dipole (B) and quadratic dipole (C) cumulant corrections are shown in panels B and C, respectively. The sum of all terms for each voltage difference are plotted in panel D as thick curves. The thin curves superimposed in panel D are the qionmp terms from panel A.

The diffusion coefficient, D(z), required to calculate the flux may be estimated in different ways, where the general requirement is a method that separates the local dissipative forces from the systematic forces arising from the PMF. We employ a method based on the Laplace transform of the velocity autocorrelation function, using an analysis of the generalized Langevin equation for an harmonic oscillator (80), with the expression for the diffusion constant given in Crouzy et al. (81). To estimate the value of the limit as the Laplace transform variable s 0, we linearly extrapolate from the range 15 s 35. Fig. 13 shows that the (symmetrized) axial ion diffusion coefficient D(z) is approximately two-thirds of bulk diffusion coefficient within the channel.

View larger version (13K): [in this window] [in a new window]   FIGURE 13  K+ ion diffusion profile. Calculated values of the axial component of the ion diffusion coefficient, for each window simulation, are drawn with a solid line. All values have been symmetrized (unlike Fig. 6 in (15)) and scaled relative to the calculated bulk value of 0.37 Å2/ps. The fit (dashed line) is a sigmoidal function.