Molecular Mechanism of H Conduction in the Single-File Water Chain of
the Gramicidin Channel
Re´gis Pome`s* and Benoıˆt Roux†
*Structural Biology and Biochemistry, Hospital for Sick Children, and Department of Biochemistry, University of Toronto, Toronto, Ontario
M5G 1X8, Canada; and †Biochemistry Department, Weill Medical College of Cornell University, New York, New York 10021 USA
ABSTRACT The conduction of protons in the hydrogen-bonded chain of water molecules (or “proton wire”) embedded in
the lumen of gramicidin A is studied with molecular dynamics free energy simulations. The process may be described as a
“hop-and-turn” or Grotthuss mechanism involving the chemical exchange (hop) of hydrogen nuclei between hydrogen-
bonded water molecules arranged in single file in the lumen of the pore, and the subsequent reorganization (turn) of the
hydrogen-bonded network. Accordingly, the conduction cycle is modeled by two complementary steps corresponding
respectively to the translocation 1) of an ionic defect (H) and 2) of a bonding defect along the hydrogen-bonded chain of
water molecules in the pore interior. The molecular mechanism and the potential of mean force are analyzed for each of these
two translocation steps. It is found that the mobility of protons in gramicidin A is essentially determined by the fine structure
and the dynamic fluctuations of the hydrogen-bonded network. The translocation of H is mediated by spontaneous (thermal)
fluctuations in the relative positions of oxygen atoms in the wire. In this diffusive mechanism, a shallow free-energy well
slightly favors the presence of the excess proton near the middle of the channel. In the absence of H, the water chain adopts
either one of two polarized configurations, each of which corresponds to an oriented donor-acceptor hydrogen-bond pattern
along the channel axis. Interconversion between these two conformations is an activated process that occurs through the
sequential and directional reorientation of water molecules of the wire. The effect of hydrogen-bonding interactions between
channel and water on proton translocation is analyzed from a comparison to the results obtained previously in a study of
model nonpolar channels, in which such interactions were missing. Hydrogen-bond donation from water to the backbone
carbonyl oxygen atoms lining the pore interior has a dual effect: it provides a coordination of water molecules well suited both
to proton hydration and to high proton mobility, and it facilitates the slower reorientation or turn step of the Grotthuss
mechanism by stabilizing intermediate configurations of the hydrogen-bonded network in which water molecules are in the
process of flipping between their two preferred, polarized states. This mechanism offers a detailed molecular model for the
rapid transport of protons in channels, in energy-transducing membrane proteins, and in enzymes.
INTRODUCTION
The control of proton fluxes across biomembranes consti-
tutes one of the most complex and fundamental properties
of living systems. To achieve biological energy conversion
(bioenergetics), membrane-spanning protein assemblies use
energy provided by photochemical or redox reactions to
pump protons against an electrochemical gradient, H.
In turn, the flow of H in the direction of the electrochem-
ical gradient drives ATP synthesis (Saraste, 1999). This
“chemiosmotic coupling” is the cornerstone of bioenerget-
ics (Mitchell, 1961). Despite the importance of proton trans-
port phenomena, detailed molecular mechanisms have re-
mained elusive. A high level of detail is required to
understand proton-pumping mechanisms (Wikstro¨m, 1998;
Sjogren et al., 2000; Lanyi, 2000). In general, it is particu-
larly challenging to understand the forces driving proton
movement in enzymes, because it requires knowledge, at
the atomic level, of three equilibrium properties that are
intimately coupled to each other: 1) the protonation state of
all titratable groups in the protein; 2) the electronic (charge)
state of the protein; and 3) the equilibrium distribution of
conformational states of the enzyme. Furthermore, these
properties must be characterized at each stage of the cata-
lytic cycle. This is a formidable undertaking for the complex
membrane-bound protein assemblies involved in energy
transduction and/or proton pumping, as illustrated in the
cases of bacteriorhodospin, a light-driven proton pump (La-
nyi, 1999), of bacterial photosynthetic reaction centers
(Okamura et al., 2000), and of cytochrome c oxidase, a
redox-coupled proton pump (Zaslavsky and Gennis, 2000).
In addition to these equilibrium properties, the structure
and function of well-defined pathways for proton translo-
cation, without which leaks resulting in the loss of proton
activity would occur, must be elucidated. The transient
events (nonequilibrium properties) involved in proton
movement present a special challenge in their own right.
Unlike that of other ions, the transport of protons does not
require the net diffusion of atomic or molecular species, but
may instead take place according to a Grotthuss relay mech-
anism involving the chemical exchange of hydrogen nuclei
along successive hydrogen-bond donor and acceptor groups
forming extensive networks and the subsequent reorganiza-
tion of these networks (Grotthuss, 1806; Nagle and Moro-
witz, 1978; Knapp et al., 1980; Agmon, 1995) (Fig. 1). Such
Submitted August 14, 2001, and accepted for publication January 22,
2002.
Address reprint requests to Dr. Re´gis Pome`s, Structural Biology & Bio-
chemistry, Hospital for Sick Children, 555 University Avenue, Toronto,
Ontario M5G 1X8, Canada. Tel.: 416-813-5686; Fax: 416-813-5022;
E-mail: pomes@sickkids.ca.
© 2002 by the Biophysical Society
0006-3495/02/05/2304/13 $2.00
2304 Biophysical Journal Volume 82 May 2002 2304–2316

