Pou
Biochemistry, Hospital for Sick
Children, 555 University
Avenue, Toronto, Ont., Canada
M5G 1X8
2Department of Biochemistry
University of Toronto, Toronto
Ont., Canada, M5S 1A8
3Department of Biochemistry
Weill Medical College of Cornell
University, 1300 York Avenue
New York, NY 10021, USA
The dielectric barrier of biological membranes rapid diffusion of potassium ions across cell
membranes via proteins called potassium channels
doi:10.1016/j.jmb.2004.08.036opposes the movement of ionic species. LivingIntroduction
Biological ion transport
using the control of ion permeation productively.1
Electric current across the membrane of cells results
from ion translocation via ion channels and trans-
porters. Understanding the permeation of ions or
neutral permeants is very significant for funda-
mental and practical reasons. For example, the*Corresponding authorsystems have evolved to explo
0022-2836/$ - see front matter q 2004 E
Abbreviations used: AQP, aquapo
dynamics; PMF, potential of mean f
uptake facilitator.
E-mail address of the correspond
pomes@sickkids.caAquaporins are an important class of membrane channels selective for
water and linear polyols but impermeable to ions, including protons.
Recent computational studies have revealed that the relay of protons
through the water-conduction pathway of aquaporin channels is opposed
by a substantial free energy barrier peaking at the signature NPA motifs.
Here, free-energy simulations and continuum electrostatic calculations are
combined to examine the nature and the magnitude of the contribution of
specific structural elements to proton blockage in the bacterial glycerol
uptake facilitator, GlpF. Potential of mean-force profiles for both hop and
turn steps of structural diffusion in the narrow pore are obtained for
artificial variants of the GlpF channel in which coulombic interactions
between the pore contents and conserved residues Asn68 and Asn203 at
the NPA signature motifs, Arg206 at the selectivity filter, and the peptidic
backbone of the two half-helices M3 and M7, which are arranged in head-
to-head fashion around the NPA motifs, are turned off selectively.
A comparison of these results with electrostatic energy profiles for the
translocation of a probe cation throughout the water permeation pathway
indicates that the free-energy profile for proton movement inside the
narrow pore is dominated by static effects arising from the distribution of
charged and polar groups of the channel, whereas dielectric effects
contribute primarily to opposing the access of HC to the pore mouths
(desolvation penalty). The single most effective way to abolish the free-
energy gradients opposing the movement of HC around the NPA motif is
to turn off the dipole moments of helices M3 and M7. Mutation of either of
the two NPA Asn residues to Asp compensates for charge–dipole and
dipole–dipole effects opposing the hop and turn steps of structural
diffusion, respectively, and dramatically reduces the free energy barrier
of proton translocation, suggesting that these single mutants could leak
protons.
q 2004 Elsevier Ltd. All rights reserved.
Keywords: proton translocation; Grotthuss mechanism; hydrogen bonded
network; molecular dynamics simulations; free energy calculations1Structural Biology andStructural Determinants of
Aquaporins
Nilmadhab Chakrabarti1,2, Benoıˆt Rit that property by
lsevier Ltd. All rights reserve
rin; MD, molecular
orce; GlpF, glycerol
ing author:roton Blockage in
x3 and Re´gis Pome`s1,2*
J. Mol. Biol. (2004) 343, 493–510controls many fundamental biological processes,
including electrical signaling in the nervous
system.2 The conduction of chloride ions through
chloride channels governs the electrical activity of
muscle cells and of certain neurons.3 The control of
proton translocation across biomembranes is an
essential aspect of biological energy transduction
d.

the translocation of any ion, including H and
K(bioenergetics).4 In energy-transducing mem-
branes, integral membrane proteins utilize photo-
chemical and/or redox reactions to pump hydrogen
ions against an electrochemical gradient. The
resulting chemiosmotic force is then used produc-
tively to synthesize ATP from ADP.5 One of the best
understood examples of proton pumps is bacterio-
rhodopsin, where the translocation of protons is
driven by the photo-isomerization of a retinal
group.6,7
In the last two decades, the elucidation of
structural details of gramicidin,8,9 the KcsA potass-
ium channel,2 and the ClC chloride channel,3 has
opened the way to a better understanding of the
molecular basis of biological ion transport. Together
with structural studies, computational studies are
helping to unveil the balance of fundamental
physical forces resulting in ion permeation.10 In
particular, these advances are providing micro-
scopic insight into the molecular properties of ion
channels compensating for the partial desolvation
of the ion in a low-dielectric environment. Selective
ion channels contain narrow water-filled pores in
which permeating ions are stabilized by short-range
interactions. In the potassium channel, backbone
carbonyl oxygen atoms lining the selectivity filter
provide multiple coordination sites for permeating
cations.2 Similarly, to facilitate the diffusion of alkali
metal ions through gramicidin, a bacterial peptide
that mediates the permeation of small monovalent
cations, backbone carbonyl groups directly solvate
the permeating ionvia charge–dipole interactions.11,12
In the chloride channel, main-chain amide NH
groups are involved in the stabilization of anions.3
In addition, long-range ion–channel interactions are
thought to contribute to lowering the chemical
potential of translocating ions. The electrostatic
field generated by the macro-dipoles of four tilted
a-helices was shown to play a role in the charge
selectivity of the KcsA channel by favoring the
presence of a cation in a water-filled cavity in the
center of the transmembrane region.13 Ion–dipole
interactions arising from side-chains of native and
modified gramicidin channels have been proposed to
modulate the affinity of these channels for protons.14
The study of proton translocation mechanisms
presents a particular challenge due to the reactivity
of HC in protic media. Protons can shuttle through
hydrogen-bonded networks through a relay process
often referred to as structural diffusion or the
Grotthuss mechanism. Because this process is very
fast, it is extremely difficult to capture experimen-
tally. Theoretical studies have begun to provide
meaningful insight into the molecular mechanism
mediating proton transport in systems of biological
relevance. The delocalized nature of charge in the
hydrated proton in a hydrogen bonded network,
which is intrinsic to the reactivity and the anom-
alous mobility of protons in water,15 has been
analyzed in detail in theoretical studies of proton
494transport in bulk water16–19 as well as in water
chains embedded in model non-polar channels or
nanotubes,20–25 and in the gramicidin channel.21,26,27OH . Nature has helped the evolution of channels
like AQPs, as it is essential for any living cell to
maintain ionic gradients across cellular membranes.
Sequence analysis, cryo-electron microscopy,
X-ray crystallography, and molecular dynamics
(MD) simulations have begun to clarify the struc-
tural basis of aquaporin function. The recent
determination of atomic-resolution structures of
Aqp1,33,34 GlpF,35 and AqpZ36 has revealed highly
conserved structural elements. Each monomer of
these homo-tetrameric channels consists of six
transmembrane helices defining an hour-glass
pore. The narrow part of the water-conducting
pathway is 20 A˚ long.33–35 The narrowest part of the
pore, which is referred to as the selectivity filter, is
located at the extracellular mouth and contains a
conserved Arg residue. Residues defining the
selectivity filter in GlpF are W48, G191, F200 and
R206 (F58, H182, C191, and R197 in bovine Aqp1).These studies have revealed how structural fluctu-
ations of hydrogen bonded networks give rise to
structural diffusion. In this process, small-scale
fluctuations in the position of heavy atoms result
in the long-range displacement of HC in nano-
second timescales. A notable aspect of the hydrated
proton that was shown to favor the solvation and
mobility of protons is the hydrogen bonded
coordination of protonated water molecules by
three acceptors.16,26 In bulk water and in water
chains, the rapid transfer or hop of proton between
water molecules is thought to be limited kinetically
by the comparatively slower reorganization of the
hydrogen bonded network.16,21,25 The latter process
is complementary to proton hopping and is some-
times referred to as the turn step of the hop-and-
turn Grotthuss relay mechanism. Together, these
two steps give rise to structural diffusion.22,28
Aquaporins
The high permeation rate of neutral solutes
such as water and linear polyalcohols has been
characterized in channels of the aquaporin (AQP)
superfamily.29–31 The assembly of aquaporins into
homo-tetrameric units leads to the formation of five
pores, one in each monomer and one central pore
that is impermeable to even water molecules.
Aquaporins transport billions of water molecules
per second to relieve osmotic imbalance across the
cell membrane. The mounting evidence of clinical
disorders arising from the physiological malfunc-
tion of aquaporins underlines the need for a
detailed understanding of these channels at the
molecular level. The analysis of human disease
states has confirmed that aquaporins are involved
in many different illnesses, including abnormalities
of kidney function, loss of vision, and onset of brain
edema and arsenic toxicity.32 In spite of the fast
permeation of water molecules, AQPs do not allow
C
Proton Blockage in AquaporinsTwo highly conserved loops, each containing the
signature Asn-Pro-Ala (NPA) motif, fold back into
the protein and meet in the center of the pore.

Approximately one half of each loop is non-helical
and defines a curvilinear conduction pathway,
constituted by backbone carbonyl groups that
point toward the channel interior and form hydro-
gen bonds with permeants. The other half of each
loop, known as M3 and M7 in GlpF (HB and HE in
Aqp1), is a-helical. These two helices are arranged
in head-to-head fashion at the center of the channel,
where two NPA motifs come close to each other. As
shown in Figure 1, in the GlpF channel, the macro-
dipoles of M3 and M7 helices point their positive
ends to NPA residues Asn68 and Asn203, respect-
ively, at the center of the channel.
Molecular dynamics studies of water movement
in Aqp1 and GlpF have provided meaningful
insight into the mechanism of water permeation.37,38
These studies have shown that up to nine water
molecules fit into the narrow pore of Aqp1 and
GlpF, where they assemble in single file and form a
hydrogen bonded chain (Figure 1). The transport of
water molecules through the pore is tightly coupled
to changes in their orientation.37,38 As water
molecules trickle through the pore, their polariz-
ation changes around the NPA motifs. As high-
lighted in Figures 1 and 2B, on average, water
ic GlpF channel is shown in ribbon representation (a) together
sidues, Arg206 in the selectivity filter, as well as Asn68 and
e highlighted along with the two half-membrane-spanning
Figure 2. Representative conformations of the three
polarization states of the chain of water molecules in the
GlpF pore. The three polar side-chains in the pore are also
shown. A, Fully polarized chain with molecular water
dipole moments oriented towards the periplasmic entry;
B, bipolar orientation; C, fully polarized chain with dipole
moments pointing towards the cytoplasm.
Proton Blockage in Aquaporins 495molecules 1–5 are polarized with hydrogen atoms
pointing toward the periplasmic vestibule, water
molecules 7–9 are polarized toward the cytoplasmic
vestibule, and water 6, which is located between
Asn203 and Asn68, donates both hydrogen atoms to
neighboring water molecules and orients perpen-
dicular to the channel axis. This bipolar organiz-
ation was ascribed to the opposed dipoles of the M3
and M7 helices33,38 as well as to hydrogen bond
Figure 1. Molecular model simulated here. The monomer
with a close-up of the pore region (b). Three conserved re
Asn203 from the signature Asn-Pro-Ala (NPA) motifs, ar
helices, M3 and M7. The nine single-file water molecules embedded in the pore are shown in their preferred bipolar
organization. Water molecules lying on the periplasmic and cytoplasmic side of the pore are shown at the top and at the
bottom, respectively. The molecular pictures in this and all subsequent Figures were generated with the program VMD.84

AquaporinsFour computational studies of the translocation
of an excess proton in AQPs were reported recently.
These studies used different models of the hydrated
proton and different model systems to compute the
free-energy profile of proton translocation. De
Groot et al. combined non-equilibrium molecular
dynamics simulations with stochastic proton jumps
using the QHOP method to determine the effective
free-energy profile opposing proton transport in
bovine aquaporin1 (bAqp1)40 using the maximum-
likelihood approach. Burykin & Warshel estimated
the free energy for transferring an excess proton to
the interior of the bAqp1 pore using the protein-
dipole–Langevin-dipole method.41 The potential of
mean force for both hop and turn steps of structural
diffusion in the narrow pore of the GlpF channel
was computed using a dissociable and polarizable
water model from equilibrium molecular dynamics
simulations with and without an excess proton,
respectively.42 Finally, Ilan et al. reported non-
equilibrium simulations of proton translocation in
the GlpF channel obtained with an empirical
valence bond model of protonated water.43
The four computational studies of proton trans-
location in AQP channels differ in their quantitative
estimate of the free-energy profile governing proton
movement and in their conclusions regarding the
physical origin and the mechanism of protondonation by the two amide groups of the NPA
motifs.38
Several proposals for the molecular origin of
proton blockage were put forward on the basis of
computational studies of the water chain. De Groot
& Grubmu¨ller proposed that both the resilient
disruption of the hydrogen bonded water chain
observed in the selectivity filter region and unfa-
vorable coulombic interactions between an incom-
ing cation and the conserved Arg residue could
prevent proton permeation through the pore of
Aqp1.37 Tajkhorshid et al. noted that the bipolar
organization of the water chain is incompatible with
the uptake of protons from either side and
speculated that the control of water orientation by
the channel would result in the exclusion of HC
from the single-file region.38 In that line of thought,
Jensen et al. used MD simulations to analyze
coulombic interactions of specific components of
the GlpF channel with water molecules and with a
probe ion in the water permeation pathway.39 They
concluded that electrostatic forces due primarily to
helices M3 and M7 and to the Asn side-chains of
NPA motifs give rise to the bipolar arrangement of
the water chain, thereby ensuring proton exclusion
from the pore. However, it should be noted that
none of the above studies considered the presence
of explicit protons in the pore.
Recent studies of proton translocation in
496exclusion. An overall free energy barrier to proton
translocation was computed as approximately 6,40
15,41 11,42 and 18 kcal/mol (1 calZ4.184 J).43 Thesedisparities reflect significant differences in the
models and methodologies used as well as in the
statistical convergence of the calculations. A com-
parison of free-energy profiles suggests that the
largest quantitative discrepancy is related to diffi-
culties in obtaining a reliable estimate of the
desolvation penalty for moving the excess proton
from bulk water to the channel interior (see
Discussion). Nevertheless, important qualitative
features consistently emerge. Most notably, all
four studies concur in the location of the barrier
top at the NPA motifs. In addition, the potential of
mean force (PMF) for proton movement shows a
shoulder at the selectivity filter in three of these
studies,40,42,43 indicating that the persistent inter-
ruption of the water chain in the absence of HC is
not sufficient to block the proton movement and
suggesting that charge repulsion with the con-
served Arg is not the principal determinant of
proton blockage. Significantly, there is a good
overall agreement between the results obtained for
the single-file region in the study of AQP1 by de
Groot et al.40 and that of GlpF by Chakrabarti et al.42
In and around the single-file region, the free-energy
profile obtained by Ilan et al.43 is similar, albeit with
a consistently larger amplitude.
Three of these studies analyze the role of
electrostatic forces in the free-energy barrier oppos-
ing proton translocation40–42 but differ in the
assessment of the relative importance of dielectric
and static effects in the mechanism of proton
exclusion. The total electrostatic profile for the
movement of a probe ion can be decomposed into
reaction field and static field contributions.44 The
former, which is due to dielectric boundaries,
reflects the penalty for desolvating the ionic charge
in the narrow pore, whereas the latter arises from
the distribution of electric charge in the protein. De
Groot et al.40 & Chakrabarti et al.42 compared the
PMF of proton transfer with Poisson–Boltzmann
continuum electrostatic calculations for the translo-
cation of a probe cation through the pore. The
similarity of the electrostatic profile with the free-
energy profile for the movement of HC supports the
dominance of the electrostatic field in the 20 A˚ long,
narrow part of the channel. Furthermore, Chakra-
barti et al.42 observed that the static field inside the
pore is consistent with the proton hop PMF both
qualitatively and quantitatively, and that these
profiles mirror the free energy profile opposing
the reorientation of the water chain, which consists
of a deep well strongly favoring the bipolar
conformation of the chain over either polarized
state. Together, these results indicate that dielectric
effects contribute significantly to the overall free
energy barrier for proton translocation and that the
distribution of polar and charged groups opposes
both hop and turn steps of structural diffusion in
the single-file region.42 This mechanism is consist-
ent with that put forward by de Groot et al.40 In
Proton Blockage in Aquaporinsaddition, the analysis of water orientation in the
absence and in the presence of an excess proton
indicates that, contrary to previous proposals,38,39,43

energy profile of ionic translocation.
cytoplasm, and water 6 donating one hydrogenPresent work
In the present study, we seek to determine the
structural origin of the barrier opposing proton
translocation through the pore of GlpF. Our
previous study suggested that it is the distribution
of charged and polar groups of the channel that
gives rise to the proton barrier inside the pore,42 but
did not provide detailed indications as to the
respective contributions of specific structural
elements. The balance of coulombic interactions
leading to the proton barrier is not trivial. For
instance, Arg206 at the periplasmic selectivity filter
does not appear, in itself, to block protons despite its
positive charge. Accordingly, a recent systematic
study of individual contributions to the electric field
for a probe ionic charge in the pore of GlpF suggests
that positive charge is offset by several negatively
charged groups in the periplasmic vestibule.39 In an
earlier study, the artificial turning off of coulombic
interactions suggested that the dipole moments of
M3 and M7 helices contribute to the bipolar
organization of the water single file.38 Here, we
extend the methodology used in a study42 to
consider the qualitative and quantitative effects of
structural elements on structural diffusion.
Equilibrium molecular dynamics simulations of
artificial variants in which coulombic interactions
between conserved structural elements of thethe bipolar orientation of the unprotonated water
chain does not in itself oppose proton translocation
through the narrow pore.42
The modulation of the PMF for proton movement
inside the pore by the charge distribution of the
channel40,42 contrasts with the main conclusion
reached by Burykin & Warshel, who argue that
the barrier to proton movement is essentially due to
dielectric effects and that it is desolvation penalties
(self-energy of a cation going through the channel)
that give rise to the free energy peak at the core of
the membrane, whose location fortuitously
coincides with the NPA motifs of aquaporins.41
However, this analysis leaves out some important
aspects of the mechanism. While it has long been
recognized that a self-energy penalty opposes the
passage of ions through the low-dielectric region of
lipid bilayers,45 it is a fact that many biological
channels have evolved to become selective to polar
and charged permeants by counterbalancing the
dielectric barrier.1 The dielectric barrier cannot by
itself lead to the overwhelming preference for the
bipolar organization of the unprotonated water
chain around the NPA motifs of aquaporins;42 such
bipolar arrangement is not observed in other
narrow membrane channels such as gramicidin,
where polarized states of the water chain prevail.21,46
In that context, it is essential to examine the relative
contributions of polar and charged groups to the free
Proton Blockage in Aquaporinschannel and the lumen contents of the pore are
turned off are used to compute the reversible
thermodynamic work, or PMF, for both hop andatom to each of these two polarized half-chains (see
Figure 1). In state C (mzw8 e A˚), all nine water
molecules point towards the cytoplasmic entry. The
PMF for the turn step of the Grotthuss mechanism is
depicted in Figure 3. In the native channel (wt),
there is an overwhelming preference for the bipolar
arrangement of water molecules, which is favored
by 7 kcal/mol and 12 kcal/mol over states A and C,
respectively.42
Artificial modifications of the charge distribution
of the channel strongly affect the organization of the
water chain. Turning off the charge of Arg206 (R
channel), leads to a dramatic stabilization of state A,
which becomes degenerate in energy with B. Turn-
ing off backbone atoms of M3 and M7 helices (H
channel) leads to the stabilization of states A and C,
both of which come closest in free energy to the
bipolar state. Turning off the partial charges of
the NPA Asn side-chains (N channel) leads to the
disappearance of the bipolar arrangement of the
single-file water chain as a stable or metastable
state; the PMF in the N channel is otherwise
identical with that of the wt channel. These featuresturn steps of structural diffusion. We chose to focus
on the role of the backbone of M3 and M7 helices
and of Asn68, Asn203, and Arg206, which are the
only polar side-chains in the narrow pore. These
artificial modifications probe the respective role of
dipole, hydrogen bonding, and charge interactions
with water and proton. The free-energy profiles
controlling the movement of HC are compared with
Poisson–Boltzmann (PB) calculations for the trans-
location of a probe cation throughout the water
permeation pathway. Results confirm that the
location of the barrier to protons at the NPA site is
due to the charge distribution of the channel and
indicate that the single most effective way to abolish
the barriers opposing structural diffusion is to omit
the macro-dipoles of helices M3 and M7, which
leads to a comprehensive model for the blockage of
structural diffusion by dipole–dipole and charge–
dipole interactions in the channel. This model is
exploited to propose single-point mutations that
could compromise the impermeability of aquapor-
ins to protons.
Results
Reorientation of the single-file water chain
The unprotonated single-file chain of water
molecules can adopt three distinct conformational
states (Figure 2). In state A, which corresponds to
mzwK8 e A˚, the dipole moments of all nine water
molecules are pointing towards the periplasmic
entry. State B, at mzwK1 e A˚, corresponds to a
bipolar organization, with five water molecules
pointing towards the periplasm, three towards the
497are retained qualitatively when NPA asparagine
residues and the two helix dipoles are turned off
together (HN channel), but the preference of

the absence of bipolar state B as a stable conformer
and a 6–7 kcal/mol preference for state A over state C.
Proton hopping
The PMF profiles for proton transfer along the
single-file water chain (hop step of the Grotthuss
mechanism) in native and modified versions of the
channel are shown in Figure 4. In the native channel
(wt), the hop profile is relatively flat around R206
and is opposed by a barrier of w4.5 kcal/mol
peaking at the NPA region.42 Turning off coulombic
interactions between the protonated water chain
and specific structural elements of the channel leads
to significant change in the free energy profiles. In
the variant and mutant forms of the channel
considered here, the barrier found in the wt channel
is replaced by a monotonic change, a flat profile, or
Proton Blockage in Aquaporins498polarized state A over C becomes stronger by
6 kcal/mol. Combining all three modifications (i.e.
turning off coulombic interactions between the
water chain and Arg206, Asn68, Asn203, and the
M3 and M7 helix backbones in the RHN channel)
results in an even stronger asymmetry of the PMF
profile than in the HN channel, further exacerbating
the preference for polarized state A, which is now
favored over state C by 21 kcal/mol.
Putting these results together, it appears that the
single most important effect of NPA Asn side-
chains on the organization of the water chain is to
allow the existence of state B, whereas the main
effect of helices M3 and M7 is to favor dramatically
the bipolar arrangement B over either polarized
states A and C (provided NPA is there to nucleate
the bipolar arrangement) and the primary effect of
Arg206 is to destroy the stability of state A relative
to the other conformations of the water chain.
Finally, both N68D and N203D mutations lead to
profiles resembling the turn PMF of channel N, with
a well. Turning off the charge of R206 in both R and
RHN channels has a profound effect on the PMF,
leading to a steep monotonic increase in free energy
between xwK2 A˚ and 7–9 A˚. Inversely, turning off
M3 and M7 helix dipoles (H) or NPA asparagine
residues (N) essentially abolishes gradients in the
PMF. These two profiles are highly similar to each
other, with a shallow well near xw0 A˚ between the
selectivity filter and the NPA region.
Turning off both Asn and helix dipole interactions
(HN channel) inverts the wt barrier, which is
replaced by a broad 3 kcal/mol well centered atFigure 3. Free energy profiles for the reorientation of
the single-file chain of water molecules (turn step of
structural diffusion) in the absence of an excess proton.
The potential of mean force is shown as a function of the
total dipole moment of the nine-water chain projected on
the channel axis z. The GlpF channel and its variants are
as described in Table 1. Labels A, B, and C refer to the
polarization states of the water chain defined in Figure 2.
The total dipole moment (mz) of the water chain was
scaled by a factor of 0.417, the fractional charge of H
atoms in the TIP3P potential,54 for consistency with the
results reported elsewhere.42Figure 4. Free energy profiles for the transfer of an
excess proton in the single-file chain of water molecules
(hop step of structural diffusion) for the wild-type GlpF
channel and its variants (see Table 1).

the NPA site (xZ6 A˚). A qualitatively similar result
is obtained in each of the two single-point mutants,
N68D and N203D, with a somewhat more pro-
nounced well of w4 kcal/mol located near the
carboxylate groups of the Asp side-chains, at xZ6 A˚
and 5 A˚, respectively. The similarity between HN,
N68D, and N203D profiles suggests that in each of
these two single-point mutants, charge–charge and
hydrogen bonding interactions between protonated
water and aspartate effectively cancel out the
combined effect of M3 and M7 helix dipole
moments and of the remaining Asn side-chain.
Comparison of hop PMF and static field profiles
The static field for the translocation of a positive
point charge in the water-conduction pathway of
native and modified forms of the GlpF pore was
computed with continuum electrostatic calculations
(see Materials and Methods). Results for the single-
file region are shown in Figure 5 together with the
PMF profiles presented in the above subsection. All
eight panels of Figure 5 show a good agreement
between the electrostatic potential and PMF pro-
files, indicating that the qualitative result obtained
earlier for the wt channel42 can be extended to the
channel variants. In particular, the gradients around
the NPA motifs (wt channel) and in the R and RHN
variants are well captured by the static field. In
energy wells in RHN, N68D, and N203D channels
are in very good agreement in both calculations,
indicating that the charge distribution dominates
the free energy profile of the excess proton in the
single-file region. By contrast, the profiles are in
relatively poorer agreement in cases where the PMF
profile is comparatively flat (i.e. in H, N, and HN
channels), suggesting that as the effective polarizing
force acting on the cation gets weaker, the role of
coordination and solvation of the excess proton
becomes increasingly important.
Static field across the permeation pathway
A direct comparison of static field profiles for the
translocation of a point charge obtained in native
and modified channels is shown in Figure 6. These
profiles are nearly identical outside the range
K30!z!30 A˚ but differ markedly in the region of
the channel pore. Results for the channel variants
can be grouped into four distinct categories. Turn-
ing off the partial charges of the two NPA Asn
amide groups has only a moderate effect on the
static field: both wt and N channels consist of a
2–3 kcal/mol barrier at the NPA site. In both
channel variants in which the charge of Arg206 is
turned off, R and RHN, the static field features a
Proton Blockage in Aquaporins 499Figure 5. Comparison of static field (broken lines) andaddition, both the location and curvature of free-hop PMF (continuous lines). The eight panels are
arranged as in Figures 2 and 3. The PMFs were shifted
vertically for a direct comparison.Figure 6. Static field across the permeation pathway.
Results obtained for different channels are as follows:
wild-type (black), and channel variants R (magenta), H
(green), N (blue), HN (cyan), RHN (red), N68D (walnut),
and N203D (gold). In this and subsequent Figures, the
magenta and blue bars highlight the location of the
selectivity filter (K5%z%K1 A˚) and of the fingerprint
NPA motifs (4.75%z%7 A˚), respectively. The barrier due
to the charge distribution for the wt channel peaks to
3 kcal/mol at the NPA site, whereas for channels R and
RHN, the barrier is replaced by a well at the selectivity
filter. The full profile for the wt channel is shown for the
region K80%z%80 A˚ as an inset together with the effect
of a membrane potential on static field in the wt GlpF
channel. The static field for the wild-type channel is
shown without any trans-membrane potential (black) and
with a transmembrane potential of K65 mV applied
along the channel axis (orange). The static field for the
H channel (no voltage) is depicted in green. The
membrane region, which corresponds to the range
K15 A˚!z!25 A˚, is highlighted with broken lines.

HN and RHN variants, with broad wells and
helices leaves the channel largely non-polar (see
Figure 1). Although it remains substantial in the R
form, the NPA barrier drops to zero in the RHN
variant; most significantly, the shoulder at the
selectivity filter is replaced by a well in both R
and RHN variants. Finally, the ESP features largely
vanish in the two ND mutants, where the amplitude
of the total ESP, which oscillates between 2 kcal/
mol and K2 kcal/mol, is smallest.
The latter result stems from the near-cancellation
of static and reaction field contributions to the total
ESP in the two ND mutants. A typical conformation
of the protonated water chain in the N68D mutant is
shown in Figure 8, and the decomposition of the
total electrostatic free energy for the translocation of
a probe cation in that mutant is depicted in Figure 9.
The static field yields a well (w6 kcal/mol) at the
NPA site, whereas the reaction field, which is
related to desolvation penalties, features a barrier
(w8 kcal/mol) peaked at the same region. The total
electrostatic free energy is the sum of static and
reaction fields within approximately 0.5–1.0 kcal/
mol. Although the ruggedness of the reaction fieldenergy minima of 6–7 kcal/mol located between
the selectivity filter and the NPA region. The
respective locations of the minima along the
channel axis reflect the fact that N203 is closer to
the periplasmic mouth of the narrow pore than N68.
A pairwise comparison of the results obtained for
wt and N, H and HN, and R and RHN channels
suggests that the static field contributions of
Arg206, the two NPA Asn side-chains, and the
backbone dipole moment of M3 and M7 helices are
roughly additive, with the extremum in the single-
file region dropping by approximately 1 kcal/mol,
3 kcal/mol, and 11 kcal/mol when the partial
charges of the Asn side-chains, of the backbone of
M3 and M7 helices, and of Arg206 are turned off,
respectively.
Total electrostatic energy of a permeating probe
cation
The total electrostatic free energy for the move-
ment of a probe cation through the water per-
meation pathway in various forms of the channel is
shown in Figure 7. The total electrostatic energy can
be decomposed as the sum of static field and
reaction field contributions, which arise from
charge–charge interactions and from dielectric
boundaries, respectively. Since the only modifi-
cations in the artificial variants of the channel
involve perturbations of coulombic interactions
between the channel contents and the charge
distribution of the channel, the only part that differs
significantly in the five artificial constructs con-
sidered here (R, H, N, HN, and RHN) is the static
field (as shown in Figure 6). The reaction field is
identical in the wt and in these five variants,
because by design, none of these modifications
induces a significant change of the conformation ofdeep energy well of 8 kcal/mol and 10 kcal/mol,
respectively, centered at the selectivity filter region.
By contrast, both channel variants in which M3 and
M7 helix dipole moments are omitted, H and HN,
show only minor deviations from the baseline,
which consists of a monotonic increase from about
K2 kcal/mol to 0 kcal/mol over the range K25 A˚
to 35 A˚ in the absence of a membrane voltage. This
confirms that turning off the dipole moments of M3
and M7 helices is the single most effective way to
abolish the electrostatic gradients opposing the
movement of a cation in the narrow pore of the
channel. Combining H and N modifications (HN
variant) cancels the barrier obtained in the wt
channel. In the presence of a transmembrane
voltage of K65 mV typical of physiological values
(K40 mV toK100 mV),1 the asymmetry of the static
field in the wt channel is inverted and reduced in
magnitude (Figure 6).
Finally, both N68D and N203D mutations result
in static field profiles that are intermediate between
500the channel. In addition, the reaction field is
retained in the isosteric single-point mutants
considered here (N68D and N203D). As a conse-quence, the ordering in the ESP profiles (see
Materials and Methods) shown in Figure 7 is the
same as that of the static field shown in Figure 6.
In the native channel (wt), the ESP profile consists
of an 11 kcal/mol barrier centered at the NPA site,
with a shoulder at the selectivity filter. The profile
drops off on both sides of the narrow pore, reaching
zero at zZK12 A˚ and 24 A˚ on the periplasmic and
cytoplasmic sides, respectively. Although the mag-
nitude of both the NPA barrier and the shoulder are
reduced, these features are retained in the N, H, and
HN variants. By contrast, neutralizing Arg206
results in strongly asymmetric ESP profiles. Turning
off R206, N68, N203, and backbone of the two
Figure 7. Total electrostatic free energy for the
translocation of a probe cation along the water per-
meation pathway in native and variant forms of the
channel. Colors are as described for Figure 6.
Proton Blockage in Aquaporinsin the pore region remains in the total ESP profile,
these fluctuations are essentially confined to the
narrow pore, where the reaction field is overly

tude in that region. Replacing Asn68 by Asp
Figure 8. Representative conformation of the proto-
nated water chain in the N68D channel, and the channel
groups forming hydrogen bonds with the single-file
water chain are shown. In this conformation, the excess
proton resides on water 6, which is located near the
carboxylate group of D68.
Proton Blockage in Aquaporinsreduces the overall ESP barrier from w10 kcal/
mol to a much smaller value of w4 kcal/mol,
suggesting that this mutant (and N203D, for
which similar results are obtained) may conduct
protons.
Discussion
In this section, the mechanism of proton blockage
in aquaporins is discussed in light of the above
results and of recent computational studies. This
analysis sheds a light into the balance ofsensitive to molecular details (see Discussion).
However, the drop in the static field as the probe
cation is taken from bulk water into the narrow pore
(approximately 4 kcal/mol and 5 kcal/mol on
periplasmic and cytoplasmic ends) roughly cancels
the desolvation penalty, which is of similar magni-Figure 9. Decomposition of the total electrostatic
energy (continuous line) for the translocation of a probe
cation along the permeation pathway of the N68D
channel into (broken line) static field and (dotted line)
reaction field contributions. The vertical broken lines
highlight the boundaries of the single-file region at
zZK6 A˚ and zZ14 A˚, respectively.fundamental physical forces at play in the biological
control of proton translocation.
Static field
The results of the present study support the
dominance of electrostatic effects in the mechanism
of proton blockage, a property uncovered in three
recent independent studies.40–42 The respective
contributions of static and reaction fields to the
electrostatic field are discussed in this and in the
next subsection.
The effect of modifications of coulombic inter-
actions between the channel and its pore contents
confirms that both hop and turn steps of structural
diffusion in the narrow pore of aquaporins are
opposed by the charge distribution of the channel.42
The essential role played by the charge distribution
in the free-energy profile for proton transfer within
the narrow pore of GlpF is supported by the static
field for the translocation of a probe cation. The
PMF for proton movement is nearly identical
(within 1 kcal/mol) with the static field in the
single-file region, both in the native channel and in
the five artificial variants as well as in the two
single-point mutants of the channel considered
above. Discrepancies between the two profiles are
due, in part, to factors neglected in continuum
electrostatic calculations, such as thermal averaging
and the delocalized nature of the charge of the
hydrated proton in a hydrogen bonded network.42
Nevertheless, the above results indicate unambigu-
ously that in the pore of aquaporins, significant
electrostatic gradients take precedence over the
details of hydrogen bonding coordination and
transfer properties of the hydrated proton. This
property is supported by the agreement between
the proton hop PMF profile throughout the single-
file region in two previous studies using different
model systems and methodologies,40,42 and is
consistent with the conclusions of three recent
computational studies.40,41,43 The thermodynamic
profile for proton movement is dictated by the
strong electrostatic pull away from the NPA motifs.
While this electrostatic effect is not specific to
protons, this does not necessarily mean that the
detail of short-range interactions involved in ionic
solvation would be irrelevant for other ions (anions
or cations), as discussed elsewhere.39,40,42
Desolvation penalties
Together with the charge distribution of the
channel, desolvation penalties arising from dielec-
tric boundaries contribute to the blockage of
protons. Importantly, the present study suggests
that the respective influence of these two effects is
different inside and outside the pore. As discussed
above, the static field is by far the largest contri-
bution to the free energy opposing proton per-
501meation inside the narrow single-file region. This
necessarily means that the reaction field should be
essentially constant inside the pore, in contradiction

crepancies between the results obtained in various
to the results obtained from Poisson–Boltzmann
calculations, whereby the reaction field varies by up
to 3 kcal/mol between the selectivity filter and the
NPA motifs (Figure 9). On the basis of the present
study we can discount altogether the two peaks of
the reaction field appearing at the selectivity filter
and the NPA motifs in the narrow pore region. The
contradiction indicates that dielectric boundaries
are an inadequate approximation to describe the
self-energy of the ion within the single-file region of
the pore, where the reaction field is overly sensitive
to small variations in the channel width. This
conclusion could not have been reached on the
sole basis of the study of the wt channel, where
static and reaction fields happen to be qualitatively
similar in the pore region.42
By contrast, outside the narrow pore, variations
in the static field, which reaches K2 kcal/mol and
0 kcal/mol near periplasmic and cytoplasmic
entrances, respectively, are relatively small (Figure
6). This suggests that in the vestibules of the
channel, it is desolvation penalties that dominate
the overall electrostatic barrier to proton transloca-
tion. With the present methodology, we obtained a
total electrostatic barrier peak of 11 kcal/mol at the
NPA site and a shoulder 5 kcal/mol below.42 If
fluctuations of the reaction field inside the pore are
discounted, the peak of the electrostatic free energy
profile drops to approximately 8 kcal/mol, flanked
by a shoulder at the selectivity filter 3.5 kcal/mol
below. Outside the single-file region, the total
electrostatic free energy drops by 4 kcal/mol.
Although crude, this estimate of the free-energy
barrier outside the pore obtained using macroscopic
continuum electrostatic calculations is consistent
with the results obtained from three previous
studies, where the free-energy profile for the
entrance and exit of HC from bulk to bulk was
computed using explicit atomic models.40–42
Our results are in good agreement with those
reported by de Groot et al.,40 despite the fact that
their free-energy profile reaches minima outside the
pore, then rises again by w2 kcal/mol and 4 kcal/
mol towards bulk water on periplasmic and
cytoplasmic sides, respectively. Between these two
minima, their PMF features a barrier of 8–10 kcal/
mol peaking at the NPA, a shoulder at the
selectivity filter w3–4 kcal/mol below the barrier
top, and drops ofw4 kcal/mol between the ends of
the narrow pore and the free-energy minima in the
vestibule regions. Thus, our single-file PMF profile
and our estimate of 4 kcal/mol for proton entrance
into the narrow pore form either side of the channel
are consistent with the magnitude obtained by de
Groot et al.39 The same magnitude, with somewhat
larger estimates of 6 kcal/mol and 7 kcal/mol for
proton entry from periplasmic and cytoplasmic
sides, was obtained by Ilan et al.,43 while in the work
of Burykin & Warshel,41 the desolvation penalty
appears to be in the range of 10 kcal/mol. A direct
502comparison of the results obtained in various
studies is hampered by the difficulty of gauging
the boundary of the single-file region, where thestudies reflect the difficulty to obtain reliable
estimates of the transfer free energy of an excess
proton between bulk water and the interior of the
narrow pore. Uncertainties in these calculations
include differences in the molecular models and in
the empirical force-fields used, as well as statistical
sampling errors. None of the models included the
lipid membrane or other monomers of the channel
explicitly. Finite-size models were used by all but
one study.40 Some models did not include the effect
of proton delocalization.40,41 In non-equilibrium
simulation approaches used in two studies,40,43
the water content of the lumen can change with
time, so that the simulations may probe hydration
states different from those in here, where the intrusion
and extrusion of water in and out of the pore were
precluded by construction.42 As the channel vesti-
bule widens, the PMF becomes more and more
difficult to calculate accurately, due to increasing
dimensionality and multiple pathways, whereas
the continuum electrostatic approximation is
expected to become better and better. In that
context, it is significant that the magnitude of the
desolvation barrier estimated by a continuum
model in our earlier42 and present studies is
consistent with free energy changes obtained with
atomistic force-fields.40,41,43 Together, these results
confirm that electrostatic effects are dominant in the
exclusion of protons, with the dielectric term
prevailing outside the narrow pore and the static
field compounding this desolvation barrier by
providing an additional barrier at the NPA site.
Structural determinants of proton blockage
We now turn to the analysis of modulations of the
static field by specific interactions and discuss how
single-point mutations can mitigate and possibly
even compensate for the dielectric barrier opposing
proton translocation. The above results enable us to
identify the respective contributions of specific
structural elements and of physical forces to the
mechanism opposing structural diffusion in the
narrow pore of the GlpF channel. The analysis of
artificial variants indicates that all three structural
elements considered are needed to enforce both the
strong preference for the bipolar organization of the
unprotonated chain and the presence of a free-
energy barrier opposing the passage of an excess
proton at the NPA site. Moreover, the detailed
comparison of hop and turn free-energy profiles
sheds a light into the relative effects of charge
(Arg206), hydrogen bonding (Asn68 and Asn203),
and dipolar (M3 and M7 backbone) interactions on
the mechanism of proton blockage.
Asn side-chains and helix dipolesfree-energy gradients are steep. Furthermore, dis-
Proton Blockage in AquaporinsThe two Asn side-chains of the NPA motifs and
helices M3 and M7 play complementary roles in the
stabilization of the bipolar conformation.

for proton blockage in aquaporin channels, which isOther polar and charged groups
The static field features a 2.5 kcal/mol well at
zwK20 A˚ in the periplasmic vestibule of the GlpF
channel. This feature persists qualitatively in all the
variant forms, including in H and HN channels,
even as the barrier to proton translocation in the
single-file region is being abolished, and culminates
in R and RHN variants, where the monotonic
increase of the static field from the periplasmic
mouth to the cytoplasmic side of the pore is most
pronounced. Similarly, the PMF profiles for the
reorganization of the unprotonated water chain
(Figure 3) show a resilient preference for state A,
which is favored byw7 kcal/mol over state C in the
wt channel, and again (by 4–21 kcal/mol) in the
variant forms of the channel. Together, these results
indicate consistently that the distribution of
charged and polar groups outside the narrow pore
of the channel is strongly asymmetric. This asym-
metry is due to the presence of negatively chargedArg206
Turning the charge of Arg206 off generates a
strong asymmetry in the PMF profiles for both hop
and turn steps of structural diffusion. The static
field profile in the R channel suggests that a cation
might reside in the periplasmic vestibule. However,
the amplitude of both static and total fields is very
large, suggesting that single-point mutations
removing the positive charge of residue 206
would not lead to cationic leakage. This result
confirms earlier conclusions that Arg206 is not the
dominant structural feature opposing structural
diffusion in the pore as it fails to polarize the
unprotonated water chain39,42 and it gives rise to
only a secondary peak in the PMF opposing proton
permeation.40–43Neutralizing Asn68 and Asn203 amide groups
prevents the formation of hydrogen bonds with
water molecules at the NPA site, which is sufficient
to eliminate the bipolar conformation of the water
chain (Figure 3). Thus, hydrogen bonding is
necessary to nucleate the bipolar conformation,
although it does not contribute as much to the free-
energy barrier for proton hopping past the NPA
motifs as the two helix dipoles (Figures 4–6). The
single most effective way to remove the free-energy
barriers opposing both hop and turn steps of
structural diffusion is to turn off the dipole
moments of helices M3 and M7 (H channel). This
modification results in the near-degeneracy of all
three polarization states of the unprotonated water
chain, and in a diffusive free-energy profile for the
excess proton in the single-file region (Figures 5–7).
Together, these two modifications (HN channel)
result in the disappearance of the static field barrier
(Figure 6).
Proton Blockage in Aquaporinsgroups near the periplasmic mouth of the pore. In a
previous study of electrostatic interactions in GlpF,
these groups were identified as E43, E152, anddepicted in Figure 10. According to this mechanism,
the preferred arrangement of water molecules in
GlpF is determined essentially by helix dipole–
water dipole interactions and the barrier to proton
movement inside the pore is dominated by adverse
interactions between the ionic charge and helix
dipole moments. Throughout the present study, the
consistent effects of structural elements on hop and
turn profiles underlines the complementarity of the
two steps of structural diffusion. As a confirmation
and a refinement of our earlier conclusions,42
Figure 10 shows how the very same structural
features leading to the stabilization of the unproto-D207.39 Negatively charged residues on the peri-
plasmic side of the pore include D114, D130, E116,
and E127. The intrinsic asymmetry of the static field
of the channel may be implicated in the effective-
ness of cationic blockage at physiological con-
ditions. Accordingly, the electrostatic profile
obtained in the presence of a membrane voltage of
K65 mV (Figure 6) shows that the static field of the
channel partly compensates for this field,
suggesting that aquaporins may have evolved to
function in the presence of an intracellular electro-
chemical gradient.
NPA mutants
Modifying one of the two NPA asparagine
residues to aspartate offers further insight into the
balance of physical forces at play in structural
diffusion. In the absence of a proton, the free-energy
profile governing water orientation in N68D and
N203D channels is similar to that of the channel
variant in which the partial charges of the two Asn
side-chains are turned off (N channel), underlining
the importance of replacing a hydrogen bond donor
by an acceptor to the structure and fluctuations of
the water chain. Replacement of CONH2 by COO
K
changes the free-energy barrier opposing proton
translocation at the NPA motifs to a well. Attractive
charge–charge interactions with HC more than
compensate for the effect of the other charged and
polar groups of the channel. Since the extra
electrostatic stabilization afforded by the presence
of a negative charge at the NPA site roughly cancels
out the desolvation penalty as well, our results
suggest that both N68D and N203D could mediate
proton translocation.
Mechanism of proton blockage
The present study shows that the single best way
of eliminating the barrier to proton transfer in the
single-file region is to turn off the dipole moments
of helices M3 and M7. This artifact also leads to the
destabilization of the bipolar conformation of the
water chain relative to polarized conformations.
These results lead to a comprehensive mechanism
503nated water chain with a turn defect at the NPA site
also lead to the destabilization of the ionic defect at
the NPA site. Although Figure 10A and B depict the

protons (ca K100 kcal/mol) is an order of magni-Figure 10. Effect of the charge distribution of the
channel on structural diffusion in aquaporin channels.
A, The bipolar organization of the unprotonated water
chain is stabilized by dipole–dipole interactions with
helices M3 and M7. B, Protonation of water near the NPA
site conserves the bipolar organization of the water chain
but is opposed strongly by repulsive forces between the
ion and M3, M7 helix dipoles. This results in a barrier
504two extrema in the thermodynamic cycle of
structural diffusion (respectively, free energy well
for the turn and barrier top for the hop), the only
difference between them is the addition of an ionic
charge on the water molecule located at the NPA
site (Figure 10B). In this conformation, both ion–
water and water–channel interactions are optimal,
so that it is solely by virtue of adverse ion–channel
interactions that the barrier peaking at the NPA
arises.
A consequence of that observation is that,
contrary to an earlier proposal,38,39,43 bipolar con-
trol of the water chain does not in itself oppose the
transfer of an excess proton in the narrow pore
region.42 If that was the case, then the location of an
excess proton at the NPA site would be a stable
conformation by virtue of the bipolar arrangement
of water molecules around the NPA site (see
Figure 10B). As a consequence, bipolar control of
the water chain would oppose proton movement
away from the NPA motifs, since the displacement
of HC in either direction displaces the bipolar
arrangement away from the NPA site.42 Thus, the
fact that Figure 10B corresponds instead to a barrier
top indicates that, despite the strong preference for
a bipolar arrangement of the water chain, bipolar
control is not the dominant mechanism of proton
blockage.
peak for the PMF of proton transport (see Figure 4)
despite favorable water-channel and ion–channel inter-
actions. C, In N68D and N203D mutants, the negative
charge of the carboxylate group compensates for the
repulsive interactions between a cation and helix dipole
moments, resulting in the stabilization of the excess
proton.tude larger than the reorientation penalty imposed
on the water chain by the channel (7–14 kcal/mol).
Thus, proton exclusion results from desolvation
penalties and from the static field of the channel,
which also dictates water orientation in the absence
of HC. In this process, water molecules play an
altogether passive role, their orientation reflecting
the static field both in the absence and in the
presence of the excess ion.
The introduction of a negative charge at the NPA
site through ND point mutations results in the
relative stabilization of the proton at the NPA site,
which counterbalances the adverse interactions of
the cation with the helix dipoles (Figure 10C). This
effect is supported by the PMF and static field
profiles of proton movement in both N68D and
N203D mutants, where the barrier opposing proton
translocation is replaced by a well. In the N68D
mutant in particular, the shape of that well mirrors
qualitatively the wt barrier, with a primary well
centered at NPA and a plateau at the selectivity
filter (Figure 5). In further support of that cancella-
tion mechanism, we note that the N68D mutation
has the same effect on the turn step of structural
diffusion as turning off the partial charges of N68
and N203 (N channel): in the absence of HC, the
bipolar conformation of the water chain disappears
completely (see Figure 3).
General implications
The above mechanism is consistent with emerg-
ing evidence for the role of short-range and long-
range charge–dipole interactions in biological ion
transport. In the mechanism leading to the per-
meation of alkali metal ions through gramicidin,11
and of potassium ions through the selectivity filter
of the potassium channel KcsA,3 backbone carbonyl
groups solvate the permeating ion directly via
charge–dipole interactions. The sensitivity of hop
and turn steps of structural diffusion to dipole
moments in the channel was uncovered in a recent
comparative theoretical study of proton transloca-
tion in native and modified gramicidin dimers,Furthermore, while the magnitude of the free-
energy barriers opposing water reorientation
(Figure 3) might appear to support a direct role
for the bipolar organization of the water chain in
excluding protons from the pore interior, such an
observation is incomplete because it leaves out the
effect of ion–water interactions on the orientation of
the chain. The spontaneous polarization of the
single-file water chain induced by the addition of
HC at either end of the pore indicates that ion–
water interactions surpass channel–water inter-
actions,42 underlining the plasticity of the water
chain in response to the static field created by
protein and excess ion. This result is not surprising
in light of the fact that the hydration free energy of
Proton Blockage in Aquaporinswhere the origin for the attenuation of proton
conductance in dioxolane-linked gramicidin chan-
nels was attributed to local structural distortions of

the peptidic backbone that project carbonyl groups
into the lumen.27
The present results and analysis indicate that
essential aspects of the free-energy profile govern-
ing proton movement through aquaporins can be
captured with macroscopic continuum electrostatic
calculations. The largest uncertainty in these calcu-
lations was identified as the difficulty of describing
the narrow pore of the channel with dielectric
boundaries. The concordance of results obtained
from atomistic and a continuum electrostatic
calculation suggests that detailed features of struc-
tural diffusion are secondary in the mechanism of
proton exclusion from aquaporins. This contrasts
with the molecular mechanism of rapid proton
translocation in gramicidin, where hydrogen-bond-
ing interactions favor the solvation and the mobility
of protonic and bonding defects.21,26 This is because
the features opposing proton transport need be only
relatively crude or coarse compared to conduction
processes, as noted elsewhere,42 and because they
are largely non-specific, in the sense that they are
not due to the arrangement and atomic fluctuations
of a particular subset of chemical groups lining the
narrow channel but rather to the macro-dipoles of
two a-helices, which are not in direct contact with
the lumen contents.
Whether the conclusions reached above regard-
ing the structural origin of proton blockage would
extend to other cations will depend on the
energetics of these ions in the narrow pore of the
channel.42 Nevertheless, it should be noted that an
advantage of controlling cation movement with
relatively long-range charge–dipole interactions, as
is achieved in the GlpF channel by the head-to-head
arrangement of helices M3 and M7, is that this effect
would apply equally to all cations small enough to
fit in the single-file region. This result mirrors the
electrostatic control of charge selectivity in the
potassium channel KcsA,13 where the combined
dipole moments of four helices help stabilize a
single positive charge in the heart of the transmem-
brane region. Inversely, two head-to-head helices
are found in the ClC chloride channel,3 although
their role in the electrostatic stabilization of chloride
ions appears to be marginal compared to the case of
KcsA.47 At any rate, the impermeability of aqua-
porins to anions arises necessarily from effects other
than non-specific charge–dipole interactions with
the two half-helices.39,40,42 The asymmetry of the
static field for cation translocation, which features a
well near the periplasmic mouth (Figure 6),
suggests that the presence of negatively charged
groups near the periplasmic mouth of the channel
would contribute to opposing the passage of
anions, as pointed out by de Groot et al.40
Our analysis indicates that in the present system,
variations in the free-energy profile due to details of
the proton relay mechanism arising from the
reactivity of protons in hydrogen bonded networks
Proton Blockage in Aquaporinsare secondary to significant electrostatic gradients.
Accordingly, the main factors resulting in the
exclusion of protons from aquaporins can beunderstood without consideration of structural
fluctuations of the hydrogen bonded network,41
although the analysis of structural diffusion helped
reach a detailed understanding of the blockage
mechanism. However, this does not imply that the
structural fluctuations are necessarily irrelevant to
understanding the mechanism of proton movement
in proteins.48 Again, this is because, from the
perspective of molecular design, blockage is much
easier to achieve than selective conduction. In
aquaporins, static and dielectric effects act in
concert to block protons, whereas in channels that
conduct protons, the two effects must necessarily
balance each other. In addition, in proton-conduct-
ing channels, dynamic effects must be considered at
the atomic level to explain high rates of permeation.
In such systems, the detailed properties of hydro-
gen bonded networks become essential because
they modulate the rate of permeation. Dynamic
modulation of free-energy profiles by confor-
mational fluctuations of dioxolane-linked gramici-
din channels revealed a coupling between channel
dynamics and proton movement on a timescale
commensurate with proton permeation and suggest
that conformational transitions modulate proton
conductance in these channels.27
Dynamic effects must a fortiori be considered in
detail in proton pumping, an inherently non-
equilibrium process where kinetic control is
required to prevent the leakage or back-flow of
HC and ensure the directionality of proton translo-
cation. In principle, the mechanism of a proton
pump can be formulated as a series of reactions
alternating blockage and conduction through local
or extended tracts of the proton pathway, and an
important challenge in understanding the molecu-
lar mechanism of proton pumps is to determine the
structural and physical basis for the fine modu-
lation of these two effects (proton switch). Thanks to
the elucidation of high-resolution atomic structures
obtained at various stages of its photocycle,
bacteriorhodopsin is currently the best-character-
ized proton pump; nevertheless, the nature of the
proton switch in that enzyme is still a matter of
debate (see Lanyi49 and references therein). Two
limiting mechanisms have been invoked to explain
the directionality of proton movement in pumps:
affinity switch and accessibility switch.50 In the
former case, proton gating arises from rapid
changes in the relative affinity of relay groups for
protons triggered by photochemical (or redox)
reactions, whereas accessibility switch models
emphasize temporary interruptions of the proton
pathway resulting from conformational changes in
the hydrogen bonded network. It should be noted
that these two models are not necessarily exclusive.
A continuum electrostatic study considering
changes in proton affinity of key titratable groups
in various conformational microstates of the photo-
cycle of bacteriorhodopsin concluded that these
505changes could drive unidirectional proton translo-
cation without the need for accessibility switches.50
Accordingly, electrostatic gradients in aquaporins

detailed analyses of structural diffusion throughout
enzymatic pumping cycles are needed a priori to
forces governing structural diffusion in biological
systems and have general implications for the
with TIP3P and with the PM6 model.21,22,55–57 The TIP3Pcontrol of proton translocation in membrane chan-
nels and energy-transducing enzymes.
Materials and Methods
The calculations reported here follow the methodology
of an earlier study of the native GlpF channel.42 Here,
these calculations are repeated in five variant and two
mutant forms of the channel.
Molecular modeldetermine the interplay of electrostatic forces and
structural fluctuations in such systems.
Conclusions
The role of ion channels is to catalyze the
transport of charged molecular species through
the low dielectric region of the membrane. Trans-
port is achieved by compensating for the unfavor-
able desolvation penalty using favorable
interactions between the permeating species and
polar or charged groups of the channel. The present
study confirms that electrostatic forces dominate
the mechanism of proton blockage in aquaporin
channels and provide detailed structural insight
into the molecular origin of blockage. The desolva-
tion barrier opposing the intrusion of HC into the
narrow pore is compounded by a static barrier
peaking at the NPA motifs, which opposes both the
movement of protons and the reorientation of water
molecules inside the pore. The static field is due
primarily to the macrodipoles of two head-to-head
a-helices converging at the NPA motifs. These
structural elements provide a robust electrostatic
basis opposing both hop and turn steps of structural
diffusion through charge–dipole and dipole–dipole
interactions, respectively. Such a mechanism
ensures the exclusion of protons and possibly
other cations. Replacement of either of the two
NPA Asn side-chains by negatively charged Asp
compensates for adverse dipolar interactions and
approximately cancels out the dielectric penalty,
suggesting that the proton impermeability of these
two single-point mutants may be compromised.
These findings shed light on the balance of physicalsuffice to block protons despite the presence of a
proton pathway.40–42 Inversely, proton blockage
may arise regardless of thermodynamic and elec-
trostatic gradients given persistent interruptions of
the hydrogen bonded pathway, as suggested by a
recent study of proton blockage by methanol in the
narrow pore of the gramicidin channel.51 An
implication for the study of proton pumps is that
506The GlpF channel is a homo-tetramer whose three-
dimensional structure was solved at 2.2 A˚ resolution by
X-ray crystallography.35 In this work, one monomericforce-field was used in the study of the turn step of the
Grotthuss mechanism (in the absence of an excess
proton), whereas the PM6 force-field was also employed
in simulations with an excess proton to estimate the
proton-hop free-energy barrier. A comparison of results
obtained with TIP3P and PM6 models in the GlpF channel
indicates that these two models consistently describe the
structure and fluctuations of the single-file water chain,42
despite essential differences in the functional form of
these two empirical force-fields. More specifically, the two
models were shown to be in good agreement regarding
the equilibrium distribution of water in the pore and the
PMF for the reorganization of the hydrogen bonded chain
(turn step). PM6 is a polarizable and dissociable model of
water that consists of O2K and HC moieties. This
empirical model, which has been used in many studies
of ionized water,58–63 has been shown to capture the
essential features of the mechanism of HC transport in
water wires.23,26 Hop and turn PMF and diffusion
constants obtained from molecular dynamics (MD)
studies of proton transport in the gramicidin A channel26
were used in a framework model to compute conductance
data.64,65 The results of these studies were shown to be
consistent with experimental data.14,65,66
A promising approach to construct an empirical
potential function of water molecules that allows proton
dissociation is based on Warshel’s empirical valence bond
(EVB) theory.67 This approach gave rise to a new family of
dissociable water models initiated by Vuilleumier &
Borgis,68 and later refined by Voth and co-workers.69
The refined model has been shown to yield a good
agreement with the properties of the hydrated proton in
bulk water. However, the only calculations of proton
conductance through biological channels based on atomic
models published to date are based on studies of
gramicidin using the PM6 model.26,65,66 Thus, PM6 is
currently the only force-field that has been tested and
shown to yield adequate performance with respect to
experimental proton conductance data in biological
channels. Until superior accuracy is demonstrated by
other empirical models through rigorous comparison to
experimental data in narrow pores, the PM6 model
remains a valuable device to help understand the
molecular basis of proton movement in biomolecular
systems.
Molecular dynamics simulations
The inner core of the molecular system, an orthorhom-
bic region of 30 A˚!14 A˚!14 A˚ centered on the NPA
motifs and aligned with the narrow water-filled pore, was
allowed to move during the MD simulations. This region
is comprised ofw1100 atoms;42 the rest of the system was
kept fixed. Additional restraints were applied to prevent
intrusion and extrusion of water molecules in the single-GlpF channel, with slabs of water molecules above and
below it, was used for the simulation (see Figure 1). The
monomeric system consists of 3839 protein atoms, nine
water molecules in the pore and 1383 water molecules in
the bulk, for a total of 8315 atoms, as described.42
The CHARMM force-field, version 2252,53 was
employed to model the protein and the TIP3P force-
field54 was used for bulk water. The nine water molecules
occupying the channel pore38 were modeled successively
Proton Blockage in Aquaporinsfile region, as described.42 The distribution and orien-
tation of water molecules in the pore in this finite-size
system were shown to be in good agreement with those

ively. Constants dsw and rsw were chosen as 0.05 A˚ andobtained from a simulation of the tetrameric channel
assembly in a hydrated lipid bilayer.38,42
We generated MD trajectories using the CHARMM
program, version 26.52 The Langevin equations of motion
were propagated at 300 K with an integration step of 1 fs
and a friction coefficient of 5 psK1 applied to all heavy
atoms. Non-bonded interactions were calculated without
any cut-off. Five artificial variants of the channel were
constructed in order to gauge the effect of specific
coulombic interactions between structural elements of
the channel and the single-file water chain. These variants
are listed in Table 1. In each of these variants, pairwise
coulombic interactions between specified groups of
channel atoms and the lumen contents were turned off,
while channel–channel interactions were conserved. This
procedure has the advantage of guaranteeing the struc-
tural integrity of the pore. Two single-point mutant forms
of the channel were generated by replacing Asn68 and
Asn203 by Asp. The free-energy calculations described
below for hop and turn steps were repeated for each of the
seven variants.
Water reorientation
The turn step of structural diffusion was studied in the
absence of an excess proton from molecular dynamics
Table 1. Systems studied
Selective control of electrostatic interactions
Channel Arg206
M3/M7
helices
backbone
a
NPA (Asn68/
203) NH2
wt On On On
R Off On On
H On Off On
N On On Off
HN On Off Off
RHN Off Off Off
Point mutation
b
N68D Modify Asn68 to Asp
N203D Modify Asn203 to Asp
a Partial charges of C]O, N–H, and CaHa atoms of residues
70–79 (M3) and 205–217 (M7).
b Modifying either of the two NPA asparagine residues.
Proton Blockage in Aquaporinstrajectories in which the configuration of the molecular
system was recorded every 5 fs. The reversible work
theorem was used to compute the PMF for the
reorientation of the single-file chain of water molecules.
The reaction coordinate for this free energy calculation
was the projection of the total dipole moment of
the lumen contents (water chain) on the z-axis
(channel axis) for each configuration as a function of
time t, i.e. mzðtÞZqO
P
Oi zOiðtÞCqH
P
Hj zHjðtÞ. Here, mz(t)
is expressed in units of e A˚, and qO, qH are the formal
charges of single-file water O and H atoms, respectively.
For the TIP3P water model,54 qO and qH are K0.834 e and
0.417 e, respectively, where e is the elementary charge.
This reaction coordinate has been used in previous
studies of the Grotthuss turn step.22,26,27,42 The PMF
profile was calculated from the equilibrium probability
distribution of mz. A harmonic biasing potential energy
function was imposed on mz to force the reorientation of
the unprotonated water chain, as described.42 In this
umbrella sampling scheme,70,71 an equilibration run of
20 ps was followed by a 60 ps production run for data
collection for each of 41 consecutive windows separated
by 0.25 e A˚ and ranging from K5.0 e A˚ to 5.0 e A˚. The
where ECI and EC are the electrostatic energy of the
channel with and without the ion, respectively, and EI is
the electrostatic energy of the ion in bulk water.13 The
total electrostatic free energy for the passage of a probe
cation through the GlpF pore was computed by solving
the linearized Poisson–Boltzmann equation. The PBEQ
module79 as implemented in CHARMM,52 version 28,1.40 A˚, as described.42
The PMF profile was computed from equilibrium
distributions of x obtained from umbrella sampling
simulations. A harmonic biasing potential energy func-
tion was imposed on x to force the sampling over the
entire pore region, from K6.5 A˚ to 14.5 A˚.42 To this end,
20 ps of equilibration followed by 60 ps of production
were generated for each of 43 consecutive windows
spaced by 0.5 A˚. The reaction coordinate, x, was recorded
every 2 fs. The total simulation time required to build the
PMF profile for HC translocation was 3.44 ns.
Continuum electrostatic calculations
The electrostatic potential f(r) is given by the linearized
Poisson–Boltzmann equation:
V3ðrÞVfðrÞKK2DHðrÞfðrÞZK4prðrÞ
where 3(r), r(r), r and KDH are the space-dependent
dielectric constant, macromolecular charge density, pos-
ition vector, and modified Debye–Hu¨ckel screening
factor, respectively. In the absence of mobile ionic charges
in the solution, f(r) can be calculated by solving Poisson’s
equation:74–78
V3ðrÞVfðrÞZK4prðrÞ
The electrostatic free energy of transfer of a probe
cation from the aqueous solution to the center of the pore
of the channel is defined as:
EZ ðECIKECKEIÞ=2PMF profile was computed by unbiasing the windows
and combining them with the weighted histogram
analysis method (WHAM).72,73 The total simulation
time used to build the PMF profile for the reorientation
of the single-file water chain was 3.28 ns.
Proton hopping
The reversible thermodynamic work or PMF for the
translocation of an excess proton along the single-file
water chain was computed from MD simulations using a
continuous and derivable reaction coordinate:
x ¼
X
Oi
WOi
!K1 X
Oi
zOiWOi
!
In this equation:
WOi Z
X
Hj
fswðrOiHjÞK2
is a weighting function, and:
fswðrOiHjÞZ 1=ð1Cexp½frOiHjKrswg=dsw
Þ
is a switching function, where Oi, Hj are oxygen and
hydrogen atoms of the single-file water chain, respect-
507was used together with a set of optimized atomic radii for
amino acids.80 All explicit water molecules were removed
from the molecular system. The free energy of ion transfer

reflects the charge distribution of the protein) and the
reaction field (which is governed by the low-dielectric
cell width of 1.0 A˚, with subsequent focusing using a
0.5 A˚ cell width. A slab perpendicular to the z-axis
represented the low-dielectric region of the membrane. It
was centered at zZ5 A˚ and chosen to be of 40 A˚ thickness.
We used dielectric constants of 2, 2, and 80 for protein (3p),
membrane (3m), and water (3w), respectively.42 We
assigned a high dielectric constant (3cZ80) to the interior
of the pore. Earlier results showed that the choice of
membrane dielectric had only a moderate effect on the
reaction field.42 The calculations were repeated for five
variants in which the coulombic interactions between the
lumen contents and specified subsets of atoms were
turned off as well as in two single-point mutant channels,
as listed in Table 1. The calculation of the static field for
the wild-type channel was repeated in the presence of
150 mM ionic strength and of a transmembrane voltage of
K65 mV. Finally, to determine the protonation state of the
two single-point mutants N68D and N203D at physio-
logical pH, we carried out pKA calculations in the
Poisson–Boltzmann framework47 for Asp68 and Asp203,
respectively. Results obtained with membrane and
protein dielectric constants varying between 1, 2 and 4
show that the pKA of these two side-chains are within
1 pKA unit from that of their aqueous reference state
(3.9 for Asp). The largest shifts in pKA,K0.78 and 0.80 for
N68D and N203D, respectively, were obtained with
protein and membrane dielectric constant set to 1
(3pZ1Z3m). As a consequence, both side-chains were
modeled in their deprotonated (carboxylate) state.
Acknowledgements
We gratefully acknowledge the Canadian Insti-
tutes of Health Research (operating grant
MOP43949) and the Ontario Center for Genomics
Computing for support. R.P. is a CRCP Chairholder.
References
1. Hille, B. (2001). Ionic Channels of Excitable Membranes
3rd edit., Sinauer Associates, Inc., Sunderland, MA.
2. Doyle, D. A., Morais-Cabral, J. H., Pfuetzner, R. A.,boundaries). The radius of the probe cation was chosen as
that of sodium (rZ1.66 A˚)80 because the hydration free
energy of that ion, K98 kcal/mol,81 is comparable to that
of hydronium, K102 kcal/mol82 to K104 kcal/mol.83 A
comparable effective radius of 1.6 A˚ was used in another
study based on the Born solvation model.41 Our earlier
study showed that the reaction field profile exhibits only a
weak dependence on the probe radius, except at the
selectivity filter.42 The linearized Poisson–Boltzmann
equation was solved on a cubic grid of 151 cells with a(ESP) was computed with a cation placed at successive
positions along the channel pore.42
In the continuum approximation, the total electrostatic
free energy of the cation in the pore can be expressed as
the sum of contributions from the static field (which
508Kuo, A., Gulbis, J. M., Cohen, S. L. et al. (1998). The
structure of the potassium channel: molecular basis of
KC conduction and selectivity. Science, 280, 69–76.3. Dutzler, R., Campbell, E. B., Cadene, M., Chait, B. T. &
Mackinnon, R. (2002). X-ray structure of a CIC
chloride channel at 3.0 A˚ reveals the molecular basis
of anion selectivity. Nature, 415, 287–294.
4. Cramer, W. A. & Knaff, D. B. (1991). Energy Transduc-
tion in Biological Membranes. A Textbook of Bioenergetics,
Springer, New York.
5. Mitchell, P. (1961). Coupling of phosphorylation to
electron and hydrogen transfer by a chemi-osmotic
type of mechanism. Nature, 191, 144–148.
6. Luecke, H., Schobert, B., Richter, H. T., Cartailler, J.-P.
& Lanyi, J. K. (1999). Structural changes in bacterio-
rhodopsin during ion transport at 2 A˚ resolution.
Science, 286, 255–260.
7. Lanyi, J. K. (1999). Bacteriorhodopsin. Int. Rev. Cytol.
187, 161–202.
8. Arseniev, A. S., Barsukov, I. L., Bystrov, V. F.,
Lomize, A. L. & Ovchinnikov, Y. A. (1985). 1H NMR
study of gramicidin-A channel: head-to-head
right-handed single stranded helices. FEBS Letters,
186, 168–174.
9. Ketchem, R. R., Roux, B. & Cross, T. A. (1997). High-
resolution refinement of a solid-state NMR-derived
structure of gramicidin A in a lipid bilayer environ-
ment. Structure, 5, 11655–11669.
10. Roux, B. (2002). Theoretical and computational
models of ion channels. Curr. Opin. Struct. Biol. 12,
182–189.
11. Tian, F. & Cross, T. A. (1999). Cation transport: an
example of structural based selectivity. J. Mol. Biol.
285, 1993–2003.
12. Roux, B. (2002). Computational studies of the
gramicidin channel. Accts Chem. Res. 35, 366–375.
13. Roux, B. & Mackinnon, R. (1999). The cavity and pore
helices in KcsA KC channel: electrostatic stabilization
of monovalent cations. Science, 285, 100–102.
14. Gowen, J. A., Markham, J. C., Morrison, S. E., Cross,
T. A., Busath, D. D., Mapes, E. J. & Schumaker, M. F.
(2002). The role of Trp side chains in tuning single
proton conduction through gramicidin channels.
Biophys. J. 83, 880–898.
15. Hynes, J. T. (1999). The protean proton in water.
Nature, 397, 565–566.
16. Tuckerman, M., Laasonen, K., Sprik, M. & Parrinello,
M. (1995). Ab initio molecular dynamics simulation of
the solvation and transport of H3OC and OH- ions in
water. J. Phys. Chem. 99, 5749–5752.
17. Vuilleumier, R. & Borgis, D. (1999). Transport and
spectroscopy of the hydrated proton: a molecular
dynamics study. J. Chem. Phys. 111, 4251–4266.
18. Schmitt, U. W. & Voth, G. A. (1999). The computer
simulation of proton transport in water. J. Chem. Phys.
111, 9361–9381.
19. Marx, D., Tuckerman, M. E., Hutter, J. & Parrinello, M.
(1999). The natures of the hydrated excess proton in
water. Nature, 397, 601–604.
20. Pome`s, R. & Roux, B. (1995). Quantum effects on the
structure and energy of a protonated linear chain of
hydrogen-bonded water molecules. Chem. Phys.
Letters, 234, 416–424.
21. Pome`s, R. & Roux, B. (1996). Structure and dynamics
of a proton wire: a theoretical study of HC transloca-
tion along the single-file water chain in the gramicidin
channel. Biophys. J. 71, 19–39.
Proton Blockage in Aquaporins22. Pome`s, R. & Roux, B. (1998). Free energy profiles for
HC conduction along hydrogen-bonded chains of
water molecules. Biophys. J. 75, 33–40.

