Multinanosecond molecular dynamics simulations of gramicidin A embedded in a dimyristoylphosphatidylcholine bilayer show a remarkable structural stability for both experimentally determined conformations: the head-to-head helical dimer and the double helix. Water permeability was found to be much higher in the double helical conformation, which is explained by lower hydrogen bond-mediated enthalpic barriers at the channel entrance and its larger pore size. Free-energy perturbation calculations show that the double helical structure is stabilized by the positive charges at the N termini introduced by the desformylation, whereas the helical dimer is destabilized. Together with the recent experimental observation that desformyl gramicidin conducts water hundredfold better than gramicidin, this suggests that desformyl gramicidin A predominantly occurs in the double helical conformation.

In this study, we test this hypothesis by extended MD simulations of the HD and DH conformations of gA, both in the native formyl form and in the desformylated form. Additionally, the free-energy changes for the desformylation are estimated from perturbation calculations for both the HD and DH conformation.

It is well known that with the transport of ions through the gA channel, a number of water molecules are co-translated across the membrane (). In a detailed molecular dynamics (MD) study (), the diffusion of water molecules through the HD conformation of the peptide has been characterized (), and the water permeability was found to be in reasonable agreement with experimental values (). Surprisingly, a hundredfold higher water permeability was recently observed for desformylated gA (). It has been speculated () that the HD conformation is destabilized by the positive charges at the N termini caused by the desformylation, because they would be in close proximity in the hydrophobic membrane center. As mechanism for the higher water permeability, it was suggested that the positive charges at the channel entrance of the DH conformation would lower the access resistance to water and thereby enhance water permeability.

More recently, the notion emerged that with environmental conditions like the length of the alkyl chains of the bilayer in which the peptide is embedded, the equilibrium can be shifted from one of the multiple conformations to the other ().

The peptide antibiotic gramicidin A (gA) is one of the best-characterized ion channels (for review, see e.g.,). It is hydrophobic and specific for monovalent cations (). Being a small peptide forming well-behaved channels, it has served as a prototypic model for more complicated ion channels in numerous experimental () and theoretical () studies. The primary structure consists of 15 alternating- and-amino acids, and both the N and C termini are capped by a formyl and ethanolamine group, respectively, such that the complete peptide is uncharged (). Despite extensive efforts to elucidate the active channel structure, an apparent controversy still exists concerning this issue (). Two major structural types are known, a head-to-head helical dimer (HD, Fig. 1 , left) () and a double-helical structure (DH, Fig. 1 , right) (), both of which have been claimed to be the major ion-conducting conformation.

Osmotic permeation coefficients (p, in cm) were derived from the equilibrium simulations as follows. For permeation driven by osmotic pressure gradients, pis defined as −Φ/Δn (), in which Φ is the net flux of water molecules across the pore in mol/s, and Δn is the solute concentration difference (mol/cm). Using the Van’t Hoff equation (Π = Δn kT, with Π the osmotic pressure gradient, kBoltzmann’s constant, and T the temperature), this yields Φ = −pΠ/kT. Using rate theory, we obtain: Φ = Φsinh(½ ΔG/kT). For small osmotic pressure differences, (i.e., ΔG ≪ kT), Φ ≈ ½ ΦΔG/kT. Φis the intrinsic flux, defined by the barrier height, and can be derived from the simulations. ΔG = ΠV is the free energy required to transfer one water molecule (molar volume v) from one side of the membrane to the other side against the osmotic pressure gradient. Combining the two expressions for Φ, allows to derive the permeation coefficient from the simulations, p≈ ½ Φv.

Hydrogen bond energies were calculated using the empirical formula for the potential hydrogen bond energy as described in. Hydrogen bond energies per water molecule were evaluated by binning the oxygen position of each water molecule in each snapshot along the pore axis in slices with a width of 0.25 Å and subsequent averaging over each bin. Water oxygen positions were evaluated relative to the peptide position by performing a least-squares fit onto the gA starting structure before the binning. Pore dimensions were evaluated using the HOLE program (). Free energy profiles G(z) for water molecules passing the gA pores were computed from the residence occupancy n(z) of the water molecules along the pore axis z as G = −kT ln n(z), evaluated in 0.1-Å bins.

