Search tips
Search criteria 


Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
Biochim Biophys Acta. Author manuscript; available in PMC 2010 June 1.
Published in final edited form as:
PMCID: PMC2701196

Computational Studies of Gramicidin Permeation: An Entryway Sulfonate Enhances Cation Occupancy at Entry Sites


The impact on the cation-transport free-energy profile of replacing the C-terminal ethanolamine in the gramicidin A channel with a taurine residue is studied using molecular dynamics simulations of gramicidin A (1JNO) embedded in a lipid bilayer (DMPC) with 1 mol/kg NaCl saline solution. The potential of mean force for ion transport is obtained by umbrella sampling. The presence of a negatively charged sulfonate group at the entrance of the gramicidin channel affects the depth and the location of the binding sites, producing a strong attraction for the cations in the bulk. The potential of mean force by the sulfonate acting directly through electrostatics and van der Waals interactions on the test ion is highly modulated by indirect effects (i.e., sulfonate effects on other components of the system that, in turn, affect the ion free-energy profile in the channel). Because the “entry” sites are located symmetrically at both entry and exit of the channel, the deeper free-energy wells should inhibit exit. Given that the channel has increased conductance experimentally, the simulation results suggest that the channel conductance is normally entry limited.

Keywords: Ion channel, Potential of mean force, Sulfonate parametrization, Water layering, Molecular dynamics simulations

1. Introduction

Molecular dynamics (MD) simulation is an important tool for studying both biological and chemical systems. For example, MD simulations have been used to study several ion channels, such as those formed by gramicidin A (gA) [14] and the influenza A M2 transmembrane channels [57].

A notion was developed that calcium interactions would be enhanced in a narrow channel with its selectivity-filter because of the competition of enthalpic electrostatic and available volume effects [8]. This notion was further developed with implicit water using Monte Carlo (MC) [911] and Brownian dynamics [12] simulations. MD simulations with explicit water with a non-polarizable force field yielded little selectivity for Ca+2 passage over Na+ passage in a calcium-like channel [13]. But, the inclusion of polarizability in an MC study yielded strong Ca+2 selectivity for the Ca+2-channel selectivity-filter and Na+ selectivity for the Na+-channel selectivity-filter [10]. Furthermore, comparison of Brownian dynamics studies and MC simulations, both using implicit solvation, demonstrated that the charged moieties in the selectivity-filter must be exposed inside the selectivity-filter rather than buried in the walls to attain the high degree of selectivity revealed through electrophysiological experiments [10, 11].

Because sodium and calcium channels differ dramatically in their physiological functions and have charged selectivity-filters [14], we wanted to explore the question of how will a charged selectivity-filter affect the channel selectivity? The simplicity of experimenting with and modeling gA analogs makes it possible to explore that question. Here we focus on the effects of single ion–charge interactions, for which the narrowness of the gramicidin channel is ideal. Our question is more fundamental than and helps lay the groundwork for the more complex question of charge–space competition in the larger cation channels mentioned above.

The gA channel is a polypeptide that dimerizes head–head to form a cation channel in lipid bilayer membranes. The gA channel is composed of two right-handed [15] single-stranded hydrogen-bonded β6.3-helices [16, 17] with the amino acid sequence: formyl-L-Val-Gly-L-Ala-D-Leu-L-Ala-D-Val-L-Val-D-Val-L-Trp-D-Leu-L-Trp-D-Leu-L-Trp-D-Leu-L-Trp-ethanolamine [18]. The gA channel has a simple well-known structure with well-defined functions [19, 20]. A minor variant introducing a charge at the channel entryway is the Taurine gramicidin A (TgA) channel [21, 22]. In the TgA channel, the ethanolamine residue at the C-terminus of the gA channel, which lies nearest to the lipid–water interface at each end of the channel, is replaced with taurine (i.e., the hydroxyl group of the C-terminus is replaced with the sulfonate group). Thus the TgA model is a good candidate for exploring how a charged selectivity-filter will affect ion–channel interactions.

In experiments with the TgA channel, the single-channel conductance for the TgA channel in 0.1–1.0 mol/kg univalent cations has been found to be higher than that for the gA channel [21]. The higher single-channel conductance in the TgA channel compared with that in the gA channel is consistent with other results using other charged gA analogs [2226]. That is, the presence of a negatively charged group at the entrance of the gramicidin channel causes an increase in the cation transport rate at low ionic concentration (< 1.0 mol/kg), whereas cationic adducts decrease cation conductance [2226]. The increased single-channel conductance has been explained as the result of an electrostatic effect of the fixed negative charges on the bath cation concentration near the entrance of the channel [23]. At higher ionic concentrations (≥ 1.0 mol/kg), the conductance difference is diminished, as might be expected if shielding of the sulfonate charge mutes the bath cation concentration enhancement.

Although the negatively charged sulfonate group at the entrance of the channel should attract the cations in the bath towards the vicinity of the channel entrance, it is not clear where this would occur, how much would involve direct contact with the sulfonate, and how an increased bath ionic strength would affect these factors. Furthermore, it is not clear whether the dominant effect of the sulfonate charge at each end of the channel would be to enhance cation entry into the channel or translocation through the channel, although it appears that the net result of these two effects must supercede electrostatic inhibition of cation exit and departure at the far end of the channel.

Here, we compare the one-dimensional one-ion free-energy profile for cation transport through the gA and the TgA channels, extracted from umbrella-sampling simulations, to explore some of these specific questions, with the more general goal of illuminating the roles of negatively charged side chains in other ion channels such as the voltage-gated sodium and calcium channels with their polyvalent selectivity-filters [27, 28].

2. Methods

2.1. Model System and Parameters

The initial system consisted of a gA dimer, with the 1JNO structure [29], embedded in 40 1,2-Dimyristoyl-sn-glycero-3-phosphocholine (DMPC) molecules as a lipid bilayer, as well as 2366 water molecules (with the TIP3P model) [30], and 1 mol/kg NaCl. The 1JNO structure is a relaxed structure, based on solution state NMR results for dimers in sodium dodecyl sulfate micelles [29], and consistently reproduces the solution and the solid-state NMR properties [31]. The terminal residues, ethanolamine and formyl, were given their partial charges and force field parameters based on other similar atomic groups in the force field. The initial crystal-like structure was constructed, using the VMD software [32], in a box with initial dimensions of 40 × 40 × 90 Å in x, y, and z dimensions, respectively.

