Free Energy Profiles for H Conduction along Hydrogen-Bonded Chains
of Water Molecules
Re´gis Pome`s and Benoıˆt Roux
Groupe de Recherche en Transport Membranaire, De´partements de Physique et de Chimie, Universite´ de Montre´al,
Succursalle Centre-Ville, Montre´al, Que´bec H3C 3J7, Canada
ABSTRACT The molecular mechanism for proton conduction along hydrogen-bonded chains, or “proton wires,” is studied
with free energy simulations. The complete transport of a charge along a proton wire requires two complementary processes:
1) translocation of an excess proton (propagation of an ionic defect), and 2) reorientation of the hydrogen-bonded chain
(propagation of a bonding defect). The potential of mean force profile for these two steps is computed in model systems
comprising a single-file chain of nine dissociable and polarizable water molecules represented by the PM6 model of Stillinger
and co-workers. Results of molecular dynamics simulations with umbrella sampling indicate that the unprotonated chain is
preferably polarized, and that the inversion of its total dipole moment involves an activation free energy of 8 kcal/mol. In
contrast, the rapid translocation of an excess H across a chain extending between two spherical solvent droplets is an
activationless process. These results suggest that the propagation of a bonding defect constitutes a limiting step for the
passage of several protons along single-file chains of water molecules, whereas the ionic translocation may be fast enough
to occur within the lifetime of transient hydrogen-bonded water chains in biological membranes.
INTRODUCTION
Hydrogen-bonded chains of water molecules, such as those
observed in the structure of the photosynthetic reaction
center (Baciou and Michel, 1995), and in the lumen-side
domain of cytochrome f (Martinez et al., 1996), are thought
to be implicated in the function of energy-transducing mem-
branes. An emerging idea is that such chains could act as
proton wires by providing an effective pathway for the rapid
translocation of protons over long distances. A similar
mechanism, with transient hydrogen-bonded water chains,
has been proposed in relation to proton leakage across lipid
bilayers (Nagle, 1987; Deamer and Nichols, 1989; Marrink
et al., 1996). A characterization of proton wires at the
microscopic level is required for a better understanding of
those complex systems.
The narrow transmembrane channel formed by gramici-
din (GA) provides a remarkable model system for which the
translocation of protons along a single-file hydrogen-
bonded chain of water molecules can be characterized quan-
titatively (Akeson and Deamer, 1991). In this system, the
high rate of diffusion of H arises because the mechanism
for the translocation of the excess proton does not require
the net molecular displacement of an ionic solute or of water
molecules along the way. Instead, protons can “hop” or
transfer along each hydrogen bond of a preexisting network
of hydrogen bonds, following a Grotthuss mechanism in-
volving structural rearrangement rather than molecular dif-
fusion (Agmon, 1995). The Grotthuss mechanism can be
described as the combination of two complementary steps
(Nagle and Tristam-Nagle, 1983). The conduction of sev-
eral protons in one direction under the influence of a proton
electrochemical potential necessitates both 1) the transloca-
tion of proton, or propagation of an “ionic defect,” and 2)
the reorientation of the chain, or propagation of a “turning
defect.” This is because, as shown in Scheme 1, the trans-
location, or “hop,” of an excess proton (II) along a chain of
suitably oriented hydrogen bonds (I) leaves the hydrogen
bonds in the opposite orientation (III). Thus the transloca-
tion of the next proton in the same direction requires the
preliminary reorientation of the chain from III to I. The
latter process is referred to as “turn” because it involves the
reorientation of each water molecule (or other proton donor-
acceptor group directly involved in the wire) in the chain.
Because this process necessitates breaking and reforming
hydrogen bonds along the chain, it is also called “propaga-
tion of a bonding defect.”
Even though putative mechanisms of biological proton
conduction by such hydrogen-bonded networks, or proton
wires, were formulated two decades ago (Nagle and
Morowitz, 1978), the molecular basis for the fast transloca-
tion of protons along hydrogen-bonded chains is not fully
understood. Experimentally, this is due to the difficulty of
detecting what is essentially a transient phenomenon. On
the theoretical front, however, several recent studies have
begun to provide meaningful insight into the mechanism of
proton transport in hydrogen-bonded water networks. Ab
initio studies of the dynamics of proton transfer in water
have revealed how the migration process of H in solution
arises from the structural reorganization of the hydrogen-
bonded network involving neighboring water molecules
(Ando and Hynes, 1994; Tuckerman et al., 1995). Molecu-
lar dynamics studies of single-file water chains in vacuo
using an empirical force field and discretized Feynman path
Received for publication 8 January 1998 and in final form 15 April 1998.
Address reprint requests to Dr. Re´gis Pome`s, Theoretical Biology and
Biophysics Group, T-10, Mail Stop K-710, Los Alamos National Labora-
tory, Los Alamos, NM 87545. Tel.: 505-665-9930; Fax: 505-665-3493;
E-mail: regis@lanl.gov.
© 1998 by the Biophysical Society
0006-3495/98/07/33/08 $2.00
33Biophysical Journal Volume 75 July 1998 33–40