pathways constitute “proton wires” mediating H displace-
ment over long distances.
Although the molecular structures of such important pro-
ton pumps as bacteriorhodopsin and cytochrome c oxidase
are known, the detailed nature of proton pathways and the
molecular mechanisms leading to proton translocation are
still unclear. In addition to proton-relaying amino acid res-
idues, energy-transducing proteins, like ion channels, ap-
pear to use water wires. Models for the relay of H by
buried water molecules have been substantiated in several
systems of bioenergetic interest. In particular, there is now
abundant evidence for the involvement of buried water
molecules in bacteriorhodopsin (for recent reviews, see
Dencher et al., 2000; Luecke, 2000; Kandori, 2000). How-
ever, high-resolution structures of intermediates in the
pumping cycle suggest that although several water mole-
cules reside in the protein interior, they may not at all times
form a continuous hydrogen-bonded chain (Lanyi, 2000).
Thus, water wires make up extended but incomplete tracts
of the proposed H pathways, underlining the importance
of conformational changes and dynamic fluctuations in pro-
ton transport. This, in turn, raises further questions: do water
chains function as passive units, or are they involved in the
controlled access, blockage, and gating of proton flow?
Understanding the properties governing proton transport in
hydrogen-bonded networks containing water molecules is a
prerequisite to achieving full description of both ion perme-
ation and energy-transduction mechanisms. To this end, the
study of simple proton wires can help in characterizing the
fundamental principles governing proton translocation.
Gramicidin constitutes a model of choice for the study of
proton conduction in much more complex proteins (Quigley
et al., 1999; Cukierman, 2000; DeCoursey and Cherny,
1999). With the notable exception provided by the potas-
sium channel KcsA (Doyle et al., 1998), it is to this day the
only ion channel for which detailed structure-function rela-
tionships have been characterized, both experimentally
(Tian and Cross, 1999) and theoretically (Roux and Kar-
plus, 1994). The relative structural and functional simplicity
of gramicidin A (gA) permits one to approach the proton
transport mechanism “in isolation.” In its active form, gA
assembles as a head-to-head homodimer of pentadecapep-
tides in lipid bilayers (Arseniev et al., 1985). Its alternating
L- and D-amino acids adopt a right-handed 6.3-helix struc-
ture, which leaves its hydrophobic side chains facing out
into the bilayer and its peptide backbone lining the interior
of a cylindrical pore 4 Å in diameter (Fig. 2). This hydro-
philic pore accommodates a single file of water molecules
and mediates the translocation of monovalent cations such
as H, Cs, Na, and K (Finkelstein, 1987). The partial
dehydration of cations upon entry into the single-file region
is partly compensated by the channel backbone. Whereas
permeation of other ions necessitates the net diffusion of the
single-file water column, the absence of streaming poten-
tials during H permeation through gramicidin (Levitt et al.,
1978) indicates that the long-range translocation of H
arises from a Grotthuss mechanism (Akeson and Deamer,
1991).
In recent years, theoretical studies of proton transport by
hydrogen-bonded networks of water molecules have opened
the way to the study of biological proton wires. The molec-
ular basis for the hop-and-turn mechanism has been inves-
tigated in bulk water (Tuckerman et al., 1995; Vuillemier
and Borgis, 1998, 1999; Schmitt and Voth, 1999), in single-
file water chains or “water wires” embedded in model
channels (Pome`s and Roux, 1995, 1996a, 1998; Pome`s,
1999; Drukker et al., 1998; Decornez et al., 1999; Mei et al.,
1998; Sadeghi and Cheng, 1999; Brewer et al., 2001), and in
gA (Pome`s and Roux, 1996b; Sagnella et al., 1996). These
studies have uncovered important aspects of the hop and
turn mechanism at the molecular level. In particular, the
respective roles of thermal fluctuations and of nuclear quan-
FIGURE 1 Schematic depiction of the hop-and-turn or Grotthuss mech-
anism for H transport in a proton wire containing water and/or hydroxyde
groups. A succession of proton-exchange (“hop”) steps along a polarized
hydrogen-bonded chain results in the translocation of an ionic defect, H,
from end to end. This relay process leaves the chain in the opposite
orientation, so that inversion of the chain (“turn”) must take place to
complete the translocation cycle. This latter process involves the direc-
tional migration of a bonding defect triggered by the reorientation of each
hydrogen-bearing group in the chain.
FIGURE 2 -helix structure of the gA dimer. The narrow cylindrical
pore is lined with peptide bonds and accommodates a single-file water
chain of water molecules, depicted here as red spheres. This picture was
generated with the VMD software (Humphrey et al., 1996).
Proton Relay in Gramicidin 2305
Biophysical Journal 82(5) 2304–2316

tum effects arising form the light mass of H, and the
modulation of proton-transport properties by a protein ma-
trix have been explored in simple, yet biologically-relevant
systems.
The Grotthuss mechanism in protonated chains of water
molecules forming a linear hydrogen-bonded array in non-
polar channels was the object of several computational
studies (Pome`s and Roux, 1995, 1996a, 1998; Pome`s, 1999;
Drukker et al., 1998; Decornez et al., 1999; Mei et al., 1998;
Sadeghi and Cheng, 1999; Brewer et al., 2001). The hop-
ping of H was found to be dominated by structural fluc-
tuations modulating donor-acceptor separations in the hy-
drogen-bonded chain that take place spontaneously at 300
K. These fluctuations drive the exchange between OH3-
like and O2H5-like species (Pome`s and Roux, 1995). The
translocation of protons across several water molecules may
occur in subpicosecond time scales (Pome`s and Roux,
1996a; Sadeghi and Cheng, 1999). Nuclear quantum effects
(zero-point energy and quantum tunneling of hydrogen nu-
clei) on the equilibrium structure of hydrogen-bonded
chains of water molecules have been studied with dis-
cretized Feynman path integral for the treatment of ex-
changing protons. These effects, although significant, do
not govern the transfer process in equilibrium conditions
(Pome`s and Roux, 1995, 1996a; Mei et al., 1998; Brewer et
al., 2001). By contrast, under nonequilibrium initial condi-
tions mimicking the effect of an external electric field
(Drukker et al., 1998) and of hydrogen-bonding partners
restricting the displacements of water molecules (Decornez
et al., 1999), nuclear tunneling and nonadiabatic transitions
may play an important role in the translocation.
Studies of bulk water (Tuckerman et al., 1995) and gram-
icidin (Pome`s and Roux, 1996b) have revealed how struc-
tural fluctuations of the hydrogen-bonded network funda-
mentally dominate the rapid, passive relay of H. More
specifically, changes in the hydrogen-bond connectivity of
water molecules control the progress of ionic translocation
in these systems. In bulk water, such changes consist of
making and breaking hydrogen bonds in the second hydra-
tion shell of H. In gA, proton hopping appears to be
limited by the migration of defects in the polarization of the
wire. These defects result from hydrogen bond interactions
between water molecules and carbonyl oxygen atoms lining
the pore interior. Accordingly, comparisons between the
proton translocation mechanism in gA and in model hydro-
phobic pores suggest that the Grotthuss mechanism is
highly sensitive to the detail of hydrogen-bonding and elec-
trostatic interactions between the water wire and the channel
(Pome`s and Roux, 1996b; Pome`s, 1999). Finally, calcula-
tions of the free energy or potential of mean force (PMF) for
both hop and turn steps of the Grotthuss mechanism in
nonpolar channels indicated that the reorientation of the
wire, unlike hopping, is a thermally activated process
(Pome`s and Roux, 1998), suggesting that the reorganization
of the wire, not the passage of H itself, limits the rate of
proton translocation in these simplified models.
In this article, the molecular mechanism governing both
ionic and bonding translocation steps in gA is investigated.
The molecular dynamics and free energy simulation ap-
proaches used previously are applied to the chain of water
molecules embedded in gA. The structure, dynamic fluctu-
ations, and thermodynamic properties of the hydrogen-
bonded network formed by the water molecules are com-
puted successively with and without an excess proton. The
present study focuses on the effect of the gramicidin chan-
nel on proton conduction. To this end, a detailed analysis of
the hop and turn process is presented and compared with the
results obtained previously in nonpolar channels.
MATERIALS AND METHODS
Molecular model
The molecular model used in the molecular dynamics simulations de-
scribed below consists of the head-to-head pentadecapeptide dimer form-
ing the gA pore, together with water molecules. The system contains three
well-defined sections: the polypeptidic dimer, the single-file water chain,
and two cylindrical caps of water molecules lying outside the mouths of the
channel. The three-dimensional structure of gA has been determined from
solid-state nuclear magnetic resonance spectroscopy (Arseniev et al., 1985;
Ketchem et al., 1997) (Fig. 2). The starting configuration of the channel
was taken from previous molecular dynamics simulations in which the
lipid membrane was modeled explicitly (Woolf and Roux, 1994). As in
previous simulations (Pome`s and Roux, 1996b), harmonic restraining
potentials (with a force constant of 0.1 kcal/mol/Å2) were imposed on the
heavy atoms of the eight Trp indole rings to preserve the overall fold of the
pore without affecting directly the dynamics fluctuations of backbone
atoms lining the pore interior. This constraint should be viewed as an
approximation of the restriction to indole mobility due to interactions with
the lipid membrane, which act as “anchors” of the indole rings (Woolf and
Roux, 1994). In addition, harmonic constraints with the same force con-
stant were also imposed on the heavy atoms of the peptide backbone during
high-temperature preparation of the unprotonated wire (see below) to
ensure conservation of the secondary structure of the pore.
The CHARMM force field, version 22 (Brooks et al., 1983; MacKerell
et al., 1998), was used to model protein-protein interactions. The 10 water
molecules, or proton wire, contained in the gA pore region, were modeled
with the PM6 force field of Stillinger and coworkers (Stillinger and David,
1978; Stillinger, 1979; Weber and Stillinger, 1982). PM6 is a polarizable
and dissociable model of water that consists of O2 and H moieties. It has
been used in several previous studies of proton wires (Pome`s and Roux,
1995, 1996a,b, 1998; Drukker et al., 1998; Decornez et al., 1999). This
empirical force field reproduces relevant properties of protonated water
chains (Pome`s and Roux, 1996a) and has been shown, in a comparison
with results obtained from ab initio simulations, to capture the essential
features of the mechanism of H transport in these systems (Mei et al.,
1998). The parameters used to model PM6-peptide interactions are de-
scribed elsewhere (Pome`s and Roux, 1996b).
The force fields used for the caps of water molecules lying outside the
pore differed in the protonated and the unprotonated systems. In the latter
case, 14 PM6 water molecules were used in each cap, whereas in the
former case, the TIP3P force field (Jorgensen et al., 1983) was used to
model caps of 36 water molecules each. The parameters governing PM6-
TIP3P interactions are given elsewhere (Pome`s and Roux, 1998). Such a
hybrid model offers the advantage of low computational cost compared
with an all-PM6 model and allows the inclusion of water caps sufficiently
2306 Pome`s and Roux
Biophysical Journal 82(5) 2304–2316

large to avoid finite-size effects that would preclude proton diffusion from
end to end of the single file (Pome`s and Roux, 1998). Because the TIP3P
model is not dissociable, the hybrid model eliminates the possibility of an
escape of H out of the channel during the course of the simulations.
Control simulations of the nonprotonated wire indicated that the equilib-
rium conformations of the water chain and the hydrogen-bonding coordi-
nation of interfacial water molecules obtained with the hybrid model are
consistent with those obtained with all-PM6 and all-TIP3P models.
The cylindrical caps of water molecules were carved from a periodic
box of TIP3P water equilibrated in the bulk, and superimposed onto the
outer turn of the gA monomers. The water molecules overlapping with
heavy atoms of the peptide were deleted, and in all the subsequent simu-
lations, the cap region was subjected to a boundary potential. In the
unprotonated (14-PM6 caps) and protonated (36-TIP3P caps) systems,
radial (cylindrical) restraints acted on the O atoms lying respectively
further than 6.8 and 6.0 Å from the main axis of the channel, and planar
constraints were imposed outside of the ranges 11.0  z  15.0 and
11.5  z  20 Å from the channel center. These restraining potentials
were quadratic with a force constant of 20 kcal/mol/Å2. The inner value of
the planar constraint was chosen so as to avoid the artifact of water
diffusion in the nonpolar region of the membrane around the outside of gA.
Similarly, planar restraints acting outside the range 11  z  11 Å (with
a force constant of 20 kcal/mol/Å2) were imposed on the 10 water mole-
cules inside the pore to prevent their diffusion into the caps. The location
of these planar boundaries was chosen so as to minimize perturbations of
the interactions between pore and cap water molecules, as determined from
unbiased simulations. No cutoff was imposed on nonbonded interactions.
Molecular dynamics simulations
The CHARMM program (Brooks et al., 1983) was used to propagate the
Langevin equation of motion. A friction coefficient of 5 ps1 was applied
to all heavy atoms in the system. After an initial equilibration, the molec-
ular systems were subjected to successive cycles of umbrella sampling
calculations comprising preparation, equilibration, and production. The
preparation stage was necessary to overcome the relatively long relaxation
times associated with the propagation of bonding defects in the gA channel,
as observed in previous simulations (Pome`s and Roux, 1996b). Uncorre-
lated configurations of the protonated wire were obtained from simulations
in which the electrostatic interactions between channel and single-file
atoms were turned off. This artificial procedure was seen to lead to the
disappearance of bonding defects and to high mobility of the ionic defect
(Pome`s and Roux, 1996b). The nonpolar channel resulting from this
procedure is similar to the model channels constructed with radial restraints
on single-file chains of water molecules, whose properties were character-
ized in previous studies (Pome`s and Roux, 1995, 1996a, 1998). In such
systems, complete translocation events from end to end of the single-file
region occur in the order of a few picoseconds so that it is easy to
equilibrate ionic defects by imposing collective reaction coordinate con-
straints (Pome`s and Roux, 1998) before turning electrostatic interactions
back on for equilibration and data collection (production). Using this
procedure helped to reach statistical convergence in the PMF calculation.
The collective reaction coordinates used to follow the progress of ionic
and bonding defects in the pore are defined as
z  
i
qizi (1)
in which the summation runs over all water atoms O and H inside the pore,
with charges qO  2e and qH  1e and cartesian coordinates zO and zH.
In the case of the unprotonated wire, z corresponds to the z component of
the total molecular dipole moment of the 10-pore water molecules, whereas
in the presence of an excess proton z/e is the position of the center of
charge along the z axis. These reaction coordinates reflect the organization
of the wire and offer the advantage of a continuous variable for translo-
cation of both ionic and bonding defects throughout the length of the pore
(Pome`s and Roux, 1998). The detail of the methodology is described below
successively for the unprotonated and protonated systems.
Unprotonated wire
The initial conformation of the (single-file) wire region was obtained by
deleting the excess proton from a previously equilibrated protonated gA
wire (Pome`s and Roux, 1996b). The water molecules were subjected to
energy minimization and thermalization (two 10-ps runs at 150 and 300 K),
whereas the peptide was held fixed, and the procedure was repeated
without fixing protein atoms. The subsequent equilibration consisted of a
60-ps simulation at 300 K, during which all single-file water-water hydro-
gen bonds pointed toward the same mouth of the channel (polarized
configuration). The calculation of the free energy profile, or PMF for the
reorientation of the wire (i.e., propagation of the bonding defect–see
Scheme 1), was obtained from a series of six molecular dynamics simu-
lations with umbrella sampling, whereby the reaction coordinate z was
constrained by a harmonic potential Vi(z)  1⁄2k (z  z,i0 )2, with a
force constant of k

 2 kcal/mol/(eÅ)2 and successive reference values of
z,i
0
 9, 8, 6, 4, 2, and 0 eÅ.
For each window i, four cycles of production were generated using the
following procedure. 1) For preparation, a 5-ps simulation was run at 400
K (with backbone harmonic constraints on the channel backbone to prevent
unraveling of the -helix). 2) For equilibration, a 5-ps simulation followed,
with T 300 K. 3) For production, the data were collected from an ensuing
20-ps simulation at 300 K. The preparation stage was required to ensure
that several distinct local minima be sampled that correspond to a given
value of z0. Given the long characteristic times of reorganization of the
hydrogen-bonded network, which was inferred to be of the order of 100 ps
from earlier simulations of the gA channel (Pome`s and Roux, 1996b),
much longer simulations would be required to reach statistical convergence
at 300 K. The higher temperature used in the preparation stage enables a
larger number of conformations to be sampled in the relatively short
simulation times amenable to our PMF calculation. For each window in the
umbrella sampling calculation, 32,000 configurations were obtained at
2.5-fs intervals from a total of 80 ps of simulations. These configurations
were used in the computation of the PMF (see below). The calculations
were repeated with TIP3P single-file water molecules (Jorgensen et al.,
1983) to gauge the influence of the water model in the calculations, as in
similar calculations performed on nonpolar channels (Pome`s and Roux,
1998; Pome`s, 1999).
Protonated wire
The starting configuration of the channel with protonated water wire was
taken from a previous simulation (Pome`s and Roux, 1996b). Two cylin-
drical caps of 36 water molecules lying outside the mouths of the channel
were added to the system. The energy was minimized (300 steps of steepest
descent) and the caps were briefly thermalized (5 ps, 300 K) with the rest
of the system held fixed. Eight series of biased simulations or windows
were then created, each with a quadratic restraint Vi(z) differing by their
reference z,i0 (see previous subsection). The reference values in the suc-
cessive windows were z,i0 0.5, 0, 0.5, 1, . . . , 2.5, and 3 eÅ. The latter
value corresponds to wires in which the excess proton is located at the
interface between bulk and single-file regions (Pome`s and Roux, 1998).
For each of these eight windows, the simulation protocol consisted of
four iterations of the following cycle. 1) For preparation of the protonated
wire, a 300-step SD energy minimization was followed by 5-ps MD
equilibration of the proton in the single-file water wire at 350 K with
wire-channel electrostatic interactions turned off. 2) For equilibration of
the channel, a 300-step SD minimization was followed by 25-ps MD at 300
K with full electrostatic interactions. 3) Production consisted of 50 ps of
data collection at 300 K. This procedure ensured the generation of four
Proton Relay in Gramicidin 2307
Biophysical Journal 82(5) 2304–2316