All MD simulations were performed with the Gromacs () simulation software. Partial charges for the formyl group were derived from the Charmm force field (). DMPC force field parameters were taken from

In both the DH structure as well as in the head-to-head HD structure, the free energy of the transition into desformyl gA was calculated by slow-growth free energy perturbation (FEP) simulations (). The carbon atom of the formyl group was perturbed into a proton, and simultaneously a second proton was created at the terminal nitrogen. At the same time, the oxygen and hydrogen atom of the formyl group were mutated into dummy atoms. Additionally, two simulations were conducted in which all peptide atoms were slowly mutated into dummy atoms without interactions with the surrounding lipid/solvent environment. The charges created in these simulations were not compensated. Although this may affect the calculated free energies, this effect can be expected to be similar for both desformylation reactions, and therefore does not qualitatively alter the thermodynamic cycle. All FEP simulations were performed using a soft-core potential with an α value of 1.0 (), a time step of 1 fs, and had a total length of 1 ns. For the FEP simulations, a charge-group based twin-range electrostatic cutoff radius of 1.0/1.8 nm (and 1.0/1.5 nm for comparison) were used. The Particle-Mesh Ewald contribution to the perturbed interactions cannot be evaluated with the present version of the gromacs simulation software ().

In all MD simulations the LINCS () and SETTLE () algorithms were used to constrain covalent bond lengths, allowing a time step of 2 fs. The temperature was kept constant at 300 K by separately coupling gA, lipids, and solvent (time constant τ = 0.1 ps) to an external temperature bath (). At this temperature the DMPC lipids are in the liquid crystalline phase. The pressure was controlled by anisotropic coupling (τ = 1 ps) to a pressure bath of 1 bar during the equilibration phase. In the production runs the volume was only allowed to change in the Z direction (perpendicular to the membrane). Electrostatic interactions between charge groups at a distance less than 1 nm were calculated explicitly, long-range electrostatic interactions were calculated using the Particle-Mesh Ewald method () with a grid width of 0.12 nm and a fourth-order spline interpolation. A cutoff distance of 1.0 nm was aplied for Lennard-Jones interactions.

The four structures so obtained were energy-minimized and subsequently equilibrated for 200 ps, applying positional restraints to the nonhydrogen atoms of gA, to allow lipid and water molecules to relax around the peptide. The resulting four structures were used as starting structures for production simulations of 10 ns each.

The DMPC starting structure was created by scaling an existing dipalmitoylphosphatidylcholine structure (coordinates after 500 ps of simulation E fromto an area per lipid of 0.596 nm), removing the last two carbon atoms of each chain, translating one leaflet 4 Å toward the other leaflet, energy minimizing, and equilibrating for 500 ps at constant area). To insert the peptides into the bilayer with as little perturbation of the bilayer as possible, we used the procedure of: first, two lipids from each bilayer leaflet were removed. Then a molecular surface representation of DH was used to define a volume from which lipids were expelled by a repulsive force normal to the surface. The peptide DH was then inserted into the resulting hole, the system minimized, and a MD simulation of 50 ps with position restraints on the nonhydrogen atoms of the peptide was carried out. This procedure was repeated for the HD conformation. The structures of desformyl gA in the HD and DH conformation were modeled based on the corresponding gA starting conformations (after equilibration) by replacing the two formyl groups by two protons each.

As starting structures for the MD simulations the experimental structures with protein database entries 1MAG () and 1AV2 () were taken for the HD and the DH conformation of gA, respectively. Both structures were embedded in a preequilibrated dimyristoylphosphatidylcholine (DMPC) bilayer of 124 lipids, solvated at both sides with 3659 simple point change water molecules (), adding up to a simulation system of 17,037 atoms. The periodic simulation box is depicted in Fig. 2

Simulation system setup, side view (a) and top view (b). All molecular dynamics simulations were carried out with full electrostatics in a periodic simulation box containing the peptide (blue) embedded within a DMPC lipid bilayer (yellow head-groups and green tails) surrounded by water (red/white). The total system consists of 17,037 atoms.