integrals have helped to characterize the nature of the fluc-
tuations governing proton translocation and the interplay
between the flexibility of the water chain and the contribu-
tion of quantum effects (Pome`s and Roux, 1995, 1996a).
The approach was also used to determine the influence of
the GA channel on the mechanism of proton translocation
(Pome`s and Roux, 1996b). The latter simulations suggested
that ionic translocation across much of the length of the pore
can occur within the picosecond time range, a time scale
supported by independent simulations of proton transloca-
tion in a polyglycin analog of the GA channel performed at
the ab initio level (Sagnella et al., 1996).
Although such studies have unveiled important properties
of the proton wire, the limitations of proton conduction
along water chains and other hydrogen-bonded chains re-
main unclear (see the discussion by DeCoursey and Cherny,
1997). This arises in part because the focus has been placed
primarily on proton transport itself, and not on the reorien-
tation of the unprotonated water chain. Yet the slow rate of
water rotational reorientation in model ion channels (Breed
et al., 1996), the relatively rare occurrence of dipole reori-
entation of a water chain forced to extend across the lipid
bilayer (Marrink et al., 1996), and the persistence of turning
defects over several hundred picoseconds in the GA channel
(Pome`s and Roux, 1996b) suggest that a comparatively
large activation energy may be involved in the reorientation
step, in both polar and nonpolar environments. This could
have important implications, particularly in the case of
transient, nonpolar water pores. Because water wires ex-
tending across the lipid bilayer are thermodynamically un-
stable, they form infrequently, and their lifetime is probably
limited to a few picoseconds (Marrink et al., 1996). Thus the
possible role of such hydrogen-bonded chains in proton
leakage depends on their kinetic competence in translocat-
ing protons before breaking up. The properties underlying
proton conduction clearly require further characterization.
As a first step, the present work extends the use of compu-
tational techniques to calculate the free energy profiles for
the Grotthuss mechanism intrinsic to a hydrogen-bonded
chain of water molecules in a model nonpolar pore.
In this paper we characterize the translocation and reori-
entation processes separately, using molecular dynamics
umbrella sampling simulations. We calculate the potential
of mean force (PMF) profiles for charge migration and for
the reorientation of the water chain in two model systems of
proton wires, and contrast the mechanisms involved. The
first model consists of an isolated linear chain of nine water
molecules, whereas in the second model, spherical droplets
of 25 water molecules are added to each end of the single
file to mimic the influence of a polar environment at either
side of the membrane. The center of charge is used as
reaction coordinate for the translocation of the excess pro-
ton along the single file, whereas the reorientation is char-
acterized by the inversion of the total dipole moment of the
water chain. Results indicate that in the presence of solvent
droplets, the protonic translocation proceeds spontaneously
along the entire length of the chain at 300 K. In contrast, the
reorientation of the unprotonated wire involves a substantial
free energy barrier of 7–8 kcal/mol. The respective mech-
anisms are discussed, and possible implications for the
mechanism of proton conduction across the biological
membrane are proposed.
METHODS
Two models were constructed: a linear chain of nine water molecules in
vacuo, and a “dumbbell” model of 59 water molecules. Both were used to
perform classical molecular dynamics simulations, successively in the
unprotonated chain and in the presence of an excess proton. The PMF for
the reorientation of the dipole of the single-file region in the unprotonated
chain and the potential of mean force for the displacement of the center of
charge along the axis of the protonated chain were calculated for each of
the two models by umbrella sampling simulations (Allen and Tildesley,
1987).
Potential energy function
The potential energy function used to model the single-file water chain, the
Polarization Model, version PM6, was designed by Stillinger and co-
workers to reproduce the properties of water and its ionic products (Weber
and Stillinger, 1982). The PM6 model consists of O2 and H ions. In
addition to pairwise-additive terms, the electronic polarizability is ac-
counted for by a many-body function. Because it is dissociable and
polarizable, the PM6 model offers a reasonably realistic method for the
simulation of proton transfer processes in water clusters and a tractable
alternative to costly ab initio simulations (Pome`s and Roux, 1996a). For a
complete description of the potential energy function, the reader is referred
to the original papers (Stillinger and David, 1978; Stillinger, 1979; Weber
and Stillinger, 1982); for a description of its use in conjunction with
molecular dynamics simulations of proton wires, see our earlier papers
(Pome`s and Roux, 1995, 1996a,b).
The first model (Model I) of the proton wire consists of a linear chain
of nine water molecules. To preserve the linearity and connectivity of the
hydrogen-bonded chain, the oxygen atoms were subjected to two con-
straints: 1) cylindrical restoring force on the oxygen atoms beyond 0.75 Å
from the chain’s axis (z axis) to maintain linearity; 2) harmonic restoring
force on the end points outside of the range z  10.00 Å to avoid
evaporation of the water molecules. The force constants of both constrain-
ing potentials were chosen as 10 kcal/mol/Å2. The above constraints mimic
Scheme 1.
34 Biophysical Journal Volume 75 July 1998

the influence of a hypothetical, perfectly cylindrical hydrophobic pore
accommodating a single file of water molecules.
The second model (Model II) was designed to characterize the influence
of bulk water lying at either end of the single file region on the H
translocation and water reorientation processes. To construct the second
model, the chain of nine PM6 water molecules was capped by spherical
droplets of 25 water molecules. These spherical droplets are in tangential
contact with the extrema of the linear chain, thus constituting a “dumb-
bell.” The TIP3P potential energy function (Jorgensen et al., 1983) was
used for the water droplets. The computational advantage from using
unpolarizable and undissociable potential energy functions such as TIP3P
for the bulk water regions arises in large systems because the PM6 model
is significantly more costly in terms of computational time, because of the
iterative calculation of the many-body polarization term.
To preserve the dumbbell geometry, the water droplets’ O atoms were
subjected to harmonic restoring forces acting beyond 5.5 Å from the
droplets’ centers, which were located, respectively, at z  17.2 Å.
Because of the constraints imposed, the interfaces between the droplets and
the single file are well defined spatially: the nine PM6 water molecules
remain in the single-file region throughout the simulations, and the two
water models cannot mix or interconvert.
The optimization of the TIP3P-PM6 water-water interactions was re-
quired for consistency at the interface between the linear chain and the
spherical droplet. The dissociation energy of two water molecules in vacuo
(H2O-HOH) into two isolated water molecules, and the dissociation of a
protonated water dimer (H2O-H
-OH2) into H2O and H3O
 were exam-
ined. The optimized geometries and energies are given in Table 1. The
PM6 O atoms were assigned Lennard-Jones parameters identical to those
of the TIP3P O atoms for the calculation of PM6-TIP3P interactions. The
Lennard-Jones parameters of PM6 H atoms were   0.0498 kcal/mol and
Rmin  1.84 Å. Finally, the Lennard-Jones parameter Rmin for the O(PM6)-
H(TIP3P) pairwise potential term was taken as 2.43 Å to match the
water-water hydrogen-bonding energies of the two possible isomers, do-
nor-acceptor and acceptor-donor. At 8.6 kcal/mol, this energy is com-
parable to the hydrogen bond energy of the TIP3P dimer, 7.0 kcal/mol.
Further optimization of the system is not an easy task, as the two potential
functions have radically different forms. Although discrepancies remain in
both the energy and the geometry of the hybrid hydrogen bonds, the
approximation is sufficient to conserve the tetrahedral coordination of PM6
water molecules in a TIP3P environment, and to allow the hydrogen bonds
to reorganize without favoring one type of TIP3P-PM6 hydrogen bond
over the other (donor-acceptor versus acceptor-donor). The energy of
dissociation of the protonated dimer is 28% smaller in the hybrid model
than in the all-PM6 O2H5
 (Table 1). Because the hybrid model is not
intended for the study of proton diffusion in the droplets, and because
TIP3P-PM6 proton transfer is precluded by construction, the preference for
proton hydration in the PM6 part of the dumbbell model should not affect
the study of proton translocation in the single-file region.
Molecular dynamics simulations
After inclusion of the PM6 potential energy function, version 22 of the
CHARMM program (Brooks et al., 1983; Mackerell et al., 1998) was used
to generate the molecular dynamics trajectories from which the PMF
profiles were calculated. The Langevin equations of motion were propa-
gated at 300 K with a friction coefficient of 1 ps1 and a time step of 0.5
fs. Each picosecond of dynamics required 30 s of CPU for the nonameric
chain, and 6 min of CPU for the dumbbell, on a R-4000 SGI workstation.
The PMF profiles were calculated as follows. The configuration of the
system was recorded every 5 fs. The projection of the total dipole moment
of the wire (PM6 waters only) on the z axis of the linear hydrogen-bonded
chain was calculated for each configuration at time t as
zt qH 
i1
NH
zHit qO 
j1
NO
zOjt (1)
in units of e  Å, where NH and NO are the total numbers of H and O nuclei,
and qH and qO are the formal charges of H and O (respectively, 1 and 2e
in the PM6 model). In the unprotonated chain (O9H18), this expression
yields the z component of the total dipole moment, whereas in the proton-
ated chain (O9H19
 ), it corresponds to the projection of the center of charge
along the z axis. The PMF was calculated from the equilibrium probability
distributions of z obtained from the simulations.
A biasing potential energy function was imposed on z to force the
sampling over the large free energy barrier for reorientation of the unpro-
tonated chain. This “umbrella sampling” calculation (Allen and Tildesley,
1987) was performed by means of harmonic biasing functions of the form
Vi(z)  (1/2)ki(z  zi)2 imposed on z. All ki coefficients were set to
20 kcal  mol1  e2  Å2, whereas zi varied in increments of 0.5, from
6.0 to 0.5 e  Å in the nonamer, and from 7.0 to 0.5 e  Å in the
dumbbell. Each window of the umbrella simulation consisted of a 60-ps
run with 10 ps of equilibration and 50 ps of data collection. The simulation
time required to build the PMF for the wire’s reorientation totaled 900 ps
for the isolated nonamer, and 950 ps for the dumbbell. The construction of
the PMF profiles from the data obtained from the various umbrella simu-
lations, or “unbiasing,” was performed with the WHAM algorithm (Kumar
et al., 1992; Roux, 1995). The same calculation was repeated for an
isolated chain of nine TIP3P water molecules to gauge the influence of
polarization of the force field on the reorientation mechanism (total sim-
ulation time of 3.0 ns).
In the simulations of the protonated wire, umbrella sampling was not
required, because there were no substantial free energy barriers opposing
the translocation of H. Consequently, the PMF for proton mobility was
computed from long unbiased simulations totaling 2.5 ns in the isolated
nonameric chain and 1 ns in the dumbbell system.
RESULTS
Propagation of a bonding defect
At equilibrium, the unprotonated chain of nine PM6 water
molecules is arranged preferentially in one of two “polar-
ized” configurations for which the total dipole moment is
aligned with the chain axis. Moreover, each configuration
corresponds to an oriented donor-acceptor pattern in the
hydrogen bonding between the adjacent water molecules of
the wire (Fig. 1). This can be rationalized in simple elec-
trostatic terms from the maximization of favorable interac-
tions between the molecular dipoles of each pair of water
molecules and from the formation of a continuous hydro-
gen-bonded chain. The reorientation of the chain is sequen-
tial, as illustrated in Fig. 1: intermediate configurations
involve the successive reorientation of water molecules, 1,
2, 3, etc., rather than the flip of individual water molecules
located randomly in the chain. Thus the interconversion
between the two polarized configurations involves the mi-
TABLE 1 PM6-TIP3P interactions
Dimer
Ebinding
(kcal/mol)
O . . . O
(Å)
O . . . H
(Å)
H2O. . . . . .HOH
TIP3P PM6 8.56 3.01 2.04
PM6 TIP3P 8.56 2.85 1.86
TIP3P TIP3P 7.04 2.75 1.76
PM6 PM6 5.55 2.90 1.91
H2O. . . . . .HOH2

PM6 PM6 32.78 2.47 1.24
TIP3P PM6 23.77 2.99 1.94
Pome`s and Roux Free Energy Profiles for H Conduction 35

gration of a single hydrogen-bonding defect from end to
end, consistently with the idealized picture of the Grotthuss
mechanism.
The completely oriented chain corresponds to deep wells
at z  7.15 e  Å in the PMF for the reorientation of the
unprotonated nonameric chain in vacuo (Fig. 2). In contrast,
the state corresponding to two oriented “half-chains” where
the central water molecule projects a dipole moment of zero
along the chain axis, and the dipoles of the outer four
molecules on each side of the central water point either
toward it or away from it (Fig. 1), corresponds to a free
energy barrier of 7.5 kcal/mol. Dipole-dipole and hydrogen-
bonding interactions of the water molecules in the linear
hydrogen-bonded chain constitute the dominant forces that
favor the orientation of the chain and oppose the passage
from one stable configuration to the other.
The free energy profile for reorienting a TIP3P water
chain is included in Fig. 2 for comparison. Although the
profile is qualitatively similar to that obtained with the PM6
water model, at 5.2 kcal/mol the barrier is smaller by 30%.
The inclusion of electronic polarizability has a significant
effect on the reorientation process: the structure of the
(unpolarizable) TIP3P water chain near z  0 exhibits less
order than its PM6 counterpart (not shown), and the free
energy penalty for diminishing the magnitude of the chain’s
dipole moment is higher when polarization is included.
Nevertheless, the TIP3P model, by favoring strongly the
totally oriented configurations over partially oriented or
random configurations, reproduces the essential features of
the reorientation process obtained with the PM6 model. In
particular, both water models indicate the existence of a
significant free energy barrier.
The PMF profile for the reorientation of the dumbbell
model, where solvent droplets were added to the extremities
of the nonameric chain, is also shown in Fig. 2. The free
energy for flipping the total dipole of the chain is very
similar to that of the isolated chain, and involves a slightly
larger barrier of 8.3 kcal/mol. Although the presence of a
polar medium capping the model pore further stabilizes the
polarized conformations of the chain, the influence of the
bulk water droplets at either end of the single file region is
not dramatic. Compared to the isolated water nonamer, the
oriented chain has a larger dipole component along the
channel axis (7.75 versus 7.15 e  Å), and the propaga-
tion of the turning defect is easier throughout the center of
the chain: the slope of the barrier decreases notably after the
one or two water molecules near the end of the single file
have reoriented. Whereas the free energy cost for rotation of
the waters at the end of the single file is largest, the relative
flatness of the PMF profile near the top of the barrier
reflects the somewhat collective nature of the reorientation
of the central water molecules in the chain.
Propagation of an ionic defect
As described in more detail elsewhere (Pome`s and Roux,
1995, 1996a), the structure of linear chains of water mole-
cules is strongly affected by the inclusion of an excess
proton. The excess charge of H induces short, strong
water-water hydrogen bonds, along which the potential en-
ergy profile for proton transfer involves small or nonexist-
ent barriers. Fluctuations in the polar environment of the
proton modulates the asymmetry of the well(s) and result in
variations in the O-O separations, thereby leading to the
propagation of the proton along the system. The proton
transfer itself is thus determined by thermal fluctuations of
the geometry of the hydrogen bonds in the water chain. An
FIGURE 1 Snapshots of the PM6 water nonamer illustrating the sequen-
tial reorientation of the hydrogen-bonded chain. Hydrogen bonding is
indicated by dashed lines. From left to right, z  7.0 (oriented chain),
3.5 (reorientation of the second water from bottom), 2.0 (reorientation
of the third water from bottom), and 0.0 e  Å (the fifth water reorients, two
oriented half-chains).
FIGURE 2 Potential of mean force for the reorientation of the unproto-
nated water chain: ——, PM6 water nonamer; – – –, dumbbell model; ,
TIP3P water nonamer. Note that in the latter case, z was scaled by 0.417,
the fractional charge of H atoms in the TIP3P potential, for consistency
with the magnitude of z obtained with the PM6 model, in which the
formal charge of H is 1 (see Eq. 1).
36 Biophysical Journal Volume 75 July 1998

analogous process involving the rearrangement of water
molecules around (OH3)
 and (O2H5)
 ionic species has
been proposed as a mechanism for proton diffusion in bulk
water that may be called “structural diffusion” (Tuckerman
et al., 1995).
In the nonameric water chain in vacuo, the excess proton
is located on average near the center of the system, and
thermal fluctuations cause its diffusion along the four cen-
tral hydrogen bonds of the chain. The excess proton does
not wander off to the extremities of the hydrogen-bonded
chain, where the charge would be more poorly solvated than
in the center. The PMF for proton mobility in the isolated
nonameric water chain is shown in Fig. 3. The free energy
profile for the motion of the center of charge along the
chain’s axis is a single well centered about the origin, which
corresponds on average to a central hydronium ion and
to two oriented half-chains whose dipoles point away from
the excess charge. This free energy well is fit well by a
quadratic function f (z)  (1/2)kz2, where k  0.55
kcal  e2  Å2.
After the addition of solvent droplets, the free energy
profile changes to a broad, flat-bottomed well (Fig. 3). This
is because the presence of the excess proton near the ex-
tremities of the single file is now electrostatically stabilized
by the presence of the polar medium provided by the water
droplets. As a result, the translocation of H takes place
spontaneously over the entire length of the wire under the
influence of thermal fluctuations. This process is very fast
and occurs within the time scale of a picosecond. It should
be noted that because the Langevin equation was used to
propagate the equations of motion, however, it is not pos-
sible to extract more precise temporal estimates from these
simulations.
As a collective reaction coordinate, the location of the
center of charge (z) reflects not only the location of the
excess proton in the (O9H19)
 chain, but also the configu-
ration of all of the water molecules in the single file, thereby
taking into account changes in the orientation of individual
water molecules as the charge moves along the chain. It is
useful to correlate that reaction coordinate to the location of
the excess proton along the hydrogen-bonded chain. The
latter is approximated by the following function:
OO 


i1
8
wi

1

i1
8
wii 4.5 (2)
where wi  e(ri  )/(1  e(ri  )),   30 Å1, 
2.58 Å, and ri is the distance separating O atoms i and i 
1 in the linear chain. wi is a smooth step-like function that
gives most weight to strong hydrogen bonds (ri
2.5 Å),
less to intermediate-strength hydrogen bonds (2.5
ri

2.7 Å), and least to weaker hydrogen bonds (ri  2.7 Å).
Thus the function OO follows the location of the strongest
bond(s) in the chain. Such a reaction coordinate makes use
of the strong coupling between short O-O separations and
the presence of H: for each configuration, the shortest O-O
“bond” in the wire is assigned the strongest weight because
it is the one most likely to host the excess proton. This is
equivalent to following the midpoint of a hypothetical
O2H5
 cluster over the course of the simulation. A more
general formulation extendable to two- and three-dimen-
sional hydrogen-bonded networks would be obtained by
replacing the dimensionless index i in the above expression
for OO by the Cartesian coordinates of the midpoint of each
hydrogen bond. In simulations of the GA channel, this
approach was shown to constitute a good approximation of
the location of the excess charge in the proton wire (Pome`s
and Roux, 1996b). As seen from Fig. 4, the two reaction
coordinates z and OO are strongly correlated throughout
the course of the simulation, and the excess proton does
indeed travel from end to end of the hydrogen-bonded
chain: at OO  3.5, water molecules 1–2 and 8–9,
respectively, have the properties of a protonated water
FIGURE 3 Potential of mean force for the translocation of an excess
proton. – – –, Isolated (O9H19)
; ——, dumbbell model.
FIGURE 4 Comparison of the two collective reaction coordinates for
proton translocation, z and OO, obtained from the dumbbell simulation.
There is a linear relationship between the location of the center of charge
on the z axis, z, and the index of the strongest hydrogen bond, given by
OO (see text).
Pome`s and Roux Free Energy Profiles for H Conduction 37

dimer, O2H5
. Because it incorporates the overall configu-
ration of the hydrogen-bonded chain, z is consistent with
the reaction coordinate chosen in the study of the unproto-
nated chain and constitutes an appropriate collective reac-
tion coordinate to follow proton translocation.
DISCUSSION
The mechanism of proton translocation is a semicollective
process involving thermal fluctuations of water-water hy-
drogen bond lengths. As described in detail elsewhere
(Pome`s and Roux, 1995, 1996a), these fluctuations govern
the migration of the excess proton along a tight core of
water molecules in the hydrogen-bonded chain. These stud-
ies also indicated that because the thermal fluctuations
governing the transport of H along hydrogen-bonded wa-
ter chains involve predominantly classical degrees of free-
dom, the classical limit is appropriate. In an isolated linear
water cluster, the tight protonated core is confined to the
center of the chain. However, under the influence of a polar
medium at the ends of the single-file chain, at 300 K the
migration of an excess proton in a nonpolar pore can occur
spontaneously from end to end of the wire with thermal
fluctuations. This result is in striking contrast to the large
activation energy required to reorient the unprotonated
chain.
In the absence of an excess proton, the thermodynamic
preference for oriented hydrogen-bonded chains is an in-
trinsic property of water in a nonpolar pore. There are, in
principle, many ways to orient the individual water mole-
cules so that the net dipole projection is identically zero. For
instance, any configuration in which the z components of
adjacent pairs of water molecules cancel each other results
in z  0. This would tend to favor such configurations
entropically. Yet in all such cases with a net dipole of zero,
the sum of all unfavorable molecular dipole-dipole interac-
tions is maximal compared to the oriented configurations. In
addition, the enthalpic cost of breaking a hydrogen bond in
the single file is large enough that only one defect is
observed at a time. The presence of a free energy barrier for
the chain reorientation is thus due primarily to enthalpic
factors. In previous computer simulations of single-file wa-
ter chains in GA (Chiu et al., 1991; Roux and Karplus,
1994, and references therein), in the lipid bilayer (Marrink
et al., 1996), and in model pores (Breed et al., 1996), the
hydrogen-bonded chain was often observed to be oriented
rather than disordered. The slow rotational relaxation of
water in model hydrophobic channels has been studied in
detail (Breed et al., 1996). The collective nature of the
reorientation process, after its initiation by turning of the
water molecules at the end points, was noted by Marrink et
al. (1996). Furthermore, the reduced barrier obtained with
the TIP3P model indicates that electronic polarizability is
important in modeling of the dynamics of single-file water
chains, in agreement with other studies (Duca and Jordan,
1997).
A motivation for this study arose from relatively long
(100–400 ps) molecular dynamics simulations of the GA
channel (Pome`s and Roux, 1996b), which suggested that the
overall translocation of protons is limited by slow rear-
rangements of the hydrogen-bond defects in the hydrogen-
bonded network of the proton wire. These rearrangements
involve changes in the relative orientation of individual
water molecules in the chain and constitute the primary
event in the propagation of a bonding defect. Because they
occur on a typical time scale of 100 ps in the GA channel,
they appear to be essentially independent of the motion of
the proton itself, which takes place on a much shorter time
scale (less than 1 ps) (Pome`s and Roux, 1996b; Sagnella et
al., 1996). Thus the time scales of the physical events
responsible for the propagation of the ionic and bonding
defects differ very significantly from each other. If the water
molecules in the pore do not form a continuous hydrogen-
bonded chain, the limitation for proton translocation is
given by the making and breaking of hydrogen bonds. Yet
the present simulations suggest that it is costly to interrupt
the hydrogen-bonded water chain in a nonpolar environ-
ment. The scarcity of hydrogen-bonding partners (at most
two per water molecule) favors the presence of hydrogen
bonds between adjacent water molecules. This factor, to-
gether with the preference for an oriented hydrogen-bonded
chain, contributes to a fast ionic translocation while hinder-
ing the propagation of the bonding defect.
This separation of time scales has direct implications for
the Grotthuss mechanism of proton transport in transient,
nonpolar water pores. To explain proton leakage across
biological membranes, a plausible mechanism has been
proposed that necessitates a very fast translocation of pro-
tons through transient water chains extending across the
lipid bilayer (Nagle, 1987; Deamer and Nichols, 1989;
Marrink et al., 1996, and references therein). This is because
the free energy of the formation of a single file of water
molecules in the adverse, hydrophobic core of the lipid
bilayer is largely unfavorable. Accordingly, computer sim-
ulations of water in an explicit lipid bilayer suggest that
such chains form rarely and have a very short lifetime,
estimated at a few picoseconds (Marrink et al., 1996). Based
on their observation of a single occurrence of turning defect
translocation during their 320-ps simulation of a water chain
constrained to extend across the lipid bilayer, Marrink et al.
(1996) conclude that the reorientation step is likely to be too
infrequent to occur during the lifetime of the water pore,
which rules out the passage of several protons. This argu-
ment is supported by the magnitude of the energy barrier
opposing the propagation of a turning defect, i.e., 10 times
the thermal energy, obtained here. Furthermore, provided
that the diffusion of H to the pore entrance is not rate
limiting, the concentration of these water pores would be
sufficient to explain measured conductance rates of protons
by the membrane only if it is assumed that a proton is
translocated through each of them before their collapse.
Assuming an activation control for the ionic translocation,
Marrink et al. (1996) estimate that the rate of passage of
38 Biophysical Journal Volume 75 July 1998