uncorrelated production cycles to improve configurational averaging. Pro-
duction times of 200 ps were generated for each window, for a total of 2560
ps of simulations including preparation, equilibration, and production.
PMF calculations
The wire configurations produced in the simulations of protonated and
unprotonated systems were used to calculate the PMF for ionic and
bonding defect translocations, respectively. To this end, the weighted
histogram analysis method (Kumar et al., 1992, Roux, 1995) was used as
in similar calculations reported previously (Pome`s and Roux, 1998). Sta-
tistical convergence of the PMF was verified by comparing the profiles
obtained from various subsets of configurations.
RESULTS AND DISCUSSION
Unprotonated wire
Equilibrium structure
The gA channel is a narrow cylindrical pore filled with an
array of water molecules arranged as a single file (Arseniev
et al., 1985; Ketchem et al., 1997; Fig. 2). The pore water
molecules form a hydrogen-bonded network comprising
both water-water and water-channel hydrogen bonds (Roux
and Karplus, 1994, references therein; Pome`s and Roux,
1996b). A water molecule is considered to belong to the
single file if 1) it makes at most two hydrogen bonds with
other water molecules and 2) at least one of these two
coordinating water molecules is itself a single-file water.
Thus defined, the single-file region contains eight water
molecules. In addition, two outlying water molecules form
an interface with water molecules outside the pore. Each of
the eight single-file water molecules is surrounded by two
adjacent water molecules with which it can either donate or
accept a hydrogen and carbonyl O atoms of the channel to
which it can donate a hydrogen atom. The single-file water
molecules do not in general accept hydrogen atoms from the
channel, because backbone NOH bonds, which are in-
volved in intramolecular hydrogen bonds defining the -he-
lix fold, do not deviate sufficiently from colinearity with the
channel axis to point toward water O atoms in the lumen.
A representative conformation of the water chain is de-
picted in Fig. 3. Most water molecules in the wire adopt the
following organization: donation of one H to a channel
backbone carbonyl, donation of the other H to a neighboring
water, and acceptance of one H from the other neighboring
water, for a total of three hydrogen bonds. Such a coordi-
nation results in a continuous chain of water-water hydro-
gen bonds. In the single-file region, water OOH bonds
involved in hydrogen bonding (donation) to the channel are
approximately perpendicular to the channel axis, and most
of the covalent OOH bonds making up the water-water
hydrogen-bonded chain point toward the same entrance of
the gA dimer: the hydrogen-bonded chain is preferentially
polarized. The projection of the dipole moment of each of
these water molecules along the z axis is approximately zi
 1 eÅ. In addition to the polarized chain, in general the