Figure 2 Simulation system setup, side view (a) and top view (b). All molecular dynamics simulations were carried out with full electrostatics in a periodic simulation box containing the peptide (blue) embedded within a DMPC lipid bilayer (yellow head-groups and green tails) surrounded by water (red/white). The total system consists of 17,037 atoms.

Results and discussion

Figure 3 RMSD of the four studied conformations/forms from their starting structures as a function of time, calculated over all backbone atoms from the dimeric gramicidin A molecules for the four simulations that were carried out. The root mean square deviation (RMSD) from the respective experimental starting structures in each of the four 10-ns simulations is depicted in Fig. 3 . In both the native and the desformylated form (red and blue curves, respectively), the DH conformation is very stable (average dimer backbone RMSD of 0.6 Å) throughout the entire simulations. Also the native, formylated, HD conformation (black curve) is stable during the simulation (with an average RMSD of 1.1 Å and a final RMSD of 1.3 Å after 10 ns). The desformylated HD conformation (green curve), however, shows a large drift and reaches an RMSD value of ∼3 Å after 10 ns. This drift is presumably due to the introduction of the two positive charges at the N termini of the peptide that meet in the hydrophobic membrane center. The electrostatic repulsion between these positive charges cause the N termini to fold toward the interior the channel region, thereby destabilizing the helical dimeric structure and interrupting the water chain through the pore.

Figure 4 Water distribution in the two simulated conformations of gA; the helical dimer (left) and the double-helical conformation (right). Fig. 4 shows snapshots of the distribution of water molecules in both the HD and DH conformations of the native gA channel after equilibration. A single-file hydrogen-bonded water column is observed in both conformations, with seven to nine water molecules present inside each of the two conformations.

Pómès and Roux, 1996 Pómès R.

Roux B. Structure and dynamics of a proton wire: a theoretical study of H+ translocation along the single-file water chain in the gramicidin A channel. Chiu et al., 1999b Chiu S.-W.

Subramaniam S.

Jakobsson E. Simulation study of the gramicidin/lipid bilayer system in excess water and lipid: II. Rates and mechanism of water transport. Figure 5 Spontaneous motion of water molecules in the pore during the simulations of the HD (top panel) and DH (bottom panel) gA conformation. Trajectories along the pore axis of all water molecules that visit the pore region are shown in Fig. 5 . As can be seen, there is a dramatic difference in water mobility between the HD and DH conformations. In the HD form, none of the water molecules that were initially placed inside the pore leave the pore, and all water molecules fluctuate around a well-defined position inside the gA matrix, in agreement with other simulation studies (). In contrast, in the native DH conformation, the water molecules show collective motion along the pore direction, and in total 17 spontaneous passages were observed in the simulation time of 10 ns. In both the HD and DH conformations, the water molecules fluctuate in a highly correlated manner, keeping the hydrogen-bonded water chain intact at all times. Principal component analysis of these fluctuations shows that this correlated fluctuation of water molecules contributes 70% to 80% to the total fluctuation in the pore direction (data not shown). This implies that the relevant event that determines the water permeation rate is the hopping of the whole water column by one water-water distance, rather than the hopping of single water molecules. Such events were encountered 3, 62, and 58 times for the native HD, the native DH, and the desformyl DH conformation, respectively, during the 10-ns simulations. For the HD conformation, this number is too small to reliably compute a rate, but the ratio between the numbers shows that the water permeation rates for the HD and DH conformations differ at least by an order of magnitude.

Chiu et al., 1999b Chiu S.-W.

Subramaniam S.

Jakobsson E. Simulation study of the gramicidin/lipid bilayer system in excess water and lipid: II. Rates and mechanism of water transport. The peptide dimer is significantly shorter (∼25 Å) than the DMPC bilayer (∼40 Å) for both the HD and DH conformations ( Fig. 2 ). As has been observed before (), also in our simulations lipid head groups interfere with water exchange at the channel entrance. In both the simulation of the native HD and DH conformations, a lipid head-group was found at each side of the channel entry during most of the simulation length, thereby significantly affecting water permeability. In the HD conformation the lipid headgroups remain more tightly associated to the peptide compared with the DH conformation. This effect causes the barrier for water molecules to cross the channel entry/exit to be higher in the HD conformation than in the DH conformation.