The initial TgA system was modeled in the same way as the gA system except for the taurine residue. Because the only difference between the TgA dimer and the gA dimer is the C-terminus residue, the two ethanolamine residues (NHCH2CH2–OH) in the 1JNO dimer structure were replaced with taurine (NHCH2CH2–SO3) to produce the TgA dimer. The taurine residue was split into two parts: the NHCH2 group and the CH2SO3 group. The first part was assigned the same parameters as the glycine residue, whereas the methylsulfonate parameters were obtained using the Extensible Computational Chemistry Environment (ECCE) [33] after performing quantum mechanics calculations using NWChem [34, 35]. A (CH3CONHCH2–CH2SO3) molecule was built in ECCE and geometrically optimized at the RHF level of theory with the 6-311++G(3df,3pd) basis set, with Pople(3df,3pd) as a polarization function, and with Pople-style as a diffusion function. The bond lengths for methyl sulfonate and its connections in the taurine residue were derived from the resulting structure. The atomic partial charges were derived using electrostatic potential fitting while constraining the two hydrogens in the CH2SO3 group and all atoms in the CH3CONHCH2CH2 group with the atomic partial charges from the empirical force field. The taurine parameters are listed (Table 1). To produce electroneutrality for the TgA system, two Cl were removed.

Table 1
Taurine Parameters

2.2. System Preparation

The total simulation time for system preparation was 25 ns, with different periodic boundary conditions and harmonic restraints (Table 2). The preparation included conjugate-gradient energy minimization of the initial configuration, followed by heating from 0 to 303.15 K, equilibration in the NVT ensemble, three-stages of annealing (40 ps of heating to 1000 K followed by 80 ps of cooling to 303.15 K), and progressive release of restraints on the channel over multiple stages of equilibration (Table 2). This preparation allowed relaxation of the channel and disordering the lipid molecules, ions, and water molecules starting from the crystal-like structure.

Table 2
Summary of System Preparation Approach

Both the lipid deuterium order parameters and the lipid surface area are sensitive to the changes in the bilayer structures due to many factors, such as the lipid hydration level, number of lipid molecules, temperature, statistical ensemble, surface tension, and saline-solution concentration [3643]. Surface tension is the force applied per unit length to a filament required to stretch a liquid film [44, 45]. Mathematically, the surface tension per interface (γ) can be calculated from the pressure tensors, averaged over the central cell containing a lipid bilayer system normal to the z direction, as follows [39]:


where Lz is the length of the unit cell, Pzz is the component of the pressure tensor normal to the bilayer surface, Pxx and Pyy are the tangential components of the pressure tensor, and the angle brackets denote an average over time. For simulating a pure lipid bilayer in pure water, the appropriate surface tension has been evaluated to be within the range of 30-45 dyn/(cm·interface) [39, 41, 46, 47]. However, with a more complex lipid bilayer system, such as one that includes salt and polypeptide, the surface tension is unknown.

We carried out a test for the suitable surface tension using the Equilibration 3 stage conditions. The test was done by fixing the channel in the original 1JNO structure to avoid channel-surface perturbation effects on the lipid bilayer properties [48], applying two planar harmonic restraints on the lipid phosphorus atoms [kf = 10 kcal/(mol·Å2)] at z = ± 17.5 Å to agree with the experimental head–head or phosphate–phosphate thickness of 35.3 Å [49, 50], and using 2.5 unit increments of the surface tension in the range of 0–75 dyn/(cm·interface). Based on accuracy and precision in the lipid molecular surface area, which was consistent with the experimental DMPC molecular surface area [50], the optimal surface tension was found to be 40 dyn/(cm·interface).

The Equilibration 3 stage was then carried out with γ= 40 dyn/(cm·interface). In the Equilibration 3 stage, the channel was fixed in the original 1JNO structure, and two planar harmonic restraints were applied on the lipid phosphorus atoms [kf = 10 kcal/(mol· Å2)] at z = ± 17.5 Å (Table 2). This stage was run for 20 ns, which was enough time to obtain stable lateral dimensions for the system under study [51]. This stage gave, in the last 10 ns, a lipid surface area of 60.5 ± 0.8 Å2 and 60.4 ± 1.0 Å2 in the gA and the TgA systems, respectively, after subtracting 250 Å2 as a gA surface area [52], which we did on the reasonable assumption that theC-terminus replacement has negligible effect on the TgA surface area compared to gA surface area. The purpose of this stage was to prepare a realistic lipid environment for the polypeptide channel.

During the subsequent equilibrations, the restraints on the channel and the lipid phosphorus atoms were progressively reduced without significant changes in the system structure. The simulations in the NPAT ensemble (Equilibrations 4–10) utilized a tetragonal box having a lateral dimension of 38.2 Å for both the gA and the TgA systems, the average value from the second 10 ns of the Equilibration 3 stage. The similarity of the lateral dimension values, obtained independently for the two systems, as well as the calculated lipid surface area values demonstrate the robustness of the equilibration procedure and the stability of the lipid-polypeptide-electrolyte structures after electrolyte homogenization and lipid melting. Both the Equilibration 10 stage and the routine 10 ns production run (done in parallel with the umbrella sampling simulations) were performed with the channel center-of-mass (COM) as the only restraint. Considering that the lipid molecules immediately adjacent to the gA or the TgA channel have different properties than those farther away from the channel due to channel-surface perturbations [48], and that 16 lipids per leaflet has previously been shown to have sufficient bulk-like lipid surroundings for the gA dimer [3], our final system has a reasonable size.

2.3. Molecular Dynamics Simulations