wire also contains exactly one water molecule that donates
both of its H atoms to the channel. Such coordination forces
the plane of that water molecule to be approximately per-
pendicular to the channel axis, with zi  0. In general, the
two water molecules adjacent to this perpendicular water
point one of their OOH bonds toward it. The resulting
inversion in the topology of the hydrogen-bonded network
defines what we shall refer to as a bonding defect.
The average hydrogen-bonding coordination of each of
the 10 water molecules in the wire was computed from a
20-ns simulation of the unprotonated wire with the TIP3P
water model (Table 1). Results obtained from shorter sim-
ulations with the PM6 model (not shown) are highly similar.
The asymmetry in the coordination of pairs of water mole-
cules with indices i and (11  i) indicates that statistical
convergence is not achieved within 20 ns. Nevertheless, this
analysis underlines important features in the organization of
the hydrogen-bonded network. The preferred arrangement
of the wire is reflected in the dominance of a hydrogen-
bonding coordination of three for the single-file water mol-
ecules (indices 2–9). The seemingly high occurrence of null
H bond coordination between water and channel (12%)
reflects the fact that a multiple choice of carbonyl O atoms
often results in bifurcated H bonds, which are not counted.
The bonding defect, characterized by donation of 2 H to the
channel, is preferentially at water 2, 3, 8, or 9, near the end
of the single file. This bonding defect is not necessarily a
FIGURE 3 Representative structure of the hydrogen-bond network of
water molecules inside gA without an excess proton. Seven of the 10 water
molecules located inside the pore are shown. All of them are tricoordinated
and donate one H atom to carbonyl O atoms of the channel, except for
water 2, which is oriented perpendicularly to the channel axis and forms
four hydrogen bonds, donating both of its H atoms to the channel. The
polarization of water molecules inverts around water 2, which defines a
bonding defect.
2308 Pome`s and Roux
Biophysical Journal 82(5) 2304–2316

hydrogen-bonding defect: although water molecules 2, 3, 8,
and 9 are more likely than the other single-file water mol-
ecules to be engaged in only one water-water bond, they are
also more likely to form a total of four H bonds. However,
even if the hydrogen-bonded chain is continuous, the local
inversion of its polarity precludes passage of protons in the
sense of a Grotthuss mechanism (see Fig. 1). Finally, inter-
facial water molecules (indices 1 and 10) often make two
hydrogen bonds with outlying water molecules, preferably
as acceptors. This preferred orientation reflects the presence
of an intervening bonding defect.
The preferred organization of the unprotonated wire as a
polarized water chain is not due to the gA channel, which,
by virtue of its -helical secondary structure and its assem-
bly as a head-to-head homodimer, possesses no net dipole
moment. Rather, the polarization is an intrinsic property of
the single-file water chain, which was attributed to the
maximization of favorable water dipole-dipole interactions
in nonpolar channels (Pome`s and Roux, 1998; Pome`s,
1999). The presence of a bonding defect in the preferred
conformations of the unprotonated water wire in gA con-
trasts with the fully polarized chains obtained in hydropho-
bic channel models. Full polarization is driven by the opti-
mization of dipole-dipole interactions between the single-
file water molecules. The preferred location of the bonding
defect, near the end of the single file, maximizes the length
of the polarized segment in gA.
Dynamic fluctuations
Pore water molecules can adopt three distinct polarization
states corresponding respectively to zi  1, 1 (polar-
ized), and 0 eÅ (bonding defect). The reorientation of water
molecules gives rise to unitary transitions between these
three states, which results in the displacement of the bond-
ing defect. This process, which occurs spontaneously in the
simulations, is depicted schematically in Fig. 4. The con-
formations of the wire include metastable states in which
one water molecule is a bonding defect, as well as transient
conformations in which two adjacent water molecules are
perpendicular to the channel axis. The unitary translocation
of a bonding defect arises from the succession of two
elementary reorientation steps. First, a water molecule ad-
jacent to the bonding defect (zi1  0) reorients from a
polarized state (zi  1) to perpendicular (zi  0). This
process replaces one water-to-water H bond by a water-to-
channel H bond, which creates a hydrogen-bonding defect,
in a configuration which is sometimes referred to as nega-
tive Bjerrum defect, whereby two adjacent water O atoms
are without an intervening H atom. Inversely, in the second
step water i  1 reorients from zi1  0 to 1, which
eliminates the hydrogen-bonding defect in the wire and
completes the unitary progression of the bonding defect.
Thus, the hydrogen-bonding coordination provided by gA
determines at once the nature of the bonding defect and the
detailed mechanism for its migration (Pome`s, 1999). It is
the alternation of the two types of bonding defects that
mediates the unitary migration of a bonding defect in the gA
channel.
The dynamics of dipole reorientation of the pore water
molecules is illustrated in Fig. 5. Initially, the 2nd water
molecule from the bottom is a bonding defect (as depicted
in Fig. 3). Sequences of reorientations trigger the migration
TABLE 1 Hydrogen-bonding coordination* of pore water
molecules averaged over a 20-ns MD simulation (in
percentage points)
Water
(#)
Water to
channel Water-water Number of H bonds
0 1 2 1 2 3 1 2 3 4
1 9 85 5 11 47 42 1 10 45 44
2 8 58 34 25 75 0 2 19 58 21
3 11 71 19 19 81 0 2 19 67 12
4 11 86 3 11 89 0 2 17 79 2
5 16 82 2 9 91 0 2 20 76 2
6 19 79 1 12 88 0 2 24 71 2
7 12 86 2 14 86 0 2 20 75 2
8 14 77 9 16 84 0 2 21 70 7
9 13 71 16 17 83 0 2 21 67 10
10 9 88 3 8 48 44 1 9 45 45
Average 12 78 10 15† 85† 0† 2† 20† 70† 6†
*Assuming that a hydrogen bond is formed whenever the donor-acceptor
separation is less than 3.2 Å and the donor-H-acceptor angle is greater than
120 degrees.
†Averages performed on the single-file water molecules (indices 2–9) only.
FIGURE 4 Schematic representation of the unitary translocation of a
bonding defect. Thick lines represent water molecules forming hydrogen
bonds (dashed lines) with each other and with the channel. (Top) A
bonding defect is on the 2nd water molecule from the right; (Middle)
reorientation of water 3 creates a hydrogen-bonding defect between water
molecules 2 and 3; (Bottom) subsequent reorientation of water 2 moves the
bonding defect to water 3.
Proton Relay in Gramicidin 2309
Biophysical Journal 82(5) 2304–2316

of the bonding defect between its two preferred locations,
i.e., near the mouths of the pore. These transitions are
infrequent and involve metastable states corresponding to
bonding defects (zi  0) located at water molecules 3
through 8. Thus, the bonding defect migrates from water 2
to water 5 at t  503 ps, remains there for 5 ps, and
proceeds to water 8, which inverts the polarization of the
chain. The reverse process takes place at 522  t  530 ps,
with intermediate stations at water molecules 7 and 6. The
time evolution of the collective reaction coordinate for the
translocation of the bonding defect, z 
i zi, is depicted
at the bottom of Fig. 5. The total dipole moment of the chain
is highly sensitive to transitions in the orientation of indi-
vidual water molecules and reflects the position of the
bonding defect. The cooperativity of water reorientation
was noted in earlier simulations (Chiu et al., 1989), and it
was shown that the process is influenced by fluctuations in
the conformation of the channel (Chiu et al., 1991). Al-
though thermal fluctuations of the chain take the bonding
defect (and the polarization of the chain) back and forth
from end to end of the single file, this process is infrequent
(it occurs only 12 times in the 20-ns simulation). This,
together with cooperativity, suggests an energy-activated
process.
PMF for the propagation of a bonding defect
The free energy profile for the translocation of a bonding
defect, as calculated with umbrella sampling, is shown in
Fig. 6. The results obtained for the two water models used,
TIP3P and PM6, are qualitatively similar to each other. The
preferred location of the bonding defect at water molecules
2 and 9 is reflected in the presence of a free energy well at
z  6.5 eÅ, which corresponds to mostly-polarized
conformations in which eight of the ten pore water mole-
cules form an oriented hydrogen-bonded chain. The inter-
conversion between these two polarized conformations in-
volves an activation energy barrier centered at z  0,
which corresponds to a hydrogen-bonding defect located
between water molecules 5 and 6 (z5  z6  0). The
barrier height is 3.8 kcal/mol with the polarizable water
model PM6 and 2.2 kcal/mol with the TIP3P model. The
PMF profile obtained from earlier simulations of a nonpolar
channel of nine PM6 water molecules (Pome`s and Roux,
1998), also shown in Fig. 6, is qualitatively similar but
involves a much larger free energy barrier. At 7.6 kcal/mol,
this barrier is about twice as large as that obtained for PM6
in the gA channel. Thus, the primary effect of the channel
appears to be to catalyze the propagation of a bonding
defect throughout the single file.
The PMF obtained with the PM6 model in gA includes
six secondary minima located at z  0.5, 2.3, and
4.6 eÅ. Each of these secondary minima corresponds to
metastable conformations in which one of the water mole-
cules between water 2 and water 9 is in the process of
flipping (see above). Similarly, the inflection points ob-
served at z   8.1 eÅ correspond to conformations in
which water molecules 1 and 10 are perpendicular to the
channel axis. The PMF profile obtained with the TIP3P
model also shows inflection points and secondary minima,
albeit not as pronounced. By contrast, the PMF for reorien-
tation of a wire in a nonpolar channel is much smoother,
which is due to the absence of stabilizing interactions on the
FIGURE 5 (Top) Time evolution of individual water dipole moment
projections on the channel (z) axis obtained from molecular dynamics
simulation of the unprotonated water wire. Each of the dipole components
fluctuates in the range 1  zi  1 eÅ; each is shifted by (i  1) eÅ
for clarity. (Bottom) z-component of the total dipole moment of the wire for
the same time window.
FIGURE 6 Potential of mean force for the reorientation of the unproto-
nated water chain in the single-file region of the GA channel: (solid)
polarizable (PM6) water wire; (dashed) TIP3P water wire. Note that in the
latter case, z was scaled by 0.417, the partial charge of H atoms in the
TIP3P model (Jorgensen et al., 1983), for direct comparison with the PM6
model, in which the formal charge of H atoms is 1. Results of the same
calculation performed with nine PM6 water molecules in an analogous
nonpolar channel (Pome`s and Roux, 1998) are also shown (dotted) for
comparison.
2310 Pome`s and Roux
Biophysical Journal 82(5) 2304–2316

channel wall. These results suggest that the physical origin
for catalysis of the turn step in the gA channel is identical to
that which gives rise to a bonding defect: by accepting two
hydrogen bonds from a water molecule, the carbonyl O
atoms lining the pore interior stabilize intermediate arrange-
ments of the water wire in which a water molecule is in the
process of reorienting or flipping between its two preferred
(polarized) states.
Protonated wire
Structure
A representative structure of the protonated water chain is
shown in Fig. 7. The presence of the excess proton modifies
the strength and the topology of water-water hydrogen
bonds. As discussed in detail elsewhere (Pome`s and Roux,
1995, 1996a,b), in a single-file environment the hydrated
proton can generally be described either as a hydronium,
OH3
, or as a “Zundel cation,” O2H5
, in which the proton is
shared by two water molecules brought together by the
excess charge in a strong hydrogen bond. At 2.4 Å, this
strong hydrogen bond is much shorter than a typical water-
water hydrogen bond (in the absence of an excess proton).
In gA, as in nonpolar water-filled pores, the polarization of
the water chain induced by the excess charge extends be-
yond the immediate vicinity of the ion. The excess proton is
surrounded on either side by two oriented chains of water
molecules pointing their O atom toward the positive charge.
By contrast to water-water interactions, the topology of
the hydrogen-bonded network formed between water mol-
ecules and the channel is not affected significantly by the
excess proton. In particular, protonated water molecules
retain a hydrogen-bonding coordination of exactly three, as
depicted in Fig. 7. In general, the hydrogen-bonded chain
does not extend through the entire length of the single-file
region of the channel. Instead, interruptions of the chain, or
bonding defects, are observed sufficiently far from the ex-
cess proton (Pome`s and Roux, 1996b). When the proton is
near the middle of the wire (as in Fig. 7), these defects are
located preferentially at either end of the single-file region,
whereas when the excess proton is closer to one of the
mouths, the defect is located near the other mouth.
Dynamic fluctuations
The rapid exchange of H in the chain of water molecules
in gA occurs spontaneously with thermal fluctuations at 300
K (Pome`s and Roux, 1996b). As depicted schematically in
Fig. 8, this exchange involves transitions between succes-
sive hydronium-like and Zundel-cation-like arrangements
of the water molecules. Accordingly, the dynamic relay of
H in the proton wire may be easily followed either by
monitoring at every snapshot of the simulation, which O
atom is closest to three H nuclei (hydronium-like water
molecule), or alternatively, by monitoring the position of
the shortest OO separation (proton-sharing water dimer).
These two reaction coordinates are complementary ways to
view the translocation process (Pome`s and Roux, 1996b).
However, neither of these representations describes the
translocation process fully, because they implicitly assume
that hopping of the excess proton is entirely determined by
the local arrangement of the one or two water molecule(s)
closest to it. Rapid cascades in the movement of the hy-
drated proton over as many as six or seven water molecules
were observed occasionally, suggesting that collective fluc-
tuations are important to long-range transport (Pome`s and
Roux, 1996b).
Although the two steps of the Grotthuss mechanism are
studied separately in the present work, it should be stressed
that the propagation of H is not independent of reorienta-
tion steps. The migration of bonding defects standing in the
FIGURE 7 Representative structure of the hydrogen-bond network of a
protonated water wire in gramicidin. In this snapshot, the excess proton
(ionic defect) is in the form of a hydronium H3O
 and the bonding defect,