For both simulations of the DH conformation, the number of partial permeation events corresponds to a water permeation coefficient of ∼1.0 × 10−13 cm3 s−1, which is some- what smaller than the experimentally determined permeation coefficient for desformyl gA of 1.1 × 10−12 cm3 s−1. The lower computational rate is presumably caused by the behavior of the lipid headgroups at the channel entries/exits: occasional diffusion of these headgroups away from the channel entry, on a time scale beyond the 10 ns that was simulated, would result in higher water transmisson rates.

Chiu et al., 1999b Chiu S.-W.

Subramaniam S.

Jakobsson E. Simulation study of the gramicidin/lipid bilayer system in excess water and lipid: II. Rates and mechanism of water transport. Figure 6 Hydrogen bond energies per water molecule along the pore axis, averaged over all water molecules that visited the pore regions at some instance of time during the simulations, for the HD conformation (top panel), the DH conformation (middle panel), and the desformylated DH conformation (bottom panel). Schematic drawings of the HD and DH structures are included at the correct scale to illustrate the location of the barriers. Fig. 6 shows average hydrogen bond energies of individual water molecules as a function of the position along the pore axis. For the HD simulation (top panel), the main enthalpic barriers for water molecules passing the channel are located at the channel entrance. In this region both the number of water-water hydrogen bonds per water molecule as well as the strength per hydrogen bond is reduced, to approximately one-fourth of their strength in bulk water. This loss can only partially be compensated by gA-water hydrogen bonds, such that a total enthalpy barrier of ∼20 kJ/mol still remains. The change of hydration environment for water molecules upon entering the channel has been suggested to be the major barrier (). This finding is confirmed in our simulations. Additionally, the barriers may in part be formed by lipid head groups, which have only limited hydrogen-bonding capabilities, that are located near the channel entry and thereby interfere with water entrance and exit from the channel. In contrast, no significant hydrogen-bond energetic barriers are observed in this case ( Fig. 6 , middle panel), despite the fact that also here, lipid head-groups are also found near the channel entrance during the simulation of the DH conformation. In the center of the channel, both gA-water and water-water hydrogen bonds are on average stronger in the HD conformation than in the DH conformation. This also lowers the total hydrogen-bond energy per water molecule in the central part of the HD channel with respect to DH. The hydrogen-bond enthalpic barrier for the DH conformation is ∼15 kJ/mol.

No significant differences are observed between the native and desformylated DH conformation ( Fig. 6 , middle and bottom panel). The positively charged amino group at the N termini in the desformylated conformation forms hydrogen bonds of similar strength compared with the native, formylated form. Results for the simulation of the desformylated HD conformation are not shown, since that simulation shows a large structural drift, and therefore, averaging over multiple snapshots as was done for the other trajectories is not meaningful here.

Hao et al., 1997 Hao Y.

Pear M.R.

Busath D.D. Molecular dynamics study of free energy profiles for organic cations in gramicidin A channels. Figure 7 Free energy profiles of water molecules along the pore axis (solid lines) and average pore radii (dashed lines), for the simulation of the HD conformation (top panel), the DH conformation (middle panel) and the desformylated DH conformation (bottom panel). Note that the free energy along the pores also depends on the gA concentration within the membrane. Therefore, to enable comparison, the ground level is chosen equal for all three panels, assuming equal concentration. No significant barriers were observed outside the shown regions. Free energy profiles were calculated for water molecules passing the pore from the occupancy distribution along the pore axis (see Materials and Methods). The well-defined positions of the seven water molecules that remain bound to the HD conformation throughout the simulation cause the distinct oscillating free energy pattern with seven minima seen in Fig. 7 (top panel). No comparable profile is seen for the hydrogen-bond energies ( Fig. 6 , top), which are the main interactions for water molecules with the gA channel, for the central part of the channel. Therefore, the positions of these water molecules in the pore (and hence the corresponding minima and maxima in the free energy profile) are not so much defined by specific interactions with the gA molecule but are predominantly determined by the free energy barriers located at the channel entry and exit. These keep the water molecules in the channel in a “locked” position, at a distance of ∼2.6 Å from each other, i.e., in the optimal hydrogen bond distance. The free energy barriers at the channel entry and exit are mainly determined by the hydrogen-bond energy profile ( Fig. 6 ) and have also been found to be the main barrier for organic cations passing the pore ().