Simulations were performed with NAMD software [53], version 2.6, and using version 31 of the CHARMM force field [5459]. Trajectory files were written every 1 ps for both the 10 ns test-ion-free production run and the umbrella sampling production runs using a time step of 2 fs and the NPAT ensemble. A smoothing function was applied to both the electrostatic and the van der Waals forces at a distance of 8 A with a switching cutoff distance of 12 Å, a pair list distance of 14 Å, and all bonds with hydrogen being rigid using the ShakeH algorithm. The non-bonded interactions parameters (switching, spherical cutoff, and pair list distances) were carefully chosen to achieve realistic macromolecular simulations [60]. The Particle Mesh Ewald (PME) algorithm, implemented to add a reciprocal space long-range electrostatic energy to the limited non-bonded component, was used with an interpolation function of order 5, and PME box grid size of 64, 64, 128 in x, y, and z direction, respectively. Langevin dynamics were used to maintain the temperature at 303.15 K with a damping coefficient of 5 ps−1, whereas the pressure in the z direction was maintained at 1 atm using the Nosé–Hoover Langevin piston method with an oscillation of 1 ps and a damping time of 0.1 ps [61, 62].

Analysis scripts were composed using VMD and Tcl [63] commands. The deuterium order parameters of the deuterium-labeled segment were calculated from [42]:


where θi is the angle between the ith C–H vector and the z axis (bilayer normal), and the angle brackets denote an average over both acyl chains from all lipid molecules, all C–H vectors associated with the ith carbon, and time. The deuterium order parameters have a range from −0.5 to 1.0 and reflect the lipid tails orientation (i.e., a −0.5 value represents fully ordered lipid tails, whereas a 1.0 value represents fully disordered lipid tails).

The root-mean-square displacement (RMSD) reflects the degree of similarity between three-dimensional protein structures. RMSD was computed by measuring the RMSD between equivalent atoms after optimal superposition of the two structures (the protein structure in trajectory frame i and the initial trajectory frame).

2.4. Umbrella Sampling Simulations

All the umbrella sampling simulations were performed using the final frame of the Equilibration 10 stage of the system under study as an initial structure. In these simulations, the channel COM was restrained to within 0.1 Å of the origin [ kfCOM=1000kcal/(mol·2)] and the test ion, selected from the bath ions, was restrained to a plane normal to the z axis [kf = 30 kcal/(mol·Å2)], with the position of the plane being varied along the transport coordinate (“reaction coordinate”). All kf symbols designated force constants (i.e., harmonic spring constants).

Here we use the term “transport coordinate” because no bond-forming or -breaking takes place during ion transport. We don’t view the restraining force on the channel COM as a biasing force for the umbrella sampling, but rather as a structural restraint force designed to compensate for missing components in the system, that is, the rest of the experimental system (planar bilayer chamber, etc), and thus to establish a reference point for the coordinate system.

During umbrella sampling with the test ion, a spherical restraint [kf = 20 kcal/(mol·Å2)] was applied on all other ions to keep them more than 14 Å away from the origin and a cylindrical restraint [kf = 20 kcal/(mol·Å2)] was applied on the test ion to keep it within 5 Å of the z axis. The purpose of the cylindrical restraint was to assure a one-dimensional free-energy profile [2, 6466], whereas the purpose of the spherical restraint was to assure adequate sampling in the bath region for the one-ion free-energy profile [2, 6466].

We used 201 windows with a window width of 0.2 Å to cover the transport coordinate −20 ≤z ≤20 Å. For each window, the z coordinate of the test ion was set to the window position and the x and y initial values of the test ion were set to zero before performing minimization for 3000 steps, heating from 0 to 303.15 K, equilibration in the NPAT ensemble for 1 ns, and production in the NPAT ensemble for 1 ns. We used the weighted histogram analysis method (WHAM) to unbias the umbrella potential on the test ion and extract the free-energy profile [67]. The mathematical details for calculating or extracting the free-energy profile were the same as described previously [68].

The decomposed potential of mean force (PMFi) is calculated by integrating the z component of the mean force left angle bracketFi(zion)right angle bracket applied by the interesting group (i) on the test ion over the transport coordinate [69],


where the PMFi(zo) is an arbitrary reference value, generally assumed to be zero at zo, zw is the window-reference position, and left angle bracketFi(zion) right angle bracket is the average equilibrium force applied by group (i) on the test ion at window w. The force Fi(zion) is extracted via the pair-interaction-calculation method available in the NAMD software package. The pair-interaction-calculation method analyzes the available dynamic trajectory (ours is sampled every 1 ps for the ensemble average) and calculates the forces, van der Waals and electrostatic, from the second group, group (i), on the first group, the test ion. The PMF components are additive, with the total PMF equaling the sum of the PMF components [69].

The potential of mean restraint (PMR) is used as a test for the WHAM results, free-energy profiles. The PMR can be extracted as follows [70]:



where zion is the test ion position in the z direction (sampled every 2 fs for the ensemble average), zw is the window-reference position, and F (zion) is the instantaneous restraint force applied in the z direction to the test ion. The F(zion) is expected to be equal, on average, but opposite to the time average force of the system on the test ion under mechanical equilibrium conditions, which apply in umbrella sampling. Therefore, the PMR should equal the PMF.

The symmetrization of our results was obtained by averaging the free-energy values at positive and negative values of the transport coordinate. The Composite Simpson’s Rule was used as a numerical integration method [71].

3. Results and Discussion

3.1. Preparation and Simulation Methods Stability

The simulating annealing, based on the deuterium order parameters [Fig. 1(a)], was sufficient to produce disordered lipid tails in TgA and gA systems. Our results agree excellently with both the experimental [7274] and the simulation [42] results, except in the region close to the head group, where the order remains high in experiments [7274], but commonly declines in simulation [42]. The RMSD results show reasonable polypeptide stability throughout the 10 ns production run time scale in both the TgA and the gA systems [Fig. 1(b)].

Fig. 1
System stability: (a) Deuterium order parameters for the DMPC tails in the TgA (black line) and the gA (grey line) systems. (b) Root-mean-square displacement versus time for all the TgA (black line) and the gA (grey line) channel atoms. (c) Averaged ...

There is reasonable head group stability, demonstrated by the averaged z-positions of the lipid-phosphorus atoms during the production run, in both the TgA and the gA systems [Figs. 1(c) and (d)]. The C-terminus stability in the TgA system is demonstrated by the z-positions of the taurine-sulfur atoms [Fig. 1(c)], whereas the C-terminus stability in the gA system is demonstrated by the z-positions of the ethanolamine-oxygen atoms [Fig. 1(d)].