highlighted by a green dot, is a hydrogen-bonding defect.
FIGURE 8 Schematic representation of the hydrogen-bonded network
coordination of both (top) hydronium and (bottom) Zundel forms of pro-
tonated water molecules in the gA channel.
Proton Relay in Gramicidin 2311
Biophysical Journal 82(5) 2304–2316

path of H is a prerequisite to complete proton transloca-
tion, because H transfer necessitates preexisting and suit-
ably-oriented hydrogen bonds. Thus, the slower rate of
transport of bonding defects apparently limits the rate of
translocation of the ionic defect in gA (Pome`s and Roux,
1996b). Conversely, it might be expected that the entry of
H at one end of an unprotonated single file would hasten
the translocation of the bonding defect via stabilization of
one polarized conformation of the wire relative to the other.
The polarization of water molecules induced by the pres-
ence of Na at the entrance of the channel was noted in an
earlier computational study of gA (Jordan, 1990).
The displacement of ionic and bonding defects occurs on
different time scales. The migration of bonding defects is
governed by water reorientation events that were observed
to occur infrequently in the range of 100 ps or longer. This
is significantly slower than the movement of H, which
takes place in the picosecond or subpicosecond time range
(Pome`s and Roux, 1996b). For this reason, the average
position of H may remain in the vicinity of its starting
point over simulations of a few hundred picoseconds. In the
present work, successive simulations for which the excess
proton was equilibrated at various positions were used to
overcome the separation of time scales between ionic and
bonding defect translocation. The collective reaction coor-
dinate, z, used in the present study reflects the polarization
of the wire by taking into account not only the position of
the excess proton but also the relative distribution of O and
H nuclei along the channel axis z (Pome`s and Roux, 1998).
Thus, this reaction coordinate incorporates the presence of
bonding defects in gA.
PMF for the propagation of an ionic defect
Because bonding defects move on a time scale comparable
with that of the simulations, it is important that the ensemble
of conformations used for the calculation of the PMF cover
not only the displacement of H but also of the bonding
defect. This proved essential to the proper convergence of
the PMF for the migration of an excess proton from end to
end of the single-file region, which is shown in Fig. 9. This
profile represents the reversible thermodynamic work for
the complete translocation of the protonic defect. At z 0,
the center of charge is located on average near the center of
the channel, whereas numerical values of z  3 eÅ
correspond to configurations in which H is located at the
extremity of the polarized water chain with the ionic defect
near z  10 Å and the total dipole moment of the wire
pointing away at 7 Å. The largely activation-less PMF
profile suggests a diffusive mechanism for proton hopping.
The shallow well centered at the origin indicates that the
excess proton is somewhat better “solvated” near the dimer
junction with a moderate preference over outlying single-
file locations reflected by a 1-kcal/mol drop in free energy.
This could be due in part to the propensity of extremal water
molecules to form bonding defects. A corollary to the po-
larization of water around the excess proton is that H is
always hosted by polarized water molecules, never by a
bonding defect. In addition, finite-size effects could lead to
an artificial bias in proton location. In studies of nonpolar
channels without water caps, the excess proton remained
localized near the center of the wire, whereas in the pres-
ence of droplets of 25 water molecules, H was evenly
distributed (Pome`s and Roux, 1998). A similar effect was
also observed in another study (Mei et al., 1998). In the
present study, increasing the size of the caps from 36 to 64
water molecules each (results not shown) did not affect
significantly the equilibrium distribution of H in the single
file. This indicates that the water caps used in the present
study are adequately large and do not underlie the bias in
proton location.
The PMF obtained for a nonpolar homolog of the channel
(Pome`s and Roux, 1998) is also depicted in Fig. 9. This
profile, which was obtained in a very similar water system
consisting of a single file of nine water molecule between
spherical droplets of 25 water molecules each (“dumbbell”
model), is a broad single well in the region 2.0  z 
2.0 eÅ. This result indicates that in the absence of hydro-
gen-bonding partners to water molecules in the pore, no
work is required to relay an excess proton from one end of
the single file water chain to the other. This result is largely
reproduced in gA, which suggests that the local environ-
ment provided by the channel is well suited for proton
mobility. Thus, based on the present results it appears that
the periodicity of groups presented by gramicidin provides
no binding site or “trap” of protons in the single-file region.
Gramicidin as a proton duct
Our choice of molecular system reflects the focus of this
work on the specific contributions made by the peptide on
FIGURE 9 Potential of mean force for the translocation of an excess
protonic charge (solid line) in the lumen of gramicidin and (dashed line) in
an analogous, nonpolar channel (Pome`s and Roux, 1998).
2312 Pome`s and Roux
Biophysical Journal 82(5) 2304–2316

proton conduction through a comparison of the molecular
mechanisms obtained in nonpolar channels and in gA. This
incremental, comparative approach presents the advantage
of mitigating the effect of systematic errors inherent to
empirical energy functions. Clearly, this approach also lim-
its the scope of our conclusions. The absence of the lipid
membrane precludes a realistic treatment of long-range
electrostatic interactions, which play a role in the translo-
cation of ions. Based on studies of proton permeation in
which Trp residues were fluorinated and replaced by Phe
residues, it has been proposed that dipole-dipole interac-
tions between water molecules and indole rings could play
a role in the preferred arrangement of water molecules in the
channel, particularly near the mouths of the pore (Phillips et
al., 1999). Our simplified representation of the bulk-water/
channel interface makes the model inadequate to the treat-
ment of such effects, as well as of the entrance and exit of
protons and of water molecules.
The present results highlight the effect of water-water and
water-channel hydrogen-bonding interactions at play in gA
(Pome`s and Roux, 1996b). The fine balance between these
forces modulates the structure and dynamic fluctuations of
the hydrogen-bonded network, which in turn govern the
mechanisms of translocation of bonding and ionic defects in
the single-file region of the channel. For a valid description
of the mechanism it is therefore essential that the relative
strengths of water-water and water-channel interactions be
adequately modeled. Quantitative discrepancies in the acti-
vation energy for the reorientation of the water wire in gA
(Fig. 6) obtained with PM6 and TIP3P reflect the different
nature of these two water models. Despite this difference,
the qualitative agreement obtained for the mechanism of
translocation of a bonding defect indicates that both models,
which exhibit identical mechanisms of reorientation in non-
polar channels (Pome`s, 1999), also respond consistently to
gA. In particular, the two models are in agreement regarding
the coordination of water molecules, the nature of the bond-
ing defect, and its preferred distribution in the unprotonated
wire.
The relationship between water coordination and hop
versus turn mobility in proton wires has been discussed
elsewhere in a comparison of the Grotthuss mechanism in
liquid water, in nonpolar channels, and in gA (Pome`s,
1999). Implications of the molecular mechanism for the
permeation of H in nonpolar channels were proposed in
terms of proton leakage via transient hydrogen-bonded
chains (Pome`s and Roux, 1998). In the absence of water-
channel hydrogen-bonding interactions, each water mole-
cule in the wire is tightly coordinated to both of its neigh-
bors. As a result, there are no bonding defects in the chain
and water molecules move in concert (i.e., cooperative
motions of the O atoms are enhanced). Both of these effects
conspire to the high mobility of H in single-file arrays of
water molecules embedded in nonpolar channels (Pome`s
and Roux, 1996b). However, the large activation energy
calculated for the reorientation of water molecules
(Pome`s and Roux, 1998; Pome`s, 1999) suggests that the
turn step of the Grotthuss mechanism would be prohibi-
tively long compared with the lifetime of single-file
hydrogen-bonded arrays of water in a nonpolar cavity.
Thus, transient water chains, which have been proposed
to mediate the leakage of protons through pure lipid
membranes (Nagle, 1987; Marrink et al., 1996; Paula et
al., 1997) and to play a role in the uptake of H in proton
pumps (Wikstro¨m, 1998), would be suited at best to the
passage of only one proton before breaking apart. In such
a mechanism, the rate-limiting step for proton transport
would be the nucleation of the wire.
By contrast to leakage, efficient proton-relaying channels
such as gA provide a different mechanism for the rapid
succession of proton conduction cycles. The analysis of this
and previous simulation studies (Pome`s and Roux, 1996b;
Pome`s, 1999) offers clues as to why gA constitutes a more
effective proton duct than hypothetical hydrophobic coun-
terparts. This is achieved by assisting the proton-relay chain
in tackling the dual and seemingly contradictory require-
ment of a proton wire: to enable strong (water-water) hy-
drogen bonds between relaying groups for the rapid transfer
and relay of H and to help weaken (and break) these
hydrogen bonds so as to facilitate the reorientation of pro-
ton-relaying groups.
Locally, ideal solvation of protons in aqueous systems is
achieved by a coordination of three hydrogen-bond accep-
tors, because that is the coordination of protonated water
molecules in OH3
 and O2H5
 ions (Agmon, 1995). In the