23. Mei, H. S., Tuckerman, M. E., Sagnella, D. E. & Klein,
M. L. (1998). Quantum nuclear ab initio molecular
dynamics study of water wires. J. Phys. Chem. B, 102,
10446–10458.
24. Brewer, M. L., Schmitt, U. W. & Voth, G. A. (2001). The
formation and dynamics of proton wires in channel
environments. Biophys. J. 80, 1691–1702.
25. Dellago, C., Naor, M. M. & Hummer, G. (2003). Proton
transport through water-filled carbon nanotubes.
Phys. Rev. Letters, 90, 105902/1–105902/4.
26. Pome`s, R. & Roux, B. (2002). Molecular mechanism of
HC conduction in the single-file water chain of the
gramicidin channel. Biophys. J. 82, 2304–2316.
27. Yu, C.-H. & Pome`s, R. (2003). Functional dynamics
of ion channels: modulation of proton movement
by conformational switches. J. Am. Chem. Soc. 125,
13890–13894.
28. Nagle, J. F. & Morowitz, H. J. (1978). Molecular
mechanisms for proton transport in membranes. Proc.
Natl Acad. Sci. USA, 75, 298–302.
29. Borgnia, M. J., Nielsen, S., Engel, A. & Agre, P. (1999).
Cellular and molecular biology of the aquaporin
water channels. Annu. Rev. Biochem. 68, 425–458.
30. Heymann, J. B. & Engel, A. (1999). Aquaporins:
phylogeny, structure and physiology of water chan-
nels. News Physiol. Sci. 14, 187–194.
31. Fujiyoshi, Y., Mitsuoka, K., de Groot, B. L., Philippsen,
A., Grubmu¨ller, H., Agre, P. & Engel, A. (2002).
Structure and function of water channels. Curr. Opin.
Struct. Biol. 12, 509–515.
32. Agre, P. & Kozono, D. (2003). Aquaporin water
channels: molecular mechanism for human diseases.
FEBS Letters, 555, 72–78.
33. Murata, K., Mitsuoka, K., Hirai, T., Walz, P., Agre, J. B.,
Heymann, A. et al. (2000). Structural determinants of
water permeation through aquaporin-1. Nature, 407,
599–605.
34. Sui, H., Han, B. G., Lee, J. K., Walian, P. & Jap, B. K.
(2001). Structural basis of water-specific transport
through the AQP1 water channel. Nature, 414,
872–878.
35. Fu, D., Libson, A., Miercke, L. J., Weitzman, C.,
Nollert, P., Krucinski, J. & Stroud, R. M. (2000). The
structure of a glycerol-conducting channel reveals the
basis for its selectivity. Science, 290, 481–486.
36. Savage, D. M., Egea, P. F., Robles-Colmenares, Y.,
O’Connell, J. D., III & Stroud, R. M. (2003). Archi-
tecture and selectivity in aquaporins 2.5 A˚ X-ray
structure of aquaporin. Z. PLoS Biol. 1, 334–340.
37. de Groot, B. L. & Grubmu¨ller, H. (2001). Water
permeation across biological membranes: mechanism
and dynamics of aquaporin-1 and GlpF. Science, 294,
2353–2357.
38. Tajkhorshid, E., Nollert, P., Jensen, M., Miercke, J. W.,
O’Connell, J. L., Stroud, R. M. & Schulten, K. (2002).
Control of the selectivity of the aquaporin water
channel family by global orientational tuning. Science,
296, 525–530.
39. Jensen, M.., Tajkhorshid, E. & Schulten, K. (2003).
Electrostatic tuning of permeation and selectivity in
aquaporin water channels. Biophys. J. 85, 2884–2899.
40. de Groot, B. L., Frigato, T., Helms, V. & Grubmu¨ller, H.
(2003). The mechanism of proton exclusion in the
aquaporin-1 water channel. J. Mol. Biol. 333, 279–293.
41. Burykin, A. & Warshel, A. (2003). What really
Proton Blockage in Aquaporinsprevents proton transport through aquaporin?
Charge self-energy versus proton wire proposals.
Biophys. J. 85, 3696–3706.42. Chakrabarti, N., Tajkhorshid, E., Roux, B. & Pome`s, R.
(2004). Molecular basis of proton blockage in aqua-
porins. Structure, 12, 65–74.
43. Ilan, B., Tajkhorshid, E., Schulten, K. & Voth, G. A.
(2004). The mechanism of proton exclusion in aqua-
porin channels. Proteins: Struct. Funct. Bioinform. 55,
223–228.
44. Roux, B., Berne`che, S. & Im, W. (2000). Ion channels,
permeation, and electrostatics: insight into the func-
tion of KcsA. Biochemistry, 39, 13295–13306.
45. Parsegian, V. A. (1969). Energy of an ion crossing a
low dielectric membrane: solutions to four relevant
electrostatic problems. Nature, 221, 844–846.
46. Chiu, S. W., Subramaniam, S. & Jakobsson, E. (1999).
Simulation study of gramicidin/lipid bilayer system
in excess water and lipid. II. Rates and mechanism of
water transport. Biophys. J. 76, 1939–1950.
47. Faraldo-Gomez, J. D. & Roux, B. (2004). Electrostatics
of ion stabilization in a ClC chloride channel
homologue from Escherichia coli. J. Mol. Biol. 339,
981–1000.
48. Eisenberg, B. (2003). Why can’t ions move through
water channels? Biophys. J. 85, 3427–3428.
49. Lanyi, J. K. (2004). Bacteriorhodopsin. Annu. Rev.
Physiol. 66, 665–688.
50. Onufriev, A., Smondyrev, A. & Bashford, D. (2003).
Affinity changes driving unidirectional proton trans-
port in the bacteriorhodopsin photocycle. J. Mol. Biol.
332, 1183–1193.
51. Pome`s, R. & Yu, C.-H. (2003). Relay and blockage or
protons in water chains. Frontiers Biosci. 8, d1288–
d1297.
52. Brooks, B. R., Bruccoleri, R. E., Olafson, B. D., States,
D. J., Swaminathan, S. & Karplus, M. (1983).
CHARMM: a program for macromolecular energy
minimization and dynamics calculations. J. Comput.
Chem. 4, 187–217.
53. MacKerell, A. D., Jr, Bashford, D., Bellott, M.,
Dunbrack, R. L., Jr, Evanseck, J. D., Field, M. J. et al.
(1998). All-atom empirical potential for molecular
modeling and dynamics studies of proteins. J. Phys.
Chem. B, 102, 3586–3616.
54. Jorgensen, W. L., Chandrasekhar, J., Madura, J. D.,
Impey, R. W. & Klein, M. L. (1983). Comparison of
simple potential functions for simulating liquid water.
J. Chem. Phys. 79, 926–935.
55. Stillinger, F. H. & David, C. W. (1978). Polarization
model for water and its ionic dissociation products.
J. Chem. Phys. 69, 1473–1484.
56. Stillinger, F. H. (1979). Dynamics and ensemble
averages for the polarization models of molecular
interactions. J. Chem. Phys. 71, 1647–1651.
57. Weber, T. A. & Stillinger, F. H. (1982). Reactive
collisions of H3O
C and OHK studied with the
polarization model. J. Phys. Chem. 8, 1314–1318.
58. Wo´jcik, M. J. & Rice, S. A. (1986). The vibrational
spectrum of the water dimer: some model based
predictions. J. Chem. Phys. 84, 3042–3048.
59. Kobayashi, C., Iwahashi, K., Saito, S. & Ohmine, I.
(1996). Dynamics of proton attachment to water
cluster: proton transfer, evaporation, and relaxation.
J. Chem. Phys. 105, 6358–6366.
60. Ojama¨e, L., Shavitt, I. & Singer, S. J. (1998). Potential
models for simulations of the solvated proton in
water. J. Chem. Phys. 109, 5547–5564.
50961. Gomez, M. A. & Pratt, L. R. (1998). Construction of
simulation wave functions for aqueous species:
D3O
C. J. Chem. Phys. 109, 8783–8789.

62. Hammes-Schiffer, S. (1998). Mixed quantum/classical
dynamics of hydrogen transfer reactions. J. Phys.
Chem. A, 102, 10443–10454.
63. Webb, S. P. & Hammes-Schiffer, S. (2000). Fourier grid
Hamiltonian multi-configurational self-consistent-
field: a method to calculate multidimensional hydro-
gen vibrational wavefunctions. J. Chem. Phys. 113,
5214–5227.
64. Schumaker, M. F., Pome`s, R. & Roux, B. (2000). A
combined molecular dynamics diffusion model of
single proton conduction through gramicidin. Bio-
phys. J. 79, 2840–2857.
65. Schumaker, M. F., Pome`s, R. & Roux, B. (2001).
Framework model for single proton conduction
through gramicidin. Biophys. J. 80, 12–30.
66. Chernyshev, A. & Cukierman, S. (2002). Thermodyn-
amic view of activation energies of proton transfer in
various gramicidin A channels. Biophys. J. 82, 182–192.
solutions and in enzymes. J. Am. Chem. Soc. 102, 6218–
6226.
empirical valence bond model for describing proton
transfer in HC(H O) clusters and liquid water. Chem.
histogram analysis method for free-energy calcu-
lations in biomolecules. 1. The method. J. Comput.
Chem. 13, 1011–1021.
73. Roux, B. (1995). The calculation of the potential of
mean force using computer simulations. Comput.
Phys. Commun. 91, 275–282.
74. Warshel, A. & Russell, S. T. (1984). Calculations of
electrostatic interactions in biological systems and in
solutions. Quart. Rev. Biophys. 17, 283–422.
75. Klapper, I., Hagstrom, R., Fine, R., Sharp, K. &
Honig, B. (1986). Focusing of electric fields in the
active site of Cu–Zn superoxide dismutase: effects of
ionic strength and amino-acid modification. Proteins:
Struct. Funct. Genet. 1, 47–59.
76. Davis, M. E. & McCammon, J. A. (1990). Electrostatics
in biomolecular structure and dynamics. Chem. Rev.
90, 509–521.
77. Levitt, D. (1991). General continuum theory for a
in biology and chemistry. Science, 268, 1144–1149.
solutions to the Poisson–Boltzmann equation.
in
510 Proton Blockage in Aquaporins2 n
Phys. Letters, 284, 71–77.
69. Day, T. J. F., Soudackov, A. V., Cˇuma, M., Schmitt,
U. W. & Voth, G. A. (2002). A second generation
multistate empirical valence bond model for proton
transport in aqueous systems. J. Chem. Phys. 117,
5839–5849.
70. Torrie, G. M. & Valleau, J. P. (1974). Monte–Carlo free
energy estimates using non-Boltzmann sampling—
application to subcritical Lennard–Jones fluid. Chem.
Phys. Letters, 28, 578–581.
71. Torrie, G. M. & Valleau, J. P. (1977). Non-physical
sampling distributions in Monte–Carlo free energy
estimation: umbrella Sampling. J. Comput. Phys. 23,
187–199.
72. Kumar, S., Bouzida, D., Swendsen, R. H., Kollman,
P. A. & Rosenberg, J. M. (1992). The weighted
(Received 28 May 2004; receivedComput. Phys. Commun. 111, 59–75.
80. Nina, M., Beglov, D. & Roux, B. (1997). Atomic Born
radii for continuum electrostatic calculations based on
molecular dynamics free energy simulations. J. Phys.
Chem. B, 101, 5239–5248.
81. Burgess, J. (1978). Metal Ions in Solution, Wiley,
New York.
82. Pearson, R. G. (1986). Ionization potentials and
electron affinities in aqueous solution. J. Am. Chem.
Soc. 108, 6109–6114.
83. Chambers, C. C., Hawkins, G. D., Cramer, C. J. &
Truhlar, D. G. (1996). Model for aqueous solvation
based on class IV atomic charges and first solvation
shell effects. J. Phys. Chem. 100, 16385–16398.
84. Humphrey, W., Dalke, A. & Schulten, K. (1996). VMD:
visual molecular dynamics. J. Mol. Graph. 14, 33–38.
Edited by G. von Heijne
revised form 6 August 2004; accepted 11 August 2004)68. Vuilleumier, R. & Borgis, D. (1998). An extended 79. Im, W., Beglov, D. & Roux, B. (1998). Continuum
solvation model: electrostatic forces from numerical67. Warshel, A. M. & Weiss, R. M. (1980). An empirical
valence bond approach for comparing reactions in
multiion channel. Biophys. J. 59, 271–277.
78. Honig, B. & Nichols, A. (1995). Classical electrostatics