Elber et al., 1995 Elber R.

Chen D.P.

Rojewska D.

Eisenberg R. Sodium in gramicidin: an example of a permion. Note that the distinct pattern of peaks in the free energy profile is derived from the single-particle occupancy distribution of water molecules along the pore axis. As described above, however, the correlated fluctuation of water molecules in the pore region, renders the relevant coordinate for describing water permeation a collective one, mainly composed of a column of several water molecules, analogous to the “permion” concept (). The relevant free energy barrier for this coordinate is the energetic cost of the hopping of this water column (in our case defined as a vertical column of six water molecules) by one water-water distance, such that each water molecule shifts into its neighbor’s position. The free energy profiles show that the height of this barrier is determined predominantly by the channel entry and exit and not so much by the (relatively smooth) channel interior. A closer analysis of the collective coordinate consisting of a column of six water molecules showed that the barrier for of shifting the water column by one water-water distance across this barrier is ∼10 to 12 kJ/mol.

−14 cm3 s−1 and 1.1 × 10−12 cm3 s−1 for gA and desformyl gA, respectively ( Pohl and Saparov, 2000 Pohl P.

Saparov S.M. Solvent drag across gramicidin channels demonstrated by microelectrodes. Saparov et al., 2000 Saparov S.M.

Antonenko Y.N.

Koeppe R.E.

I I.

Pohl P. Desformylgramicidin: a model channel with an extremely high water permeability. Compared with the HD conformation, a much smoother free energy profile was obtained from the simulation of the DH conformation ( Fig. 7 , middle panel). The associated free energy barrier is estimated to be in the order of 2.5 kJ/mol at room temperature. The difference in the calculated barrier heights is in good agreement with the observed difference in permeation rates for both conformations (1.6 × 10cmand 1.1 × 10cmfor gA and desformyl gA, respectively ()). Note that the difference between the (collective) free energy barrier heights for HD and DH is much larger than the (single molecule) hydrogen-bond enthalpic barriers due to the highly correlated motion of water molecules. Not only the hydrogen bond energy profile, but also the free energy profile for the desformylated DH conformation is similar to the native, formylated form ( Fig. 7 , bottom panel), as is the observed permeation rate.

In addition to the hydrogen-bond energy term, sterical effects may also play a role in the lower observed permeation rate for the HD conformation as compared with the DH conformation. In particular, the average pore radius of the DH conformation is ∼0.2 Å wider than that of the HD conformation (1.6 and 1.4 Å, respectively, see Fig. 7 , dashed lines). Also the minimal pore radius is ∼0.2 Å wider in the DH conformation (1.5 and 1.3 Å). As was also the case for the hydrogen-bond energies, the desformylated DH conformation is similar to the native, formylated form (bottom panel). Apart from a direct (van der Waals) energetic barrier that is accompanied by this steric effect, also an entropic effect can be expected for passing water molecules: the wider pore as observed for DH allows passing water molecules more rotational and translational freedom, thereby lowering the entropic barrier with respect to the HD conformation.

Saparov et al., 2000 Saparov S.M.

Antonenko Y.N.

Koeppe R.E.

I I.

Pohl P. Desformylgramicidin: a model channel with an extremely high water permeability. Taken together, the simulations show that the DH conformation is at least 10 times more permeable to water than the HD conformation. The insertion of a positive charge at the N terminus by desformylation hardly affects the behavior of water molecules in the DH conformation. Therefore, the recent hypothesis () that the positive charge at the channel entry might be involved in lowering the access resistance for water to enter the channel and thereby increasing water permeability is not confirmed. Instead, a shift in equilibrium between the HD and DH conformation toward DH is proposed here to cause the higher observed water permeation rate in desformyl gA.

Figure 8 Thermodynamic cycle of the calculated free energies. The vertical, solid arrows depict the simulations that were actually performed, with the corresponding numbers the associated free energy change (in kilojoule/mole) including long-range electrostatics (the numbers in parenthesis denote the free energy change including only short-range electrostatics, until 1.5 nm). The desformylation was achieved by mutating the native formyl-capped N termini into positively charged amino groups (see text). The horizontal, dashed arrows denote the transitions of interest, and estimates for the associated free energies are shown. “0” denotes a state in which the peptide is completely vanished. This is achieved by slowly mutating all gA atoms into dummy atoms, without interactions with the surrounding bilayer and solvent. To investigate the effect of the desformylation on the relative stability of the channel in more detail, free energy perturbation simulations were carried out. In these simulations, the formyl group was slowly “mutated” into a positively charged amino group for both the HD and DH conformation (for details, see Materials and Methods). The results are summarized in Fig. 8 . As already observed from the RMSD values ( Fig. 3 ), desformylation leads to a structural destabilization of the HD conformation. Energetically, however, introduction of the positive charges at the N termini of the peptide is favorable, even (at first sight unexpectedly) for the HD conformation. When only short-range electrostatic interactions (up to 15 Å) are taken into account, the energy required for the mutation (values in brackets) is positive, due to the two positive charges in close proximity. Therefore, the favorable decrease of this energy is a long-range electrostatic effect (presumably caused by interactions with dipolar lipid head-groups and water molecules). Note that this free energy change corresponds largely to the transfer of two charges from vacuum to the lipid environment. In contrast, the transfer of an ion from bulk water to the membrane, as would occur in forming the desformyl gA HD conformation, is a highly unfavorable process. Further note that the end structure of the simulation is structurally unstable (see also Fig. 3 ) due to short-range electrostatic repulsion. The conformational response to the perturbation continues at the end of the FEP simulation, rendering the desformyl HD conformation (and the corresponding free energy) only a snapshot along the pathway toward the equilibrium desformyl gA structure.

For the DH conformation, in contrast, also the short-range contribution to the free energy of desformylation is negative, rendering the total energy for the mutation more favorable than in the HD conformation. This can be explained from the fact that in the DH conformation the two positive charges are located on both sides of the membrane and thus can be fully solvated by water molecules and polar lipid head-groups.

Unfortunately, the conformational transition between the HD and DH conformation is too complex and, therefore, cannot be simulated directly. This renders a direct comparison of the stability of the HD and the DH conformation problematic, because the free energy difference between the native HD and DH structures, which are required to close the thermodynamic cycle, cannot be directly computed. To overcome this problem, we defined a nonphysical “empty” state (denoted “0”) with all peptide atoms mutated to dummy atoms, i.e., without interactions with the surrounding lipids or with water. Because this state is identical for both structural forms, the transitions to this state connect the left- and right-hand side of the thermodynamic circle and hence the “0” state is a suitable reference structure for calculating the free energy differences of interest; those between the HD and DH conformation in both the native and desformyl form. Note that the “mutation” of the complete peptide into dummy atoms is still a rather drastic perturbation, which cannot be expected to be fully equilibrated in the nanosecond simulation time. Yet, full closure of the membrane is observed in the simulations, which justifies to interpret the free energies associated with the transitions to the empty “0” state at least qualitatively.

The calculated free energy change to the empty “0” state is almost identical for the two structural forms. Because the left and right half of the thermodynamic cycle are now connected, one can conclude that the conformational transition from the HD to the DH conformation for desformyl gA must be associated with a highly favorable free energy change. Hence, desformylation shifts the natural equilibrium between the two structural forms toward a higher population of the DH conformation as compared with native gA.

Interestingly, the very similar free energy values obtained for the simulations to the “0” state for both conformations imply that the HD and DH conformation are about equally stable in the simulation setup. Therefore, the HD and DH dimer can be expected to coexist in DMPC bilayers. Note, however, that from these simulations no accurate assessment of the population of both conformations can be obtained.

Saparov et al., 2000 Saparov S.M.

Antonenko Y.N.

Koeppe R.E.

I I.

Pohl P. Desformylgramicidin: a model channel with an extremely high water permeability. The experimentally observed lower cation conductance for desformyl gA as compared with gA () can be rationalized by the electrostatic repulsion introduced by the positively charged channel entry and exit in desformyl gA. This repulsion can be expected to create an additional free energy barrier for ions passing the pore. Explicit simulations of cation permeation, accompanied by free energy calculations would certainly provide more detailed information on this phenomenon.