protons is on the order of one per nanosecond, which leads
them to question the validity of such a mechanism. How-
ever, the activationless nature of the PMF for proton trans-
location and the high rate of ionic diffusion observed in this
work suggest that transient water chains would be kineti-
cally competent to permit proton leakage. Recent measure-
ments of the permeability of the phospholipid bilayer to
protons were consistent with the transient pore mechanism
for proton translocation in lipids with acyl chains shorter
than 20–22 carbon atoms, corresponding to an30-Å-thick
hydrophobic region (Paula et al., 1996). The present study
provides elements to reconcile the transient hydrogen-
bonded chain mechanism with a detailed molecular picture
of water wires in nonpolar channels.
The model systems studied here represent proton wires
that are not subjected to any electrostatic interactions with
the wall of the pore. As such, they constitute an idealized
model of hydrophobic water channels. In the presence of
local electrostatic interactions between the water wire and
hydrogen-bonding or polar groups in a biological channel,
one would expect particular configurations of the water
chain to be stabilized, thus modifying the free energy pro-
file. For example, similar to the results obtained for other
cations (Woolf and Roux, 1997), PMF profiles may reveal
“hydrated-proton binding sites” in the GA channel. Like-
wise, polar interactions involving the single-file region,
while favoring pore hydration, are also likely to affect the
free energy profile for the wire’s reorientation after the
passage of a proton. For instance, for an efficient proton
transport mechanism it is possible that polar transmembrane
proton conductors such as the GA channel facilitate the
reorientation process relative to nonpolar water pores
through the stabilization of intermediate configurations.
Thus it may be expected that the balance of water-water and
water-channel interactions modulates the relative ease of
ionic and bonding translocations.
CONCLUSIONS
In the above study we investigated the dynamics and PMF
profiles for the translocation of an excess proton and for the
reorientation of a linear chain of nine water molecules,
successively in vacuo and with the addition of two droplets
of water (dumbbell) to gauge the influence of bulk solvent
at either end of the single file. The translocation of an excess
proton along the single file region is confined to the center
of the isolated chain for energetic reasons, whereas in the
dumbbell model the stabilization provided by the water
droplets allows the proton to readily migrate from one end
of the chain to the other. This process involves no signifi-
cant free energy barrier and thus occurs spontaneously with
thermal energy.
In contrast, the PMF for the reorientation of the dipoles of
the unprotonated hydrogen-bonded water chain involves a
free energy barrier of 8 kcal/mol. This barrier does not
change significantly upon the inclusion of solvent droplets
and can be attributed to the intrinsic enthalpy for unfavor-
able dipole-dipole interactions in the chain. Thus the present
study suggests that the reorientation process, not the passage
of a proton itself, constitutes the limiting step for the pas-
sage of several protons through water chains in nonpolar
pores.
The results support kinetic competence of a hopping
mechanism for the proton conductance by nonpolar water
pores such as those proposed to form transiently in lipid
bilayers (Nagle, 1987; Marrink et al., 1996; Paula et al.,
1996). In such systems, the lifetime of the oriented hydro-
gen-bonded chain would be long enough to allow the rapid
passage of one proton, though too short for the reorientation
to take place. Although more detailed studies are needed to
further our understanding of the factors governing proton
translocation across biological membranes, this study dem-
onstrates the feasibility of free energy simulations with
collective reaction coordinates and their usefulness in the
study of the detailed molecular properties of hydrogen-
bonded chains relevant to both steps of the Grotthuss
mechanism.
This work was supported in part by a grant from the Medical Research
Council of Canada and by the U.S. Department of Energy through the Los
Alamos National Laboratory LDRD-CD grant for Bioremediation. BR is a
research fellow of the Fonds de la Recherche Scientifique du Que´bec.
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.
Allen, M. P., and D. J. Tildesley. 1987. Computer Simulations of Liquids.
Clarendon Press, Oxford.
Ando, K., and J. T. Hynes. 1994. Ionization of acids in water. In Structure
and Reactivity in Aqueous Solution. C. J. Cramer and D. G. Truhlar,
editors. ACS Books, Washington, DC. 143–153.
Baciou, L., and H. Michel. 1995. Interruption of the water chain in the
reaction center from Rb. sphaeroides reduces the rates of the proton
uptake and of the second electron transfer to QB. Biochemistry. 34:
7967–7972.
Breed, J., R. Sankararamakrishnan, I. D. Kerr, and M. S. P. Sansom. 1996.
Molecular dynamics simulations of water within models of ion channels.
Biophys. J. 70:1643–1661.
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.
Chiu, S.-W., S. Subramaniam, E. Jacobsson, and J. A. McCammon. 1991.
Water and polypeptide conformations in the gramicidin channel. Bio-
phys. J. 56:253–261.
Deamer, D. W., and J. W. Nichols. 1989. Proton flux mechanisms in model
and biological membranes. J. Membr. Biol. 107:91–103.
DeCoursey, T. E., and V. V. Cherny. 1997. Deuterium isotope effects on
permeation and gating of proton channels in rat alveolar epithelium.
J. Gen. Physiol. 109:415–434.
Duca, K. A., and P. C. Jordan. 1997. Ion-water and water-water interac-
tions in a gramicidinlike channel: effects due to group polarizability and
backbone flexibility. Biophys. Chem. 65:123–141.
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.
Pome`s and Roux Free Energy Profiles for H Conduction 39

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. 1. The method. J. Comp. Chem.
13:1011–1021.
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, III, 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.
Martinez, S. E., D. Huang, M. Ponomarev, W. A. Cramer, and J. L. Smith.
1996. The heme redox center of chloroplast cytochrome-f is linked to a
buried 5-water chain. Protein Sci. 5:1081–1092.
Nagle, J. F. 1987. Theory of passive proton conductance in lipid bilayers.
J. Bioenerg. Biomembr. 19:413–426.
Nagle, J. F., and H. J. Morowitz. 1978. Molecular mechanisms for proton
transport in membranes. Proc. Natl. Acad. Sci. USA. 75:298–302.
Nagle, J. F., and S. Tristam-Nagle. 1983. Hydrogen bonded chain mech-
anisms for proton conduction and proton pumping. J. Membr. Biol.
74:1–14.
Paula, S., G. Volkov, N. van Hoek, T. H. Haines, and D. W. Deamer. 1996.
Permeation of protons, potassium ions, and small polar molecules
through phospholipid bilayers as a function of membrane thickness.
Biophys. J. 70:339–348.
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.
Roux, B. 1995. The calculation of the potential of mean force using
computer simulations. Comp. Phys. Comm. 91:275–282.
Roux, B., and M. Karplus. 1994. Molecular dynamics simulations of the
gramicidin channel. Annu. Rev. Biophys. Biomol. Struct. 23:731–761.
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.
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.
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.
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.
Woolf, T. B., and B. Roux. 1997. The binding site of sodium in the
gramicidin A channel: comparison of molecular dynamics with solid-
state NMR data. Biophys. J. 72:1930–1945.
40 Biophysical Journal Volume 75 July 1998