single-file region of gA, carbonyl O atoms of the peptide
provide the extra hydrogen-bond acceptor that is missing in
nonpolar channels (Fig. 8). Thus, the effect of gA on proton
solvation (i.e., on the stabilization of an ionic defect) is
similar to that observed above regarding the two forms of
bonding defects (Fig. 4): the hydrogen-bond structure sta-
bilizes both forms of the ionic defect (and by the same
token, of larger protonated clusters of the form OnH2n1

also found in gA). Furthermore, it is because both primary
forms of the two types of defects are locally stabilized that
the migration of ionic and bonding defects, whose elemen-
tary event consists of interconversions between these two
forms (Figs. 4 and 8), is facilitated. In addition, mobility of
protons benefits from the equivalence of proton-relay
groups (i.e., the absence of proton-binding sites) in the pore.
Thanks to the periodicity of the channel backbone and to the
polarization and equal propensity of water molecules point-
ing left and right, protonic solvation is approximately as
good anywhere along the single file. Likewise, the bonding
defect is locally stabilized throughout the length of the pore.
Thus, the hydrogen bond coordination offered by the gA
channel achieves an effective compromise for the stabiliza-
tion and the mobility of both ionic and bonding defects.
Proton Relay in Gramicidin 2313
Biophysical Journal 82(5) 2304–2316