The effect of the taurine terminal compared with that of the ethanolamine terminal is shown by the molality profiles in the z direction [Figs. 2(a) and (b)]. Ignoring the high occupancy at (9 ≤z ≤11 Å), there is a higher Na+ occupancy near the mouth of the TgA channel (12 ≤|z|≤22 Å) compared with that of the gA channel. The high occupancy at (9 ≤z ≤11 Å) in the long simulation was spontaneous, the channel having started empty, whereas it was not the case for the shorter umbrella sampling runs, which started and ended with no occupancy of the binding sites by mobile cations.

Fig. 2
(a) Na+ molality profiles in the z direction for both the TgA (black line) and the gA (grey line) systems. (b) Cl molality profiles in the z direction for both the TgA (black line) and the gA (grey line) systems.

The consistency of the results for the TgA and the gA systems (Figs. 1 and and2)2) demonstrates the robustness of the equilibration procedure and the stability of the lipid-polypeptide-electrolyte structures after electrolyte structural homogenization and lipid melting.

Could the interaction between the taurine terminals and the Na+ near the mouth of the channel lead to a deepening of the energy wells at each end of the channel? What is the single-occupancy translocation free-energy barrier in the TgA system compared with that in the gA system? The answer to these questions are discussed in the next section.

3.2. C-terminus Effects on the Na+ PMF Profile

The Na+ PMF profiles in the TgA and the gA systems are shown [Fig. 3(a)]. Generally, the Na+ PMF profile in the TgA system is lower than that in the gA system except at the central-channel region (−4 ≤z ≤4 Å), where it is higher. The TgA channel has binding sites at ±10.5 Å with well depths of −3.0 kcal/mol and a central-barrier height of 3.5 kcal/mol, whereas the gA channel has binding sites at ±9.1 Å with well depths of −1.9 kcal/mol and a central-barrier height of 1.2 kcal/mol. Indeed, the free-energy wells are deeper in the TgA system compared with that in the gA system, but the central-barrier height is higher.

Fig. 3
One-dimensional Na+ PMF profiles as a function of the transport coordinate. (a) The symmetrized Na+ PMF profiles for both the TgA (black line) and the gA (grey line) systems. (b) The non-symmetrized Na+ PMF profiles for the TgA system including different ...

Free-energy profile convergence depends on many factors, such as adequacy of sampling of test ion positions in the umbrella windows, variations in surrounding structures on the umbrella window time scale (which, in turn, could be different for usage of serial processions along the transport coordinate rather than uniform starting structures for each umbrella window), and statistical fluctuations affecting relative positioning of PMF components from adjacent windows. We examined free-energy profile convergence by comparing free-energy profiles obtained with increasingly longer portions of the umbrella window trajectories and examining recovery of the initial reference free-energy at the end of the transport coordinate. The non-symmetrized PMF profiles [Figs. 3(b) and (c)], in TgA and gA systems, show good convergence. That is, PMF profiles based on increasingly complete fractions of the umbrella sampling trajectory (0.5, 0.75, and 1.0 ns) converged to within 0.5 kcal/mol over the entire extent of the integration, and to within 0.5 kcal/mol of the reference energy (0 kcal/mol at z = −20 Å) at the opposite z position (z = 20 Å) [Figs. 3(b) and (c)]. In the demonstration of convergence, the non-symmetrized PMF profiles [Figs. 3(b) and (c)] suggest that the energy uncertainty due to structural variations in the two identical halves of the channel is on the order of ≈ 1 kcal/mol.

One way to illuminate the influence of replacing the terminal ethanolamine residues with the taurine residues on the Na+ PMF profile is by decomposing it into its components, Na+ PMFi [Eq. (3)]. The Na+ PMFi profiles [Figs. 4(a) and (b)] are calculated with respect to the reference z = ±20 Å. The choice of extracting the Na + PMFi profile with respect to this reference, z = ±20 Å, is based on the use of this reference in previous work [3, 51, 65, 66, 75]. Even though our calculations with respect to the reference (z = ±20 Å) yield flat (i.e., bulk-like) Na+ PMF profiles [Fig. 3(a)] at the symmetrical reference points, the Na+ PMFi profiles [Figs. 4(a) and (b)] are not flat at these symmetrical reference points. Therefore, the lack of flatness at the reference points may be considered to represent an error, probably by fixed and canceling amounts, in the Na+ PMFi profiles.

Fig. 4
One-dimensional Na+ PMFi profile as a function of the transport coordinate. (a) The Na+ PMFchannel profiles for both the TgA (black line) and the gA (grey line) systems from the force of the channel on the test ion. (b) The Na+ PMF H2O profiles for both ...

Indirect effects of force field changes can be complex and the compensation among the major system components in producing the PMF profile can be unpredictable. We define the indirect effect as an effect of the modified force field on the free-energy profile mediated by induced changes in the positions or dynamics of other atoms in the system. We illustrate this principle here, focusing on two of the four major components of the system, the channel and the water [Figs. 4(a) and (b)].

There is a higher central-barrier height for the Na+ PMFTgA profile, as well as deeper and shifted free-energy wells compared with those in Na+ PMFgA profile [Fig. 4(a)]. Furthermore, comparing water contributions for the two systems [Fig. 4(b)], the Na+ PMF H2O profile has higher values outside the channel region. The high Na+ PMF H2O profile barrier beyond z = ±10 Å counteracts or partially shields the effect of the SO3 group. The shielding effect of water molecules is more dominant in the TgA system than in the gA system.

These two examples demonstrate the complexity of indirect effects and show that it is difficult to predict the free-energy barriers to transport from intuition based on electrostatics. Therefore, the final height of the central-energy barrier of Na+ PMF profile in the TgA system [Fig. 3(a)] may be influenced by the indirect water, polypeptide, or lipid effects induced by the lack of the polarizability effect in our simulations [2, 3, 69, 76, 77], or other unintended consequences of empirical force field parameters or applied restraints.