CONCLUSIONS AND PERSPECTIVES
We studied the complete translocations of a protonic defect
and of a bonding defect by water molecules in the single-file
region of the gA channel. We found that the mobility of an
excess proton in gA is essentially determined by the fine
structure and the dynamic fluctuations of the hydrogen-
bonded network. The translocation of H is mediated by
spontaneous (thermal) fluctuations in the relative positions
of oxygen atoms in the wire. In this diffusive mechanism, a
shallow free-energy well slightly favors the presence of the
excess proton near the middle of the channel. The unproto-
nated water chain adopts either one of two polarized con-
formations corresponding to oriented donor-acceptor hydro-
gen-bond pattern along the channel axis. Interconversion
between these two conformations is an activated process
that occurs via the sequential reorientation of water mole-
cules of the wire. The effect of hydrogen bonding between
channel and water on proton conduction was analyzed by
comparison to the results obtained previously in a study of
model nonpolar channels, in which such interactions were
missing (Pome`s and Roux, 1998). This analysis revealed
that hydrogen-bond donation from water to the backbone
carbonyl oxygen atoms lining the pore interior not only
offers a compromise for the solvation and the mobility of an
excess proton, it also enables the emergence of a bonding
defect and facilitates its transport.
Among the challenges facing the study of biological ion
transport mechanisms, particularly in the case of proton
movement, are the difficulty to detect ion movement exper-
imentally and to relate structure and function. The present
study has opened the way to a better understanding of ion
channel structure-function relationships by allowing the
first confrontation of a molecular-level proton translocation
mechanism with experimental measurements. Based on the
results presented here, a framework model was developed
for single-proton transport flux through gramicidin (Schu-
maker et al., 2000, 2001). That model describes the trans-
port of H (ionic defect) and of a bonding defect by poten-
tial of mean-force profiles (Figs. 6 and 9) and diffusion
coefficients obtained from molecular dynamics simulations.
Proton entrance and exit in and out of the single-file region,
which were not studied by molecular dynamics in the
present model, are represented parametrically. A reasonable
choice of these parameters yields a good fit to experimental
proton conductance data (Eisenman et al., 1980).
Theoretical approaches have begun to provide unprece-
dented and meaningful insight on the transient events gov-
erning proton relay mechanisms. Whereas the detailed stud-
ies performed to date on biological systems have been
limited to the water chain embedded in the gramicidin
channel, the results of these studies provide both a frame-
work for understanding the basic physical principles at play
in water wires and an impetus for the study of relay mech-
anisms in more complex proton wires. What is emerging
from this and preceding studies of water wires is that the
control of proton movement in biomolecular systems may
be achieved with subtle local structural fluctuations of hy-
drogen bonded networks of low dimensionality, which are
themselves determined by the fine hydrogen-bonding coor-
dination, arrangement, and topology of proton-relay groups.
The low dimensionality of a biological proton wire offers
the possibility of a tight control of the geometry and the
topology of hydrogen-bond interactions. The low-amplitude
fluctuations of a flexible array (the proton wire) are thus
mastered by a comparatively rigid environment presenting a
few “judiciously disposed” hydrogen-bond partners.
Efficient proton conduction requires a wire that can al-
ternatively form strong hydrogen bonds and break them
relatively easily with thermal fluctuations. Water molecules
are well suited for that, because water is ubiquitous in
biological systems and forms highly modulable hydrogen-
bonded networks. gA has harnessed those properties to
assist both hop and turn steps of the Grotthuss mechanism.
By the same token, it may be expected that somewhat
similar features in the coordination of hydrogen-bond net-
works will be observed in other biomolecules that provide
efficient ducts for the passive long-range relay of H. In
particular, based on the present study the detailed charac-
terization of the structure of hydrogen-bonded arrays of
water molecules embedded in channels and enzyme cavities
might offer a clue as to their ability to mediate rapid proton
transport. Inversely, it is likely that other proteic matrices
make use of the coordination properties (both static and
dynamic) of water molecules and of other titratable groups
to control the long-range transport of protons in other ways:
for the selectivity, delay, gating, pumping, and blockage of
protons.
Stimulating discussions with Samuel Cukierman, Mark Schumaker, and
Ching-Hsing Yu are gratefully acknowledged. We also thank C.-H. Yu for
his help in the preparation of Table 1. This work was supported by a grant
from the Canadian Institutes of Health Research. R.P. is a CRCP
Chairholder.
REFERENCES
Agmon, N. 1995. The Grotthuss mechanism. Chem. Phys. Lett. 244:
456–462.
Akeson, M., and D. W. Deamer. 1991. Proton conductance by the gram-
icidin water wire. Biophys. J. 60:101–109.
Arseniev, A. S., I. L. Barsukov, V. F. Bystrov, A. L. Lomize, and Y. A.
Ovchinnikov. 1985. 1H-NMR study of gramicidin-A transmembrane ion
channel: head-to-head right handed single stranded helices. FEBS Lett.
186:168–174.
Brewer, M. L., U. W. Schmitt, and G. A. Voth. 2001. The formation and
dynamics of proton wires in channel environments. Biophys. J. 80:
1691–1702.
Brooks, B. R., R. E. Bruccoleri, B. D. Olafson, D. J. States, S. Swami-
nathan, and M. Karplus. 1983. CHARMM: a program for macromolec-
ular energy minimization and dynamics calculations. J. Comput. Chem.
4:187–217.
2314 Pome`s and Roux
Biophysical Journal 82(5) 2304–2316

Chiu, S.-W., E. Jacobsson, S. Subramaniam, and J. A. McCammon. 1991.
Time-correlation analysis of simulated water motion in flexible and rigid
gramicidin channels. Biophys. J. 60:273–285.
Chiu, S.-W., S. Subramaniam, E. Jacobsson, and J. A. McCammon. 1989.
Water and polypeptide conformations in the gramicidin channel. Bio-
phys. J. 56:253–261.
Cukierman, S. 2000. Proton mobilities in water and in different stereoiso-
mers of covalently linked gramicidin A channels. Biophys. J. 78:
1825–1834.
Decornez, H., K. Drukker, and S. Hammes-Schiffer. 1999. Solvation and
hydrogen-bonding effects on proton wires. J. Phys. Chem. A 103:
2891–2898.
DeCoursey, T. E., and V. V. Cherny. 1999. An electrophysiological com-
parison of voltage-gated proton channels, other ion channels, and other
proton channels. Isr. J. Chem. 39:409–418.
Dencher, N. A., H. J. Sass, and G. Bu¨ldt. 2000. Water and
bacteriorhodopsin: structure, dynamics, and function. Biochim. Biophys.
Acta. 1460:192–203.
Doyle, D. A., J. M. Cabral, R. A. Pfuetzner, A. Kuo, J. M. Gulbis, S. L.
Cohen, B. T. Chait, and R. MacKinnon. 1998. The structure of the
potassium channel: molecular basis of K conduction and selectivity.
Science. 280:69–77.
Drukker, K., S. W. de Leeuw, and S. Hammes-Schiffer. 1998. Proton
transport along water chains in an electric field. J. Chem. Phys. 108:
6799–6808.
Eisenman, G., B. Enos, J. Ha¨gglund, and D. Sandblom. 1980. Gramicidin
A as an example of a single filing ionic channel. Ann. N. Y. Acad. Sci.
329:8–20.
Finkelstein, A. 1987. Water movement through lipid bilayers, pores, and
plasma membranes. In Theory and Reality. John Wiley & Sons, New
York.
Grotthuss, C. J. T. de. 1806. Me´moire sur la de´compostion de l’eau et des
corps qu’elle tient en dissolution a` l’aide de l’e´lectricite´ galvanique. Ann.
Chim.58:54–74.
Humphrey, W., A. Dalke, and K. Schulten. 1996. VMD: visual molecular
dynamics. J. Mol. Graphics. 14:33–38.
Jordan, P. C. 1990. Ion-water and ion-peptide correlations in a gramicidin-
like channel. Biophys. J. 58:1133–1156.
Jorgensen, W. L., J. Chandrasekhar, J. D. Madura, R. W. Impey, and M. L.
Klein. 1983. Comparison of simple potential functions for simulating
liquid water. J. Chem. Phys. 79:926–935.
Kandori, H. 2000. Role of internal water molecules in bacteriorhodopsin.
Biochim. Biophys. Acta. 1460:177–191.
Ketchem, R. R., B. Roux, and T. A. Cross. 1997. High resolution refine-
ment of a solid-state NMR-derived structure of gramicidin A in a lipid
bilayer environment. Structure. 5:11655–11669.
Knapp, E.-W., K. Schulten, and Z. Schulten. 1980. Proton conduction in
linear hydrogen-bonded systems. Chem. Phys. 46:215–229.
Kumar, S., D. Bouzida, R. H. Swendsen, P. A. Kollman, and J. M.
Rosenberg. 1992. The weighted histogram analysis method for free-
energy calculations in biomolecules: I. The method. J. Comp. Chem.
13:1011–1021.
Lanyi, J. K. 1999. Bacteriorhodopsin. Int. Rev. Cytol. 187:161–202.
Lanyi, J. K. 2000. Molecular mechanism of ion transport in
bacteriorhodopsin: insights from crystallographic, spectroscopic, ki-
netic, and mutational studies. J. Phys. Chem. B. 104:11441–11448.
Levitt, D. G., S. R. Elias, and J. M. Hautman. 1978. Number of water
molecules coupled to the transport of sodium, potassium, and hydrogen
ions via gramicidin, nonactin, or valinomycin. Biochim. Biophys. Acta.
512:436–451.
Luecke, H. 2000. Atomic resolution structures of bacteriorhodopsin pho-
tocycle intermediates: the role of discrete water molecules in the func-
tion of this light-driven ion pump. Biochim. Biophys. Acta. 1460:
133–156.
MacKerell, A. D., Jr., D. Bashford, M. Bellot, R. L. Dunbrack, J. D.
Evanseck, M. J. Field, S. Fischer, J. Gao, H. Guo, S. Ha, D. Joseph-
McCarthy, L. Kuchnir, K. Kuczera, F. T. K. Lau, C. Mattos, S. Mich-
nick, T. Ngo, D. T. Nguyen, B. Prodhom, W. E. Reither, I. I. I., B. Roux,
M. Schlenkrich, J. C. Smith, R. Stote, J. Straub, M. Watanabe, J.
Wio´rkiewicz-Kuczera, and M. Karplus. 1998. All-atom empirical po-
tential for molecular modeling and dynamics studies of proteins. J. Phys.
Chem. B. 102:3586–3616.
Marrink, S. J., F. Ja¨hnig, and H. J. C. Berendsen. 1996. Proton transport
across transient single-file water pores in a lipid membrane studied by
molecular dynamics simulations. Biophys. J. 71:632–647.
Mei, H. S., M. E. Tuckerman, D. E. Sagnella, and M. E. Klein. 1998.
Quantum nuclear ab initio molecular dynamics study of water wires.
J. Phys. Chem. B. 102:10446–10458.
Mitchell, P. 1961. Coupling of phosphorylation to electron and proton
transfer by a chemi-osmotic type of mechanism. Nature. 191:144–148.
Nagle, J. F. 1987. Theory of passive proton conductance in lipid bilayers.
J. Bioener. Biomemb. 19:413–426.
Nagle, J. F., and H. J. Morowitz. 1978. Molecular mechanisms for proton
transport in membranes. Proc. Natl. Acad. Sci. U. S. A. 75:298–302.
Okamura, M. Y., M. L. Paddock, M. S. Graige, and G. Feher. 2000. Proton
and electron transfer in bacterial reaction centers. Biophys. Biochim.
Acta. 1458:148–163.
Paula, S., G. Volkov, N. van Hoek, T. H. Haines, and D. W. Deamer. 1997.
Permeation of protons, potassium ions, and small polar molecules
through phospholipid bilayers as a function of membrane thickness.
Biophys. J. 70:339–348.
Phillips, L. R., C. D. Cole, R. J. Hendershot, M. Cotten, T. A. Cross, and
D. D. Busath. 1999. Noncontact dipole effects on channel permeation:
III. Anomalous proton conductance effects in gramicidin. Biophys. J.
77:2492–2501.
Pome`s, R. 1999. Theoretical studies of the Grotthuss mechanism in bio-
logical proton wires. Isr. J. Chem. 39:387–395.
Pome`s, R., and B. Roux. 1995. Quantum effects on the structure and
energy of a protonated linear chain of hydrogen-bonded water mole-
cules. Chem. Phys. Lett. 234:416–424.
Pome`s, R., and B. Roux. 1996a. Theoretical study of H translocation
along a model proton wire. J. Phys. Chem. 100:2519–2527.
Pome`s, R., and B. Roux. 1996b. Structure and dynamics of a proton wire:
a theoretical study of H translocation along the single-file water chain
in the gramicidin A channel. Biophys. J. 71:19–39.
Pome`s, R., and B. Roux. 1998. Free energy profiles for H conduction
along hydrogen-bonded chains of water molecules. Biophys. J. 75:
33–40.
Quigley, E. P., P. Quigley, D. S. Crumrine, and S. Cukierman. 1999. The
conduction of protons in different stereoisomers of dioxolane-linked
gramicidin A channels. Biophys. J. 77:2479–2491.
Roux, B. 1995. The calculation of the potential of mean force using
computer simulations. Comp. Phys. Commun. 91:275–282.
Roux, B., and M. Karplus. 1994. Molecular dynamics simulations of the
gramicidin channel. Annu. Rev. Biophys. Biomol. Struct. 23:731–761.
Sadeghi, R. R., and H. P. Cheng. 1999. The dynamics of proton transfer in
a water chain. J. Chem. Phys. 111:2086–2094.
Sagnella, D. E., K. Laasonen, and M. L. Klein. 1996. Ab initio molecular
dynamics study of proton transfer in a polyglycine analog of the ion
channel gramicidin A. Biophys. J. 71:1172–1178.
Saraste, M. 1999. Oxidative phosphorylation at the fin de sie`cle. Science.
283:1488–1493.
Schmitt, U. W., and G. A. Voth. 1999. The computer simulation of proton
transport in water. J. Chem. Phys. 111:9361–9381.
Schumaker, M. F., R. Pome`s, and B. Roux. 2000. A framework model for
single proton conductance through gramicidin. Biophys. J. 79:
2840–2857.
Schumaker, M. F., R. Pome`s, and B. Roux. 2001. A combined molecular
dynamics and diffusion model of single proton conduction through
gramicidin. Biophys. J. 80:12–30.
Sjogren, T., M. Svensson-Ek, J. Hajdu, and P. Brzezinski. 2000. Proton-
coupled structural changes upon binding of carbon monoxide to cyto-
Proton Relay in Gramicidin 2315
Biophysical Journal 82(5) 2304–2316

chrome cd(1): a combined flash photolysis and X-ray crystallography
study. Biochemistry. 39:10967–10974.
Stillinger, F. H. 1979. Dynamics and ensemble averages for the polariza-
tion models of molecular interactions. J. Chem. Phys. 71:1647–1651.
Stillinger, F. H., and C. W. David. 1978. Polarization model for water and
its ionic dissociation products. J. Chem. Phys. 69:1473–1484.
Tian, F., and T. A. Cross. 1999. Cation transport: an example of structural
based selectivity. J. Mol. Biol. 285:1993–2003.
Tuckerman, M., K. Laasonen, M. Sprik, and M. Parrinello. 1995. Ab initio
molecular dynamics simulation of the solvation and transport of H3O

and OH ions in water. J. Phys. Chem. 99:5749–5752.
Vuillemier, R., and D. Borgis. 1998. Quantum dynamics of an excess
proton in water using an empirical valence-bond Hamiltonian. J. Phys.
Chem. 102:4261–4264.
Vuillemier, R., and D. Borgis. 1999. Transport and spectroscopy of the
hydrated proton: a molecular dynamics study. J. Chem. Phys. 111:
4251–4266.
Weber, T. A., and F. H. Stillinger. 1982. Reactive collisions of H3O and
OH studied with the polarization model. J. Phys. Chem. 86:
1314–1318.
Wikstro¨m, M. 1998. Proton translocation by bacteriorhodopsin and heme-
copper oxidases. Curr. Opin. Struct. Biol. 8:480–488.
Woolf, T. B., and B. Roux. 1994. Molecular dynamics simulations of the
gramicidin channel in a phospholipid bilayer. Proc. Natl. Acad. Sci.
U. S. A. 91:11631–11635.
Zaslavsky, D., and R. B. Gennis. 2000. Proton pumping by cytochrome
oxidase: progress, problems and postulates. Biochim. Biophys. Acta.
1458:164–179.
2316 Pome`s and Roux
Biophysical Journal 82(5) 2304–2316