The higher central-energy barrier in the TgA means a lower single-channel conductance would be anticipated in the TgA channel compared with that in the gA channel in disagreement with experimental studies [2126]. In fact, the high central-energy barrier in the PMF profile for the TgA system compared with that in the gA system [Fig. 3(a)] might be exaggerated from applying the spherical restraint. That is, the spherical restraint excludes all mobile cations except for the test ion from a sphere of radius 14 Å where the channel is in the center of this sphere. No bath cations are allowed to approach the channel, which prevents the natural phenomenon where by the cations screen the negatively charged group. Hence, when using the spherical restraint, alleviating the influence of the negatively charged group (SO3 group) is only achieved by orientation of water molecules with the water hydrogen atoms pointing towards the SO3 group, and prevents possible shielding effects of the bath ions, which could influence the height of the central barrier. Therefore, the extra height of the central barrier should be viewed with caution until further confirmed. We point out, however, that if the central barrier is higher, it would isolate the entry process from the exit process and lead to the prediction that the observed conductance enhancements are due primarily to improved channel entry rather than improved translocation through the channel. From this, one would predict that heterodimers with TgA at the entry and gA at the exit would have enhanced conductance, whereas those with TgA at the exit should have unchanged or diminished conductance [22].

In the process of analyzing these simulations, we wondered whether the sulfonate group influenced water orientation around the channel entry. In the process we noticed a much broader phenomenon: the water is layered near the bilayer with or without the sulfonate on the peptide channel. This did not seem to influence the free-energy profiles shown above [Fig. 3(a)], but for completeness, and because it may prove to have broad ramifications, we display this interesting structural phenomenon in the last section below.

3.3. Water Molecule Orientation

The orientations of groups of water molecules, both inside and outside the channel, are shown [Figs. 5(a) and (b)]. The color coded average orientation of all molecules in a 0.5 Å slab, averaged over the course of an umbrella sampling trajectory, is plotted as a pixel with a y coordinate equal to the center of the slab and x coordinate equal to the umbrella potential position (i.e., the desired test ion position along the transport coordinate). The angle θ is the angle between the product of the two O–H vectors and the z axis.

Fig. 5
Normalized water molecules orientation as a function of the transport coordinate and the z axis of a particular window frame for both (a) the TgA and (b) the gA systems.

There is a layering of external water molecules in the z direction with the innermost layer being slightly more intense in the TgA system [Fig. 5(a)] than that in the gA system [Fig. 5(b)], presumably as a result of replacing the hydroxyl group with the sulfonate group. Such oscillations are reminiscent of large oscillations in water molecules density observed in simulations within 1 nm of a rigid surface [78], and may play a role in the decay in dielectric coefficient observed as an atomic force microscope tip gets within ≈ 10 nm of a mica surface [79].

In future studies, it may be interesting to consider the possible effects of this structured water on membrane functions, such as ion transport. In particular, structured water near the membrane surface may play a role in the interesting phenomenon of lateral diffusion by protons observed with biological membranes [8082].

These results suggest that acidic groups in cation selective channels, such as voltage-gated sodium and calcium channels, would deepen the nearby free-energy wells, but because of indirect effects the effect would be substantially muted compared to that predicted from electrostatics alone.

4. Conclusions

The presence of the negatively charged taurine residues instead of the electroneutral ethanolamine residues at the entrance of the gA channel affects both the depth and the location of the binding sites producing a strong attraction for the cations in the bulk towards the channel mouth.

The direct and indirect effects of changing the C-terminus on the test ion are interwoven in a complicated way to produce the PMF profile, such that the well depth enhancement is less significant than expected from the electrostatics of fixed charges and the translocation barrier is enhanced rather than diminished.

The high central free-energy barrier in the TgA system compared with that in the gA system might be partly a result of the spherical restraint. The deepening of the binding sites would explain the experimental observation of shieldable increase in channel conductance. If the central barrier is enhanced, it leads qualitatively to the prediction that the entry monomer, rather than the exit monomer, is responsible for the conductance increase.


We are grateful for the financial support of a Roland K. Robbins Graduate Research Fellowship. A generous allotment of computer time by the Ira and Marylou Fulton Supercomputing Center at Brigham Young University is acknowledged with thanks. We are grateful for helpful comments from Jacob Durrant, Richard Swenson and Ryan Frei. This project was supported in part by NIH grant AI23007.


Publisher's Disclaimer: This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final citable form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

Contributor Information

Morad Mustafa, Department of Chemistry and Biochemistry, Brigham Young University, Provo, Utah 84602.

Douglas J. Henderson, Department of Chemistry and Biochemistry, Brigham Young University, Provo, Utah 84602.

David D. Busath, Department of Physiology and Developmental Biology, Brigham Young University, Provo, Utah 84602, March 15, 2009.


1. Roux B, Karplus M. Ion Transport in the Gramicidin Channel: Free Energy of the Solvated Right-Handed Dimer in a Model Membrane. J Am Chem Soc. 1993;115(8):3250–3262.
2. Allen TW, Andersen OS, 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. Biophys J. 2006;90(10):34473468. [PubMed]
3. Baştuğ T, Kuyucak S. Free Energy Simulations of Single and Double Ion Occupancy in Gramicidin A. J Chem Phys. 2007;126(10):105103. [PubMed]
4. Mustafa M, Busath DD. The Gramicidin Channel Ion Permeation Free-Energy Profile: Direct and Indirect Effects of CHARMM Force Field Improvements. In press. [PMC free article] [PubMed]
5. Sansom MSP, Kerr ID. Influenza Virus M2 Protein: A Molecular Modelling Study of the Ion Channel. Protein Eng. 1993;6(1):65–74. [PubMed]
6. Wu YJ, Voth GA. A Computational Study of the Closed and Open States of the Influenza A M2 Proton Channel. Biophys J. 2005;89(4):2402–2411. [PubMed]
7. Mustafa M, Henderson DJ, Busath DD. Free-Energy Profiles for Ions in the Influenza M2-TMD Channel. Proteins: Struct, Funct, Bioinf. In press. [PMC free article] [PubMed]
8. Nonner W, Catacuzzeno L, Eisenberg B. Binding and Selectivity in L-Type Calcium Channels: A Mean Spherical Approximation. Biophys J. 2000;79(4):1976–1992. [PubMed]
9. Boda D, Valiskó M, Eisenberg B, Nonner W, Henderson D, Gillespie D. The Effect of Protein Dielectric Coefficient on the Ionic Selectivity of a Calcium Channel. J Chem Phys. 2006;125(3):034901. [PubMed]
10. Boda D, Nonner W, Valiskó M, Henderson D, Eisenberg B, Gillespie D. Steric Selectivity in Na Channels Arising from Protein Polarization and Mobile Side Chains. Biophys J. 2007;93(6):1960–1980. [PubMed]
11. Boda D, Nonner W, Henderson D, Eisenberg B, Gillespie D. Volume Exclusion in Calcium Selective Channels. Biophys J. 2008;94(9):3486–3496. [PubMed]
12. Chung SH, Corry B. Conduction Properties of KcsA Measured Using Brownian Dynamics with Flexible Carbonyl Groups in the Selectivity Filter. Biophys J. 2007;93(1):44–53. [PubMed]
13. Yang Y, Henderson D, Busath DD. Calcium Block of Sodium Current in a Model Calcium Channel: Cylindrical Atomistic Pore with Glutamate Side Chains. Mol Simul. 2004;30(23):75–80.
14. Hille B. Ion Channels of Excitable Membranes. 3. Sinauer Associates Inc; Sunderland, MA: 2001.
15. Arseniev AS, Barsukov IL, Bystrov VF. NMR Solution Structure of Gramicidin A Complex with Caesium Cations. FEBS Lett. 1985;180(1):33–39.
16. Urry DW. The Gramicidin A Transmembrane Channel: A Proposed π(L, D) Helix. Proc Natl Acad Sci USA. 1971;68(3):672–676. [PubMed]
17. Urry DW, Goodall MC, Glickson JD, Mayers DF. The Gramicidin A Transmembrane Channel: Characteristics of Head-to-Head Dimerized π(L, D) Helices. Proc Natl Acad Sci USA. 1971;68(8):1907–1911. [PubMed]
18. Sarges R, Witkop B, Gramicidin AV. The Structure of Valine- and Isoleucine-Gramicidin A. J Am Chem Soc. 1965;87(9):2011–2020. [PubMed]
19. Koeppe RE, II, Anderson OS. Engineering the Gramicidin Channel. Annu Rev Biophys Biomol Struct. 1996;25:231–258. [PubMed]
20. Kelkar DA, Chattopadhyay A. The Gramicidin Ion Channel: A Model Membrane Protein. Biochim Biophys Acta Biomembr. 2007;1768(9):2011–2025. [PubMed]
21. Roeske RW, Hrinyo-Pavlina TP, Pottorf RS, Jin BTXZ, Busath D. Synthesis and Channel Properties of [Tau16]gramicidin A. Biochim Biophys Acta. 1989;982(2):223–227. [PubMed]
22. Capone R, Blake S, Restrepo MR, Yang J, Mayer M. Designing Nanosensors Based on Charged Derivatives of Gramicidin A. J Am Chem Soc. 2007;129(31):9737–9745. [PubMed]
23. Apell HJ, Bamberg E, Alpes H, Läuger P. Formation of Ion Channels by a Negatively Charged Analog of Gramicidin A. J Mol Biol. 1977;31(1):171–188. [PubMed]
24. Bamberg E, Apell HJ, Alpes H, Gross E, Morell JL, Harbaugh JF, Janko K, Läuger P. Ion Channels Formed by Chemical Analogs of Gramicidin A. Fed Proc. 1978;37(12):2633–2638. [PubMed]
25. Cifu AS, Koeppe RE, II, Andersen OS. On the Supramolecular Organization of Gramicidin Channels. The Elementary Conducting Unit is a Dimer. Biophys J. 1992;61(1):189–203. [PubMed]
26. Woolley GA, Jaikaran ASI, Zhang Z, Peng S. Design of Regulated Ion Channels Using Measurements of Cis-Trans Isomerization in Single Molecules. J Am Chem Soc. 1995;117(16):4448–4454.
27. Heinemann SH, Terlau H, Stühmer W, Imoto K, Numa S. Calcium Channel Characteristics Conferred on the Sodium Channel by Single Mutations. Nature. 1992;356(6368):441–443. [PubMed]
28. Bahinski A, Yatani A, Mikala G, Tang SQ, Yamamoto S, Schwartz A. Charged Amino Acids Near the Pore Entrance Influence Ion-Conduction of a Human L-type Cardiac Calcium Channel. Mol Cell Biochem. 1997;166(12):125–134. [PubMed]
29. Townsley LE, Tucker WA, Sham S, Hinton JF. Structures of Gramicidins A, B, and C Incorporated into Sodium Dodecyl Sulfate Micelles. Biochemistry. 2001;40(39):11676–11686. [PubMed]
30. Jorgensen WL, Chandrasekhar J, Madura JD, Impey RW, Klein ML. Comparison of Simple Potential Functions for Simulating Liquid Water. J Chem Phys. 1983;79(2):926–935.
31. Allen TW, Andersen OS, Roux B. Structure of Gramicidin A in a Lipid Bilayer Environment Determined Using Molecular Dynamics Simulations and Solid-State NMR Data. J Am Chem Soc. 2003;125(32):9868–9877. [PubMed]
32. Humphrey W, Dalke A, Schulten K. VMD: Visual Molecular Dynamics. J Mol Graphics. 1996;14(1):33–38. [PubMed]
33. G. Black, J. Daily, B. Didier, T. Elsethagen, D. Feller, D. Gracio, M. Hackler, S. Havre, D. Jones, E. Jurrus, T. Keller, C. Lansing, S. Matsumoto, B. Palmer, M. Peterson, K. Schuchardt, E. Stephan, L. Sun, K. Swanson, H. Taylor, G. Thomas, E. Vorpagel, T. Windus, C. Winters, ECCE, A Problem Solving Environment for Computational Chemistry, Software Version 4.5.1, Pacific Northwest National Laboratory, Richland, Washington 99352-0999, USA, 2007.
34. Kendall RA, Apra E, Bernholdt DE, Bylaska EJ, Dupuis M, Fann GI, Harrison RJ, Ju JL, Nichols JA, Nieplocha J, Straatsma TP, Windus TL, Wong AT. High Performance Computational Chemistry: An Overview of NWChem a Distributed Parallel Application. Comput Phys Commun. 2000;128(1–2):260–283.
35. E. J. Bylaska, W. A. de Jong, K. Kowalski, T. P. Straatsma, M. Valiev, D. Wang, E. Aprà, T. L. Windus, S. Hirata, M. T. Hackler, Y. Zhao, P.-D. Fan, R. J. Harrison, M. Dupuis, D. M. A. Smith, J. Nieplocha, V. Tipparaju, M. Krishnan, A. A. Auer, M. Nooijen, E. Brown, G. Cisneros, G. I. Fann, H. Früchtl, J. Garza, K. Hirao, R. Kendall, J. A. Nichols, K. Tsemekhman, K. Wolinski, J. Anchell, D. Bernholdt, P. Borowski, T. Clark, D. Clerc, H. Dachsel, M. Deegan, K. Dyall, D. Elwood, E. Glendening, M. Gutowski, A. Hess, J. Jaffe, B. Johnson, J. Ju, R. Kobayashi, R. Kutteh, Z. Lin, R. Littlefield, X. Long, B. Meng, T. Nakajima, S. Niu, L. Pollack, M. Rosing, G. Sandrone, M. Stave, H. Tay-lor, G. Thomas, J. van Lenthe, A. Wong, Z. Zhang, NWChem, A Computational Chemistry Package for Parallel Computers, Version 5.0, Pacific Northwest National Laboratory, Richland, Washington 99352-0999, USA, 2006.
36. Chiu SW, Clark M, Balaji V, Subramaniam S, Scott HL, Jakobsson E. Incorporation of Surface Tension into Molecular Dynamics Simulation of an Interface: A Fluid Phase Lipid Bilayer Membrane. Biophys J. 1995;69(4):1230–1245. [PubMed]
37. Zhang Y, Feller SE, Brooks RW, Pastor RW. Computer Simulation of Liquid/Liquid Interfaces. I. Theory and Application to Octane/Water. J Chem Phys. 1995;103(23):10252–10266.
38. Feller SE, Zhang Y, Pastor RW. Computer Simulation of Liquid/Liquid Interfaces. II. Surface Tension-Area Dependence of a Bilayer and Monolayer. J Chem Phys. 1995;103(23):10267–10276.
39. Feller SE, Pastor RW. Constant Surface Tension Simulations of Lipid Bilayers: The Sensitivity of Surface Areas and Compressibilities. J Chem Phys. 1999;111(3):1281–1287.
40. Koenig BW, Strey HH, Gawrisch K. Membrane Lateral Compressibility Determined by NMR and X-ray Diffraction: Effect of Acyl Chain Polyunsaturation. Biophys J. 1997;73(4):1954–1966. [PubMed]
41. Takaoka Y, Pasenkiewicz-Gierula M, Miyagawa H, Kitamura K, Tamura Y, Kusumi A. Molecular Dynamics Generation of Nonarbitrary Membrane Models Reveals Lipid Orientational Correlations. Biophys J. 2000;79(6):3118–3138. [PubMed]
42. Högberg CJ, Lyubartsev AP. A Molecular Dynamics Investigation of the Influence of Hydration and Temperature on Structural and Dynamical Properties of a Dimyristoylphosphatidylcholine Bilayer. J Phys Chem B. 2006;110(29):14326–14336. [PubMed]
43. Petrache HI, Tristram-Nagle S, Harries D, Ku cerka N, Nagle JF, Parsegian VA. Swelling of Phospholipids by Monovalent Salt. J Lipid Res. 2006;47(2):302–309. [PMC free article] [PubMed]
44. Israelachvili JN. Intermolecular and Surface Forces. Academic Press; San Diego, CA, USA: 1992.
45. Errington JR, Kofke DA. Calculation of Surface Tension via Area Sampling. J Chem Phys. 2007;127(17):174709. [PubMed]
46. Feller SE, Pastor RW. On Simulating Lipid Bilayers with an Applied Surface Tension: Periodic Boundary Conditions and Undulations. Biophys J. 1996;71(3):13501355. [PubMed]
47. Sankararamakrishnan R, Weinstein H. Surface Tension Parameterization in Molecular Dynamics Simulations of a Phospholipid-bilayer Membrane: Calibration and Effects. J Phys Chem B. 2004;108(31):11802–11811.
48. Chiu SW, Subramaniam S, Jakobsson E. Simulation Study of a Gramicidin/Lipid Bilayer System in Excess Water and Lipid. I. Structure of the Molecular Complex. Biophys J. 1999;76(4):1929–1938. [PubMed]
49. Harroun TA, Heller WT, Weiss TM, Yang L, Huang HW. Experimental Evidence for Hydrophobic Matching and Membrane-Mediated Interactions in Lipid Bilayers Containing Gramicidin. Biophys J. 1999;76(2):937–945. [PubMed]
50. Ku cerka N, Liu Y, Chu N, Petrache HI, Tristram-Nagle S, Nagle JF. Structure of Fully Hydrated Fluid Phase DMPC and DLPC Lipid Bilayers Using X-ray Scattering from Oriented Multilamellar Arrays and from Unilamellar Vesicles. Biophys J. 2005;88(1):2626–2637. [PubMed]
51. Allen TW, Baştuğ T, Kuyucak S, Chung S-H. Gramicidin A Channel as a Test Ground for Molecular Dynamics Force Fields. Biophys J. 2003;84(4):2159–2168. [PubMed]
52. Woolf TB, Roux B. Structure, Energetics, and Dynamics of Lipid-Protein Interactions: A Molecular Dynamics Study of the Gramicidin A Channel in a DMPC Bilayer. Proteins: Struct, Funct Genet. 1996;24(1):92–114. [PubMed]
53. Phillips JC, Braun R, Wang W, Gumbart J, Tajkhorshid E, Villa E, Chipot C, Skeel RD, Kalé L, Schulten K. Scalable Molecular Dynamics with NAMD. J Comput Chem. 2005;26(16):1781–1802. [PMC free article] [PubMed]
54. MacKerell AD, Jr, Bashford D, Bellott RL, Dunbrack RL, Jr, Evanseck JD, Field MJ, Fischer S, Gao J, Guo H, Ha S, Joseph-McCarthy D, Kuchnir L, Kuczera K, Lau FTK, Mattos C, Michnick S, Ngo T, Nguyen DT, Prodhom B, Reiher WE, III, Roux B, Schlenkrich M, Smith JC, Stote R, Straub J, Watanabe M, Wiórkiewicz-Kuczera J, Yin D, Karplus M. All-Atom Empirical Potential for Molecular Modeling and Dynamics Studies of Proteins. J Phys Chem B. 1998;102(18):3586–3616. [PubMed]
55. MacKerell AD, Jr, Feig M, Brooks CL., III Extending the Treatment of Backbone Energetics in Protein Force Fields: Limitations of Gas-Phase Quantum Mechanics in Reproducing Protein Conformational Distributions in Molecular Dynamics Simulations. J Comput Chem. 2004;25(11):1400–1415. [PubMed]
56. MacKerell AD, Jr, Feig M, Brooks CL., III Improved Treatment of the Protein Backbone in Empirical Force Fields. J Am Chem Soc. 2004;126(3):698–699. [PubMed]
57. Beglov D, Benôit R. Finite Representation of an Infinite Bulk System: Solvent Boundary Potential for Computer Simulations. J Chem Phys. 1994;100(12):9050–9063.
58. Schlenkrich M, Brickmann J, MacKerell AD, Jr, Karplus M. An Empirical Potential Energy Function for Phospholipids: Criteria for Parameter Optimization and Applications. In: Merz KM, Roux B, editors. Biological Membranes: A Molecular Perspective from Computation and Experiment. Birkhäuser; Boston, MA: 1996. pp. 31–81.
59. Feller SE, MacKerell AD., Jr An Improved Empirical Potential Energy Function for Molecular Simulations of Phospholipids. J Phys Chem B. 2000;104(31):7510–7515.
60. Steinbach PJ, Brooks BR. New Spherical-Cutoff Methods for Long-Range Forces in Macromolecular Simulation. J Comput Chem. 1994;15(7):667–683.
61. Martyna GJ, Tobias DJ, Klein ML. Constant Pressure Molecular Dynamics Algorithms. J Chem Phys. 1994;101(5):4177–4189.
62. Feller SE, Zhang Y, Pastor RW, Brooks BR. Constant Pressure Molecular Dynamics Simulation: The Langevin Piston Method. J Chem Phys. 1995;103(11):4613–4621.
63. Welch BB, Jones K, Hobbs J. Practical Programming in Tcl and Tk. 4. Prentice Hall PTR; Upper Saddle River, NJ: 2003.
64. Roux B. Statistical Mechanical Equilibrium Theory of Selective Ion Channels. Biophys J. 1999;77(1):139–153. [PubMed]
65. Allen TW, Andersen OS, Roux B. Energetics of Ion Conduction through the Gramicidin Channel. Proc Natl Acad Sci USA. 2004;101(1):117–122. [PubMed]
66. Allen TW, Andersen OS, Roux B. Molecular Dynamics — Potential of Mean Force Calculations as a Tool for Understanding Ion Permeation and Selectivity in Narrow Channels. Biophys Chem. 2006;124(3):251–267. [PubMed]
67. Grossfield A. An Implementation of WHAM: The Weighted Histogram Analysis Method. available from:, ????
68. Souaille M, Roux B. Extension to the Weighted Histogram Analysis Method: Combining Umbrella Sampling with Free Energy Calculations. Comput Phys Commun. 2001;135(1):40–57.
69. Roux B, Karplus M. Ion Transport in a Model Gramicidin Channel. Structure and Thermodynamics. Biophys J. 1991;59(5):961–981. [PubMed]
70. Trzesniak D, Kunz APE, van Gunsteren WF. A Comparison of Methods to Compute the Potential of Mean Force. ChemPhysChem. 2007;8(1):162–169. [PubMed]
71. Burden RL, Faires JD. Numerical Analysis. 7. Brooks Cole; Pacific Grove, CA: 2001.
72. Petrache HI, Dodd SW, Brown MF. Area per Lipid and Acyl Length Distributions in Fluid Phosphatidylcholines Determined by 2H NMR Spectroscopy. Biophys J. 2000;79(6):3172–3192. [PubMed]
73. Otten D, Brown MF, Beyer K. Softening of Membrane Bilayers by Detergents Elucidated by Deuterium NMR Spectroscopy. J Phys Chem B. 2000;104(51):12119–12129.
74. Aussenac F, Laguerre M, Schmitter JM, Dufourc EJ. Detailed Structure and Dynamics of Bicelle Phospholipids Using Selectively Deuterated and Perdeuterated Labels. 2H NMR and Molecular Mechanics Study. Langmuir. 2003;19(25):10468–10479.
75. Baştuğ T, Kuyucak S. Energetics of Ion Permeation, Rejection, Binding, and Block in Gramicidin A from Free Energy Simulations. Biophys J. 2006;90(11):3941–3950. [PubMed]
76. Baştuğ T, Kuyucak S. Test of Molecular Dynamics Force Fields in Gramicidin A. Eur Biophys J. 2005;34(5):377–382. [PubMed]
77. Baştuğ T, Patra SM, Kuyucak S. Molecular Dynamics Simulations of Gramicidin A in a Lipid Bilayer: From Structure-Function Relations to Force Fields. Chem Phys Lipids. 2006;141(12):197–204. [PubMed]
78. Crozier PS, Henderson D, Rowley RL, Busath DD. Model Channel Ion Currents in NaCl-Extended Simple Point Charge Water Solution with Applied-Field Molecular Dynamics. Biophys J. 2001;81(6):3077–3089. [PubMed]
79. Teschke O, Ceotto G, de Souza EF. Interfacial Water Dielectric-Permittivity-Profile Measurements Using Atomic Force Microscopy. Phys Rev E. 2001;64(1):011605. [PubMed]
80. Kell DB. On the Functional Proton Current Pathway of Electron Transport Phosphorylation: An Electrodic View. Biochim Biophys Acta. 1979;549(1):55–99. [PubMed]
81. Cherepanov DA, Feniouk BA, Junge W, Mulkidjanian AY. Low Dielectric Permittivity of Water at the Membrane Interface: Effect on the Energy Coupling Mechanism in Biological Membranes. Biophys J. 2003;85(2):1307–1316. [PubMed]
82. Antonenko YN, Pohl P. Microinjection in Combination with Microfluorimetry to Study Proton Diffusion along Phospholipid Membranes. Eur Biophys J. 2008;37(6):865–870. [PubMed]