Structure, Vol. 12, 65–74, January, 2004, 2004 Elsevier Science Ltd. All rights reserved. DOI 10.1016/j .str .2003.11.017
Molecular Basis of Proton
Blockage in Aquaporins
and their impaired function may lead to pathological
situations, such as nephrogenic diabetes insipidus and
congenital cataract of the eye (Borgnia et al., 1999; Ko-
Nilmadhab Chakrabarti,1,2 Emad Tajkhorshid,3
Benoıˆt Roux,4 and Re´gis Pome`s1,2,*
1Structural Biology and Biochemistry
Hospital for Sick Children zono et al., 2002). All members of the AQP superfamily
act as passive transporters of water. Aquaglyceroporins,555 University Avenue
Toronto, Ontario M5G 1X8 a subfamily of AQPs, have the ability to conduct small,
linear polyols stereoselectively in addition to water mol-Canada
2 Department of Biochemistry ecules (Fu et al., 2000; Borgnia et al., 1999; Jensen et
al., 2002; Unger, 2000). The most studied aquaglycero-University of Toronto
Toronto, Ontario porin is the E. coli glycerol uptake facilitator GlpF, a
channel used by organisms to absorb glycerol from theCanada
3 Theoretical and Computational Biophysics Group environment for use in metabolism when there is little
glycerol available.Beckman Institute
University of Illinois at Urbana-Champaign Atomic resolution models of human aquaporin-1
(AQP1) based on electron microscopy revealed the tet-405 North Mathews
Urbana, Illinois 61801 rameric architecture and folding of AQPs (Murata et al.,
2000; Ren et al., 2001). GlpF was the first AQP whose4 Department of Biochemistry
Weill Medical College of Cornell University high-resolution X-ray crystallographic structure (2.2 A˚)
was determined (Fu et al., 2000). GlpF and AQP1 exhibit1300 York Avenue
New York, New York 10021 a remarkable structural similarity (Unger, 2000). The first
GlpF structure included three glycerol molecules, a nat-
ural substrate for GlpF, in each monomer. Molecular
dynamics (MD) simulation of glycerol-saturated GlpFSummary
shed light into the details of protein-substrate interac-
tion and conduction pathway (Jensen et al., 2001). InWater transport channels in membrane proteins of
each monomer, the transmembrane segment of thethe aquaporin superfamily are impermeable to ions,
channel is formed by six helices and two half-membraneincluding H
and OH
. We examine the molecular basis
spanning loops that meet each other at the center offor the blockage of proton translocation through the
the channel. These two loops are related by quasi 2-foldsingle-file water chain in the pore of a bacterial aqua-
symmetry and are structurally and functionally importantporin, GlpF. We compute the reversible thermody-
(Murata et al., 2000; Fu et al., 2000; Ren et al., 2001;namic work for the two complementary steps of the
Zhu et al., 2001; Jensen et al., 2001; Sui et al., 2001).Grotthuss “hop-and-turn” relay mechanism: consecu-
Half of each loop is nonhelical and defines a curvilineartive transfers of H
along the hydrogen-bonded chain
conduction pathway constituted by backbone carbonyl(hop) and conformational reorganization of the chain
groups that are exposed into the channel interior (Jen-(turn). In the absence of H
, the strong preference for
sen et al., 2001). The other half of each loop (known asthe bipolar orientation of water around the two Asn-
M3 and M7 in GlpF, HB and HE in AQP1) is  helical.Pro-Ala (NPA) motifs lining the pore over both unidirec-
These helices terminate at the channel center, wheretional polarization states of the chain precludes the
two conserved NPA (Asn-Pro-Ala) motifs meet eachreorganization of the hydrogen-bonded network. In-
other (Fu et al., 2000; Murata et al., 2000; Sui et al.,versely, translocation of an excess proton in either
2001). Internal hydrogen bonds within the NPA motifsdirection is opposed by a free-energy barrier centered
lead to a configuration in which one of the amide hydro-at the NPA region. Both hop and turn steps of proton
gen atoms of each asparagine points toward the pore.translocation are opposed by the electrostatic field of
Together, these two NH2 groups provide two hydrogen-the channel.
bonding sites for the permeants at the center of the
channel. Three-dimensional structures of the glycerol-
Introduction free GlpF channel have also been solved by X-ray crys-
tallography (Tajkhorshid et al., 2002), revealing the pres-
Aquaporins (AQPs) are membrane water channels pres- ence of a single file of water molecules. Recent MD
ent in all life forms, and more than 100 of them have studies have addressed the molecular mechanism of
been characterized (Agre et al., 1998; Borgnia et al., water transport in AQP1 and GlpF (de Groot and Grub-
1999) since their discovery more than a decade ago. mu¨ller, 2001; Tajkhorshid et al., 2002; Zhu et al., 2002;
Aquaporin-1 was first identified as an integral membrane Jensen et al., 2003).
protein in red blood cells and renal proximal tubules AQPs transport water efficiently, but they are imper-
(Denker et al., 1988, Preston and Agre, 1991), where it meable to charged species. Such selectivity is of crucial
functions as a water-selective membrane pore (Preston importance to the proper function of cell membranes,
et al., 1992). Eleven human AQPs have been identified, as the flow of water in response to osmotic stress must
not disturb the electrochemical properties of the mem-
brane. Blockage of protons is essential in maintaining*Correspondence: pomes@sickkids.ca

Structure
66
the transmembrane proton gradient in the cell (Mitchell,
1961). However, this property is somewhat unexpected
given the presence of a chain of water molecules in the
pore. Water chains embedded in biological channels are
known to mediate the long-range movement of H via a
hop-and-turn Grotthuss relay mechanism (de Grotthuss,
1806; Nagle and Morowitz, 1978; Agmon 1995; Pome`s
and Roux, 1996, 1998, 2002). In this process of structural
diffusion, the excess proton need not diffuse through
the pore. Instead, successive transfers (hops) exchange
hydrogen nuclei along hydrogen bonds forming an ex-
tended network (proton wire); because this process in-
verts the polarity of proton-relaying groups, proton hop-
ping is complemented by the reorientation (turn) of each
group in the chain (Pome`s and Roux, 1998, 2002). To-
gether, these two steps, which correspond respectively
to the propagation of an ionic and a bonding defect,
result in the long-range displacement of H in proteins
(Nagle and Morowitz, 1978; Pome`s and Roux, 1996).
Detailed theoretical studies of water wires embedded in
native and chemically modified forms of the gramicidin
channel have uncovered the balance of factors at play
in the molecular mechanism of proton translocation
(Pome`s and Roux, 1996, 2002; Yu and Pome`s, 2003).
In particular, these studies have shown how structural
fluctuations of the hydrogen-bonded network formed
by the water chain and the channel mediate the rapid, Figure 1. The Molecular System Simulated in the Present Study
passive transport of H. The monomeric GlpF channel is shown in -carbon trace. Water
In the present work, we study the properties of the molecules on either side of the pore are shown, and the nine water
molecules in the 20 A˚ long single-file pore (
6  z  14 A˚) arehydrogen-bonded network in the single-file region of
highlighted in space-filling representation. Also highlighted are sa-AQPs in order to determine the molecular basis for the
lient features of the lumen: Arg206 at the selectivity filter in theexclusion of protons. Previous studies of the water con-
periplasmic entry, the two Asn side chains of the conserved NPAduction pathway have been used to infer two possible motifs, carbonyl groups protruding into the pore, and  helices M3
mechanisms of proton blockage, respectively by control and M7. This and the following depictions of molecular structures
of the orientation of the water chain and by interruption were generated with the VMD program (Humphrey et al., 1996).
of the connectivity of the hydrogen-bonded network.
MD simulations have shown that permeating water mol-
hydrogen-bonded network in the single-file region of theecules adopt a preferential orientation in the single-file
GlpF channel successively with and without an excessregion of the GlpF pore (Tajkhorshid et al., 2002). One
proton. The reversible work for the turn step of structuralwater molecule accepts two hydrogen bonds from the
diffusion is computed by modulating the arrangementtwo NPA motifs and orients perpendicularly to the chan-
of water dipoles in the GlpF pore in the absence of annel axis, and the water molecules on either side form
excess proton. We then characterize the hop step ofpolarized hydrogen-bonded chains extending in oppo-
the Grotthuss mechanism by computing the potentialsite directions (Tajkhorshid et al., 2002). This bipolar
of mean-force for forcing an excess proton throughorganization was ascribed to the opposed dipoles of
the GlpF pore. Analysis of the free energy profiles forthe M3 and M7 helices (Murata et al., 2000; Tajkhorshid
the translocation of bonding and ionic defects in theet al., 2002) as well as to hydrogen-bond donation by
water chain indicates that proton blockage is essentiallythe NPA motifs (Tajkhorshid et al., 2002). Because the
achieved by the electrostatic field of the channel.bipolar arrangement, with OH bonds pointing outwards,
is incompatible with the uptake of an excess proton into
Resultsthe pore from either end, it was proposed that proton
blockage results from orientational control (Tajkhorshid
The simulated molecular system is shown in Figure 1.et al., 2002). De Groot et al. (2001) postulated that proton
The 20 A˚ long constriction region of the GlpF monomerexclusion arises from the disruption of the single-file
accommodates nine water molecules in single-file ar-water chain in the narrowest part of the pore, which is
rangement. Hydrogen-bonding partners of these waterformed by an arginine and by two aromatic residues at
molecules include the backbone carbonyl groups (CO)the periplasmic mouth. However, the translocation of a
of residues 64–67 and 199–202, which protrude intobonding defect has not been considered, and neither
the channel. The other hydrogen-bonding groups liningof these two postulates has been tested in the presence
the channel consist of the positively charged residueof an excess proton in the lumen.
Arg206, near the periplasmic mouth, and the two con-To gain a better understanding of the thermodynamic
served asparagines (Asn68 and Asn203) of the NPA mo-basis of proton exclusion in AQPs, we use molecular
simulations to examine the structural fluctuations of the tifs (Murata et al., 2000; Fu et al., 2000; Sui et al., 2001;

Proton Blockage in Aquaporins
67
et al., 2002). Both O and H distributions obtained from
our finite-size model reproduce the preferred localiza-
tion of water atoms obtained by Tajkhorshid et al. (2002).
These preferred positions result from hydrogen bonding
with the channel. Protruding carbonyl groups stabilize
permeating water molecules along the conduction path-
way. The amide groups of Asn68 and Asn203 donate
one hydrogen bond each to the sixth water molecule
from the periplasm. The bipolar orientation of water mol-
ecules identified in previous studies (Tajkhorshid et al.,
2002; de Groot and Grubmu¨ller 2001) is reflected in the
relative arrangement of O and H atoms. In particular,
the symmetric distribution of two H atoms around z
6 A˚
confirms that water 6 remains oriented perpendicular
to the channel. Dynamic fluctuations are in very good
agreement except in the region of the selectivity filter.
In both simulations, a gap in O peaks 2 and 3 at z

2 A˚
points to a resilient interruption in the hydrogen-bonded
water chain. However, water molecules near the peri-
plasmic end of the pore are somewhat less mobile in
the present simulations, which reflects the fact that no
net permeation of water was allowed by construction.
The distribution of O atoms obtained in the absence of
restraints (Tajkhorshid et al., 2002) suggests that water
permeation gives rise to the occasional migration of this
hydrogen-bonding defect.
The potential of mean force (PMF) for the turn step
Figure 2. Equilibrium Distribution of Oxygen and Hydrogen Atoms of the Grotthuss mechanism is depicted in Figure 3 to-
of the Single-File Chain of Water Molecules in the Channel Pore gether with three representative snapshots of the water
The results obtained in the present study (red) are compared with chain. The PMF profiles computed successively with
the results obtained from a simulation by Tajkhorshid et al. (2002) the TIP3P and PM6 force fields of water are in very good
(green). Top: representative conformation of the nine water mole-
qualitative and quantitative agreement with each othercules in the single file, shown with their hydrogen-bonding partners
despite fundamental differences between the two waterlining the channel; middle: relative atomic density of oxygens, O(z);
bottom: relative atomic density of hydrogens, H(z). In this and sub- models. The free energy profile for the reorientation of
sequent figures, the location of the selectivity filter (
5  z 
1 A˚) water dipoles in the pore indicates that the preference
and of the NPA motifs (4.75  z  7 A˚) are indicated by pink and for the bipolar organization of the water chain is very
blue bars, respectively. strong. At z

1 e·A˚, the bipolar conformation (B) is
overwhelmingly favored over the two metastable states
of the chain, A and C (respectively at z 
8 and
Nollert et al., 2001). Together with Arg206, two aromatic 8 e·A˚). The latter two conformers correspond to fully
residues (Trp48 and Phe200) form the selectivity filter. polarized chains, with water dipole moments pointing
This bottleneck has been proposed to work as a size- respectively toward the periplasmic and cytoplasmic
exclusion filter and was observed to disrupt the hydro- vestibules. The free energy barrier to move the bonding
gen-bonded water chain (de Groot et al., 2001). The defect from water 6 to the periplasm (B→C transition)
other groups lining the channel (data not shown) are is 12–16 kcal/mol with PM6 and TIP3P water models,
hydrophobic. The amphipathic nature of the pore along whereas that for migration to the cytoplasmic vestibule
the conduction pathway is of crucial importance to the (B→A transition) is 7–8 kcal/mol. These two transloca-
orientational control of water and glycerol permeants in tion processes correspond to the reorientation of water
the lumen (Tajkhorshid et al., 2002). Two half helices, molecules 1–6 and 6–9, respectively.
M3 and M7, were shown to contribute to the bipolar The PMF for the translocation of an excess proton
orientation of the water chain around the NPA region (hop) is depicted in Figure 4. The top four panels illus-
(Tajkhorshid et al., 2002). trate representative configurations of the protonated
The structure and fluctuations of the water chain ob- water chain as H is forced from the periplasmic mouth
tained in the absence of an excess proton confirm the of the channel toward the cytoplasmic mouth. The free
results of previous simulation studies. The equilibrium energy profile consists of a barrier that starts at 
2 A˚
distribution of the nine single-file water molecules is and peaks to 4.5 kcal/mol at the NPA site (  6 A˚). The
shown in Figure 2. The top panel depicts the preferred corresponding snapshot (configuration III) shows that,
orientation of the water chain, whereas the middle and among all nine single-file water molecules, water 6 is
bottom panels show the distribution of oxygen and hy- least likely to host the excess proton. Accordingly, water
drogen atoms, respectively. Overall, the results com- molecules 4–7 are pulled close to each other, reflecting
puted from a 2 ns simulation with the present model are the relative delocalization of the excess proton in the
in good agreement with those reported from a simulation NPA region. This may be due in part to the lack of
hydrogen-bond acceptors at the NPA site, which makesof the channel tetramer in a lipid membrane (Tajkhorshid

Structure
68
Figure 3. Preferential Organization of the Water Chain in the Ab-
sence of an Excess Proton
(Top) Three representative conformations of the hydrogen-bonded
water chain in the GlpF pore. (A) Fully-polarized water chain with
OH bonds oriented toward the periplasmic entry (z 
8 e·A˚); (B)
preferred bipolar arrangement of the water chain (z 
1e·A˚), with
a bonding defect on water 6 at the NPA site; (C) fully-polarized water
chain with OH bonds oriented toward the cytoplasmic entry (z 
7.5 e·A˚). (Bottom) Potential of mean force (PMF) for the reorientation
of the single-file water chain obtained with water models TIP3P
(green) and PM6 (red). The primary minimum corresponds to confor-
mation B and the two secondary minima to conformations A and
C. In the case of the TIP3P model, z was scaled by a factor of
0.417, the fractional charge of H atoms in that potential (Jorgensen
et al., 1983), for consistency with the magnitude of z obtained with
the PM6 model, in which the formal charge of hydrogen is 1.
these water molecules relatively unfavorable proton
hosts (Pome`s and Roux, 2002). The relative flatness of
the 2.5 A˚ wide barrier top may also reflect the delocaliza-
Figure 4. Structure and Energy of the Protonated Water Chaintion of the excess proton. The proton translocation free
energy profile is also essentially flat in the region of the (Top) Representative conformations of the protonated water chain
in the GlpF pore, with the excess proton (I) around Arg206 at theselectivity filter
6    0 A˚). Thus, the proximity of
periplasmic mouth, (II) between the selectivity filter and the NPApositively charged Arg206 does not oppose proton
region, (III) around NPA region, and (IV) at the cytoplasmic mouth.translocation. Snapshots I and II, at  
5.5 and
2 A˚,
(Bottom) Potential of mean force (PMF) for the translocation of an
correspond to the excess proton hosted respectively by excess proton through the single-file chain of water molecules. The
water molecules 1 and 3, on either side of the selectivity free energy barrier peaks at the NPA site.
filter. Conformer IV, with the excess proton on water 8,
lies 2.5 kcal/mol below the barrier top. In both conforma-
tions II and IV, the water molecule hosting H is coordi- (Pome`s and Roux, 1996, 2002). Throughout the translo-
cation process, the polarization of all water moleculesnated by three hydrogen-bond acceptors, two from
neighboring single-file water molecules and one from a in the chain is dictated by the location of the excess
proton, with water dipoles pointing away from the centerprotruding channel carbonyl group (of Phe200 and
His66, respectively), mimicking the optimal solvation of of charge. Thus, when the excess proton is at the edge
of the constriction pore, the water chain is polarized,protonated water observed in bulk water (Agmon, 1995;
Tuckerman et al., 1995) and in the gramicidin channel similar to conformations A and C (see Figure 3) of the

Proton Blockage in Aquaporins
69
static field contribution to the ESP for the translocation
of a (point) charge through the pore also exhibits these
two features (Figure 5, bottom). Although the total elec-
trostatic barrier is dominated by effects arising from
dielectric boundaries (reaction field), the agreement be-
tween PMF, ESP, and static field profiles indicates that
the free energy for the movement of H in the single-
file region is essentially determined by the distribution
of charge of the channel. The continuing drops in total
electrostatic energy beyond the single-file region (below
z 
6 A˚ and above z  14 A˚) suggest that the electro-
static barrier opposing the passage of a cation in the
channel originates well outside the narrow pore region.
This barrier is not very sensitive to the size of the ionic
probe. The ESP profile is identical for a smaller probe
ion (radius of 1.4 A˚) and does not change appreciably
for a larger ion such as potassium (radius of 2.03 A˚),
except at the bottleneck formed by the selectivity filter,
where dielectric effects due to size exclusion result in
the emergence of an additional peak centered at Arg206.
Discussion
Electrostatic Control
We have computed the free energy profiles for both hop
and turn steps of the Grotthuss mechanism in the single-
file region of the GlpF channel. The two PMF profiles,
an asymmetric well for the turn, and an asymmetric
barrier for the hop, mirror each other in two important
respects. First, both extrema are located at the NPAFigure 5. Electrostatic Energy Profiles
site. Second, the free energy differences between the(Top) Electrostatic potential energy (ESP) for the movement of a
NPA site and the pore extremities vary in the same pro-probe cation (with a radius of 1.66 A˚) along the GlpF channel (red)
portion in both PMF profiles. At about 7 kcal/mol, theshown with the same profile obtained without the low-dielectric
membrane environment (green) and with the PMF for the transfer barrier for the translocation of a bonding defect from
of an excess proton (black) in the single-file region (see Figure 4). NPA to the cytoplasmic end of the pore is roughly half
The PMF profile was shifted for comparison. (Bottom) Effect of the the size of the barrier for migration from NPA to the
ionic probe radius on the ESP profile. The three red curves show
periplasmic extremity of the single-file region (14 2total ESP profiles obtained with r  1.66 A˚ (solid), 1.4 A˚ (dotted),
kcal/mol). Similarly, the barrier facing proton movementand 2.03 A˚ (dashed). The static field contribution to the ESP is
from the cytoplasmic end of the pore to NPA (2.5 kcal/shown in blue, and the PMF is again shown in black and shifted for
comparison. mol) is roughly half the size of that opposing the ionic
translocation from the periplasmic end (4.5 kcal/mol).
Together, these features strongly suggest that theunprotonated water chain; with H on water 6, the bipo-
forces opposing both hop and turn steps of the Grot-lar orientation of water molecules 1–5 and 7–9 is identical
thuss mechanism have the same physical origin. Theto that of conformation B.
overall agreement between the hop PMF and the energyTo gauge the effect of the electrostatic field on the
profile obtained from continuum electrostatic calcula-translocation of cations throughout the channel, we
tions (Figure 5) indicates that the free energy profilecomputed the electrostatic potential (ESP) energy pro-
opposing proton transfer trough the pore is essentiallyfile for the movement of a cationic probe along the pore
determined by electrostatic forces. Thus, the electro-in the continuum limit of the solvent (Figure 5). The Born
static field governs not only the structure and fluctua-radius of this ion (1.66 A˚) was chosen based on the
tions of the unprotonated water chain, but also the freehydration free energy of hydronium (see Experimental
energy barrier for an excess proton.Procedures). The ESP profile consists of a 11 kcal/mol
barrier centered at the NPA site (z  5.5 A˚) and drops
off on both sides of the narrow pore, reaching zero at Hierarchy of Physical Interactions
In GlpF, the strong electrostatic field of the channelz 
12 and 24 A˚. Neglect of the low-dielectric mem-
brane environment results in a 2 kcal/mol downshift of keeps the excess charge away from the NPA site and
forces the unprotonated chain into a bipolar organiza-the barrier but does not alter the shape of the ESP. In
the single-file region, the ESP for the movement of a tion. The above results suggest that the following hierar-
chy of molecular interactions is at play in the control ofpositive charge is similar to the PMF for the transfer of
an excess proton. Both profiles peak at the same loca- structural diffusion in the GlpF channel: ion-channel

ion-water
water-channel
water-water. The first twotion and feature a shoulder at the selectivity filter, about
4 kcal/mol below the barrier top (Figure 5, top). The inequalities are demonstrated by an analysis of the PMF

Structure
70
for the transfer of an excess proton. Even though ion- have severe shortcomings to represent the electrostatic
response of hydrogen-bonded water molecules in nar-water and water-channel interactions are optimized in
configuration III, where the bipolar organization of water row pores (Roux et al., 2000). Probable mitigating factors
in the present case include the relative rigidity of theprovides adequate solvation of the positive charge and
satisfies the preferred arrangement of the unprotonated GlpF channel (Jensen et al., 2001) and the fact that
overwhelming features governing ionic blockage, ratherwater chain (B), ion-channel interactions make it the
least likely conformation of the protonated chain. The than the subtle balance of effects involved in ionic solva-
tion and permeation through ion channels, dominate theprecedence of ion-water interactions over water-chan-
nel interactions is demonstrated in conformations I and present results. PB calculations also suffer from the
uncertain choice of dielectric constants for the mem-IV, where the intrusion of an excess proton at either end
of the single-file region is sufficient to polarize the water brane and the protein interior. The electrostatic gradi-
ents obtained with different dielectric representationschain fully despite the intrinsic preference for bipolar
conformation B over A and C in the absence of H. Water of the membrane (i.e., membrane  2 and 80) are in excellent
agreement in the pore region (see Figure 5), showingresponds to the ionic charge, whereas the comparatively
rigid of the channel does not, so that the ionic PMF is that the effect of the membrane low dielectric is moder-
ate in the narrow pore. Although this is only a roughdominated by interactions with the channel. Finally, the
third inequality is evidenced by the preference for the estimate, the magnitude of the overall electrostatic bar-
rier (over 10 kcal/mol) would be more than sufficient tobipolar orientation of the unprotonated water chain over
the fully polarized conformations, which have been block proton transport effectively at physiological pH.
shown to be intrinsically favored over all other states in
water chains confined in nonpolar environments (Pome`s
Protons versus Other Ionsand Roux, 1998). In summary, these results highlight the
The dominance of electrostatic forces suggests that theplasticity of the hydrogen-bonded water chain and show
physical basis for the exclusion of protons may alsothat the orientational polarization of water molecules
apply to other cations. However, while electrostatic gra-responds first to the ionic charge, then to the permanent
dients apply equally, in the first approximation, to ionscharge distribution of the channel, and last to fluctuating
of like charge, this effect may be expected to be modu-dipoles (i.e., other water molecules).
lated by differences in the fine coordination of ionic
species, particularly in single-file environments. Ion
channels compensate for the partial dehydration of ions,Limitations of the Model
The overall agreement between the ESP and PMF pro- and harness structural differences in the first hydration
shell of ions by mimicking the aqueous coordinationfiles lends confidence to the results of both MD simula-
tions and PB calculations, despite the limitations arising of permeating ions in narrow pores. Solvation of the
aqueous proton by water molecules involves tricoordi-from the approximations underlying these two ap-
proaches. Approximations in the MD calculations re- nation of protonated water molecules in hydronium or
Zundel cations by three hydrogen-bond acceptors (Ag-ported here include the use of empirical force fields
and the finite size of the model, together with restraints mon, 1995; Tuckerman et al., 1995). In the gA channel,
this very coordination is provided by one backbone car-precluding water diffusion and limiting the dynamic part
of the system to the close vicinity of the single-file pore. bonyl group and two neighboring water molecules to
each water O atom (Pome`s and Roux, 1996, 2002). Pro-However, the analysis of Figure 2 indicates that our use
of conformational restraints and of a finite-size model truding carbonyl groups of GlpF provide adequate tri-
coordination to water molecules in much of the single-only has a moderate effect on the structure and fluctua-
tions of the pore contents. In addition, the consistency file region of the pore, whereas hydrogen-bond donation
by Arg206 and the two Asn side chains of the NPA motifsbetween the properties of the water chain with the re-
sults of an earlier study of the channel tetramer in a compromise the ideal coordination of protonated water.
The latter factor probably explains the relative delocal-hydrated membrane (Tajkhorshid et al., 2002) validates
the neglect of lipid and of other monomers in the present ization of the excess proton near the NPA site, but does
not appear to contribute significantly to the free energystudy. In further support of our finite-size model, the
analysis of a simulation of the tetrameric GlpF in a lipid profile of proton translocation. Thus, based on the pres-
ent study, local solvation effects seem to play a second-membrane (Jensen et al., 2003) shows that electrostatic
interactions of water molecules in the single-file region ary role in the exclusion of H from the GlpF channel.
However, the hydration of alkali metal ions, with a largerwith the lipid molecules and the other monomers in the
tetramer are negligible. Moreover, limitations inherent number of water molecules in the first hydration shell,
is very different from that of protons. For example, eightto the potential energy function used to describe unpro-
tonated water are mitigated by the agreement between carbonyl groups stabilize K ions in the narrow selectiv-
ity filter of the KcsA potassium channel, providing “sur-the results obtained with very different force fields (Fig-
ures 2 and 3). Thus, our molecular model captures the rogate solvation” to the permeating cation (Morais-
Cabral et al., 2001). Thus, it is likely that the scarcity ofstructure and fluctuations of the water chain despite
differences in the size of the model and in the potential carbonyl groups, together with the absence of other
groups suitable to the solvation of a cation in the amphi-energy function. PB calculations are plagued by the
neglect of thermal fluctuations and the continuum de- pathic pore, would contribute much more significantly
to the relative destabilization of other cations. By thescription of the membrane and of the water solvent. The
continuum dielectric approximation has been shown to same token, in the case of anions it is likely that adverse

Proton Blockage in Aquaporins
71
solvation effects would take precedence over the elec- ion channel (Pome`s and Roux, 1996, 2002). In its con-
ducting form, gA forms a 24 A˚ long pore containing atrostatic potential energy. NH dipoles solvate permeat-
ing anions in the ClC chloride channel (Dutzler et al., single-file chain of up to eight water molecules. Contrary
to that found in AQPs, this water chain mediates the2002). In AQPs, only the NH2 groups of Arg206 and of
the two NPA Asn side chains can stabilize anions. Both relay of protons very effectively (Akeson and Deamer,
1991; Schumaker et al., 2000). The dipole moments ofthe scarcity and the alternance of CO and NH2 groups
pointing into the lumen discourage the coordination of four indole rings located at each of the two lipid-water
interfaces give rise to a moderate electrostatic fieldanions and cations alike.
which has been proposed to increase the proton affinity
of water in the lumen relative to bulk (Gowen et al., 2002).Selectivity Filter
It was proposed that interruption of the hydrogen- However, this field is not strong enough to overcome the
intrinsic preference of the unprotonated water chain forbonded chain prevents the passage of the excess pro-
ton through the selectivity filter of AQPs (de Groot and a fully (or nearly fully) polarized state (Pome`s and Roux,
1998, 2002). In the absence of H, membrane voltage,Grubmu¨ller, 2001). The water distribution obtained by
Tajkhorshid et al. (2002), however, suggests that this and strong electrostatic gradients in the gA channel,
there is no preference for one polarized conformationdefect is transient in GlpF, appearing in turn between
water molecules 1 and 2, and 2 and 3 (see Figure 2). over the other. Contrary to GlpF, in gA the structure
and fluctuations of the hydrogen-bonded network areOur results indicate that chain interruption does not have
a significant effect in the presence of an excess proton. governed by the properties of the water wire and modu-
lated by short-range interactions with the channel, bothFurthermore, the conserved Arg in the selectivity filter
constitutes only a shoulder in the electrostatic energy hydrogen bonding, charge-dipole, and dipole-dipole (Yu
and Pome`s, 2003). Each of the single-file water mole-profile at the periplasmic vestibule region of GlpF, and
does not oppose proton hopping in the single-file region. cules is a good proton host thanks to the proper coordi-
nation by carbonyl groups lining the pore of the gAThe lack of a barrier for proton movement past Arg206
(Figure 4) is somewhat surprising considering the strong channel (Pome`s and Roux, 2002). Accordingly, the PMF
for the translocation of an ionic defect in the pore of gAcoulombic repulsion between two positive charges, but
can be related to the contribution of other charged is nearly flat, with a slight (1 kcal/mol) preference for
the middle of the channel (Pome`s and Roux, 2002). Ingroups of the protein, such as Glu43, Glu152, and
Asp207, to the electrostatic field of the protein in this addition, hydrogen-bond donation to carbonyl groups
assists the translocation of a bonding defect by stabiliz-region (Jensen et al., 2003). In fact, neither the transfer
of an excess proton nor the preferential polarization of ing partially oriented conformations of the unprotonated
water chain. Thus, gA solvates both ionic and bondingwater molecules in the selectivity filter are determined by
the charge on Arg206. By itself, the guanidinium group of defects (Pome`s and Roux, 2002). The sensitivity of the
relay mechanism to the charge distribution of the chan-Arg206 would tend to impose a bimodal orientation to
nearby water molecules; in particular, water molecules nel, which is evidenced in the present study, was uncov-
ered in a recent study of dioxolane-linked gA, where the3, 4,… would point in the opposite direction to that
observed in this and other studies (de Groot and Grub- movement of proton was shown to be modulated by
conformational fluctuations exchanging locally dis-mu¨ller, 2001; Tajkhorshid et al., 2002; Jensen et al.,
2003). This shows that the charge-charge repulsion be- torted states of the channel on a (ns) time-scale com-
mensurate with proton permeation (Yu and Pome`s,tween a cation and Arg206 is more than compensated
by the electrostatic gradient culminating in the NPA re- 2003). In these distorted states, the dipole moments of
two carbonyl groups act as conformational switchesgion, in accord with the analysis of Jensen et al. (2003).
Further studies will help clarify the role of specific com- modifying the preferred organization of the unproton-
ated chain and confining the excess proton throughponents of the channel in the control of structural diffu-
sion and ion blockage in general. dipole-dipole and charge-dipole interactions.
Mechanism of Proton Blockage
ConclusionsBecause our PMF calculation for the excess proton is
By using a combination of methods comprising umbrellaconfined to the single-file region of the pore, we cannot
sampling simulations with a polarizable/dissociable modelconclude on the precise magnitude of the overall free
and macroscopic continuum electrostatic calculations,energy barrier to proton transfer. Nevertheless, it should
we were able to quantify the energetics of proton trans-be noted that the approach of an excess proton into
port through GlpF. The detailed picture of structuralthe narrow pore is opposed by the cost of polarizing
diffusion and the reversible work for the translocationthe unprotonated water chain, which, as indicated by the
of ionic and bonding defects provided by the presentturn PMF, is about 14 and 7 kcal/mol for the approach of
study clarify the physical basis of proton blockage byH from the periplasm and from the cytoplasm, respec-
AQPs. Tajkhorshid et al. (2002) noted the bipolar organi-tively (Figure 3). Thus, the free energy barriers opposing
zation of the water chain and remarked that the orienta-hop and turn steps each contribute to the molecular
tion of water molecules is incompatible with the penetra-mechanism of exclusion of protons from the channel.
tion of H into the pore from either side. The present
study confirms the stringent control of water orientationBlockage versus Conduction
but also indicates that, although it participates in theThe above mechanism contrasts with that obtained for
the rapid, passive transport of H in the gramicidin cat- exclusion of H from the pore, the bipolar organization

Structure
72
qH are the formal charges of single-file water O and H atoms. Forof the chain does not in itself preclude proton hopping
the TIP3P water model, qO and qH are
0.834 and 0.417 e, whereasin the single-file region. Rather, both turn and hop steps
for the PM6 model, qO and qH are
2 and 1 e, respectively. Thisof structural diffusion are opposed by electrostatic gra-
reaction coordinate was used in previous studies of the turn step
dients around the NPA site. Together, these two effects (Pome`s and Roux, 1998, 2002; Yu and Pome`s, 2003). The PMF
conspire to block proton permeation through the chan- profile was calculated from the equilibrium probability distributions
of z obtained from the simulations. A biasing potential energy func-nel. In view of the highly conserved features of the nar-
tion was imposed onz to force the reorientation of the unprotonatedrow pore region (Unger, 2000), it is likely that the electro-
water chain. Harmonic biasing functions of the form Vi(z)  (1/2) kistatic origin of proton blockage extends to other AQPs.
(z
zi)2 were used to carry out this umbrella sampling calculationElectrostatic forces were also shown to play an impor-
(Torrie and Valleau, 1974, 1977). For all the windows the harmonic
tant role in the charge selectivity of the KcsA potassium ki coefficients were set to 20 kcal·mol
1·e
2·A˚
2. For the TIP3P water
channel (Roux and MacKinnon, 1999) and to modulate chain, the reference zi varied in increments of 0.25 e·A˚, from
5.0
to 5.0 e·A˚. Each window of the umbrella simulation consisted of 10the Grotthuss relay mechanism in chemically modified
ps of equilibration and 40 ps of production run for data collection.variants of the gramicidin channel (Yu and Pome`s, 2003).
In the case of the PM6 water chain, zi varied in increments of 0.25Although further studies will help clarify the structural
e·A˚ from
11.0 to 11.0 e·A˚, and the same protocol was used for eachorigin of permeant selectivity, this work provides mean-
of the 89 windows. The PMF profile was computed by unbiasing and
ingful insight into the physical basis for the control of combining windows with the weighted histogram analysis method
H transport in membrane proteins, one of the most (WHAM) (Kumar et al., 1992; Roux, 1995).
Proton Hoppingimportant reactions in biology.
We define a new reaction coordinate to follow the displacement of
an excess proton along the water chain as follows:   (OiWOi)
1(OiExperimental Procedures
zOiWOi), where WOi  Hj fsw(rOiHj)
2 is a weighting function, fsw(rOiHj) 
1/(1  exp[{rOiHj
rsw}/dsw]) is a switching function, and Oi, Hj areMolecular Model
respectively oxygen and hydrogen atoms of the single-file waterThe initial conformation of GlpF was taken from a previous simula-
chain. Based on earlier studies using the PM6 model (Pome`s andtion study of the homotetrameric channel in a hydrated lipid mem-
Roux, 1996), we chose dsw  0.05 A˚ and rsw  1.4 A˚. This continuousbrane (Tajkhorshid et al., 2002). One monomer, together with slabs
and derivable reaction coordinate keeps track of the water mole-of water molecules nearby, was retained, and the rest of the system
cule(s) that are close to three hydrogen nuclei, yielding the locationwas discarded (Figure 1). The monomeric system consists of 3839
of the oxygen atom for a hydronium ion OH3 and that of the sharedprotein atoms, 9 water molecules in the pore, and 1383 water mole-
proton in a symmetric Zundel ion (O2H5). The PMF profile of transfercules in the bulk, for a total of 8315 atoms. The inner core, an
of an excess proton in the water chain was computed from equilib-orthorhombic region of 30
14
14 A˚3 and consisting of
1100
rium distributions of  obtained from umbrella sampling simulations.atoms, was allowed to move during the MD simulations; the rest of
A biasing potential energy function was imposed on  to force thethe system was kept fixed.
sampling over the entire pore region. 43 successive sampling win-The CHARMM force field, version 22 (Brooks et al., 1983; MacKer-
dows with harmonic biasing functions of the form Vi()  (1/2) ki (
ell et al., 1998), was used to model the protein. Bulk water was
i)2 were used, with ki  0.75 kcal·mol
1A˚
2 and i ranging from
6.5represented by the TIP3P force field (Jorgensen et al., 1983) and
to 14.5 A˚ in increments of 0.5 A˚. Each window of umbrella samplingthe nine water molecules in the channel pore were modeled succes-
consisted of 20 ps of equilibration followed by 60 ps of production.sively with TIP3P and with the PM6 model (Stillinger and David,
The reaction coordinate was recorded every 2 fs. The total simula-1978; Stillinger, 1979; Weber and Stillinger, 1982; Pome`s and Roux,
tion time required to build the PMF profile for H translocation was1996, 1998). Both force fields were used in the study of the turn
3.44 ns.step of the Grotthuss mechanism, whereas the PM6 force field was
employed in simulations with an excess proton. PM6 is a polarizable
and dissociable model of water that consists of O2
and H moieties. Electrostatic Calculations
This empirical model has been shown to capture the essential fea- We computed the total electrostatic potential (ESP) for the passage
tures of the mechanism of H transport in biological proton ducts of a cation through the GlpFpore by solving the linearized Poisson-
(Mei et al., 1998; Pome`s and Roux, 2002). Boltzmann (PB) equation. The PBEQ module of CHARMM version
Geometric constraints were applied on the single-file chain of 28 was used together with the set of optimized atomic radii for
nine water molecules and on the bulk water to prevent the entrance amino acids (Nina et al., 1997). All explicit water molecules were
or exit of water into or from the single-file region of the pore. Planar removed from the molecular system. The ESP was computed with
restraints were imposed on water molecules at the interface of single a cation placed at successive positions along the channel pore. In
file and bulk. To check the consistency of this finite-size model, we the continuum approximation, the total electrostatic free energy of
computed the distribution of oxygen and hydrogen atoms (O, H) the cation in the pore can be expressed as the sum of contributions
of PM6 water molecules from a 2 ns simulation and compared it from the static field (which reflects the charge distribution of the
with earlier results obtained from simulations of the full tetrameric protein) and from the reaction field (which is governed by the low-
channel immersed in an explicit lipid bilayer with periodic boundary dielectric boundaries). The radius of the cation was chosen as that
conditions (Tajkhorshid et al., 2002). The finite-size model offers the of sodium (Nina et al., 1997) because the hydration free energy of
advantage of computational efficiency as well as allowing the study that ion,
98 kcal/mol (Burgess, 1978), is comparable to that of
of structural diffusion in the absence of water permeation. hydronium,
102 (Pearson, 1986) to
104 kcal/mol (Chambers et
al., 1996). The linearized PB equation was solved on a cubic grid
of 151 cells with a cell width of 1.0 A˚, with subsequent focusingMolecular Dynamics Simulations
The program CHARMM (Brooks et al., 1983) was employed to gener- using a 0.5 A˚ cell width. The low-dielectric region of the membrane
was represented by a slab perpendicular to the z axis and centeredate MD trajectories. The Langevin equations of motion were propa-
gated at 300 K with an integration step of 1 fs and a friction coeffi- at z  5 A˚. The thickness of the membrane slab was chosen as 40 A˚
based on a previous simulation study (Jensen et al., 2001). Thecient of 5 ps
1 applied to all heavy atoms. Nonbonded interactions
were calculated without any cut off. dielectric constants used in these calculations were respectively 2,
2, and 80 for protein, membrane, and water. A high-dielectric con-Water Reorientation
The PMF profiles were computed as follows. The configuration of stant (  80) was assigned to the interior of the pore. No trans-
membrane voltage was applied. The calculation was repeated withthe molecular system was recorded every 5 fs. The projection of
the total dipole moment of the lumen content (water chain) on the a membrane dielectric constant of 80 to gauge the sensitivity of the
results to the description of the membrane. Finally, the calculationz axis (channel axis) was calculated for each configuration at time
t as z(t )  qO Oi ZOi(t )  qH Hj ZHj (t ) in units of e·A˚, where qO and was repeated with smaller (1.4 A˚) and larger (2.03 A˚) cation radii to

Proton Blockage in Aquaporins
73
gauge the sensitivity of the results to the size of the ionic probe; and Klein, M.L. (1983). Comparison of simple potential functions for
simulating liquid water. J. Chem. Phys. 79, 926–935.the latter radius corresponds to a K ion (Nina et al., 1997).
Kozono, D., Yasui, M., King, L.S., and Agre, P. (2002). Aquaporin
Acknowledgments water channels: atomic structure and molecular dynamics meet clin-
ical medicine. J. Clin. Invest. 109, 1395–1399.
We thank Klaus Schulten and Wonpil Im for useful discussions. We Kumar, S., Bouzida, D., Swendsen, R.H., Kollman, P.A., and Rosen-
gratefully acknowledge the Canadian Institutes of Health Research berg, J.M. (1992). The weighted histogram analysis method for free-
(operating grant MOP43949) and the Ontario Centre for Genomics energy calculations in biomolecules, 1. The Method. J. Comput.
Computing for support. R.P. is a CRCP Chairholder. Chem. 13, 1011–1021.
MacKerell, A.D., Jr., Bashford, D., Bellott, M., Dunbrack, R.L., Jr.,
Received: August 23, 2003 Evanseck, J.D., Field, M.J., Fischer, S., Gao, J., Guo, H., Ha, S., et
Revised: September 29, 2003 al. (1998). All-atom empirical potential for molecular modeling and
Accepted: October 5, 2003 dynamics studies of proteins. J. Phys. Chem. B 102, 3586–3616.
Published: January 13, 2004
Mei, H.S., Tuckerman, M.E., Sagnella, D.E., and Klein, M.L. (1998).
Quantum nuclear ab initio molecular dynamics study of water wires.
References J. Phys. Chem. B 102, 10446–10458.
Mitchell, P. (1961). Coupling of phosphorylation to electron andAgmon, N. (1995). The Grotthuss mechanism. Chem. Phys. Lett.
hydrogen transfer by a chemi-osmotic type of mechanism. Nature224, 456–462.
191, 144–148.
Agre, P., Bonhivers, M., and Borgnia, M.J. (1998). The aquaporins,
Morais-Cabral, J.H., Zhou, Y., and Mackinnon, R. (2001). Energeticblueprints for cellular plumbing systems. J. Biol. Chem. 273, 4659–
optimization of ion conduction rate by the K selectivity filter. Nature4662.
414, 37–42.
Akeson, M., and Deamer, D.W. (1991). Proton conductance by the
Murata, K., Mitsuoka, K., Hirai, T., Walz, T., Agre, P., Heymann, J.B.,gramicidin water wire. Biophys. J. 60, 101–109.
Engel, A., and Fujiyoshi, Y. (2000). Structural determinants of water
Borgnia, M.J., Nielsen, S., Engel, A., and Agre, P. (1999). Cellular
permeation through aquaporin-1. Nature 407, 599–605.
and molecular biology of the aquaporin water channels. Annu. Rev.
Nagle, J.F., and Morowitz, H.J. (1978). Molecular mechanisms forBiochem. 68, 425–458.
proton transport in membranes. Proc. Natl. Acad. Sci. USA 75,
Brooks, B.R., Bruccoleri, R.E., Olafson, B.D., States, D.J., Swamina-
298–302.
than, S., and Karplus, M. (1983). CHARMM: a program for macro-
Nina, M., Beglov, D., and Roux, B. (1997). Atomic Born radii formolecular energy minimization and dynamics calculations. J. Com-
continuum electrostatic calculations based on molecular dynamicsput. Chem. 4, 187–217.
free energy simulations. J. Phys. Chem. B 101, 5239–5248.Burgess, J. (1978). Metal Ions in Solution (New York: J. Wiley).
Nollert, P., Harries, W.E.C., Fu, D., Miercke, L.J., and Stroud, R.M.Chambers, C.C., Hawkins, G.D., Cramer, C.J., and Truhlar, D.G.
(2001). Atomic structure of a glycerol channel and implications for(1996). Model for aqueous solvation based on class IV atomic
substrate permeation in aqua (glycero) porins. FEBS Lett. 504,charges and first solvation shell effects. J. Phys. Chem. 100, 16385–
112–117.16398.
Pearson, R.G. (1986). Ionization potentials and electron affinities inDenker, B.M., Smith, B.L., Kuhadja, F.P., and Agre, P. (1988). Identifi-
aqueous solution. J. Am. Chem. Soc. 108, 6109–6114.cation, purification, and partial characterization of a novel Mr 28,000
Pome`s, R., and Roux, B. (1996). Structure and dynamics of a protonintegral membrane protein from erythorcytes and renal tubules. J.
wire: a theoretical study of H translocation along the single-fileBiol. Chem. 263, 15634–15642.
water chain in the gramicidin channel. Biophys. J. 71, 19–39.de Groot, B.L., and Grubmu¨ller, H. (2001). Water permeation across
Pome`s, R., and Roux, B. (1998). Free energy profile for H conduc-biological membranes: mechanism and dynamics of aquaporin-1
tion along hydrogen-bonded chains of water molecules. Biophys.and GlpF. Science 294, 2353–2357.
J. 75, 33–40.de Grotthuss, C.J.T. (1806). Me´moire sur la de´compostion de l’eau
Pome`s, R., and Roux, B. (2002). Molecular mechanism of H conduc-et des corps qu’elle tient en dissolution a` l’aide de l’e´lectricite´ galva-
tion in the single-file water chain of the gramicidin channel. Biophys.nique. Ann. Chim. 58, 54–74.
J. 82, 2304–2316.Dutzler, R., Campbell, E.B., Cadene, M., Chait, B.T., and Mackinnon,
Preston, G.M., and Agre, P. (1991). Isolation of the cDNA for erythro-R. (2002). X-ray structure of a CIC chloride channel at 3.0 angstrom
cyte integral membrane protein of 28 kilodaltons: member of anreveals the molecular basis of anion selectivity. Nature 415, 287–294.
ancient channel family. Proc. Natl. Acad. Sci. USA 88, 11110–11114.Fu, D., Libson, A., Miercke, L.J., Weitzman, C., Nollert, P., Krucinski,
Preston, G.M., Carroll, T.P., Guggino, W.B., and Agre, P. (1992).J., and Stroud, R.M. (2000). The structure of a glycerol-conducting
Appearance of water channels in Xenopus oocytes expressing redchannel reveals the basis for its selectivity. Science 290, 481–486.
cell CHIP28 protein. Science 256, 385–387.Gowen, J.A., Markham, J.C., Morrison, S.E., Cross, T.A., Busath,
Ren, G., Reddy, V.S., Cheng, A., Melnyk, P., and Mitra, A.K. (2001).D.D., Mapes, E.J., and Schumaker, M.F. (2002). The role of Trp
Visualization of a water-selective pore by electron crystallographyside chains in tuning single proton conduction through gramicidin
in vitreous ice. Proc. Natl. Acad. Sci. USA 98, 1398–1403.channels. Biophys. J. 83, 880–898.
Roux, B. (1995). The calculation of the potential of mean force usingHumphrey, W., Dalke, A., and Schulten, K. (1996). VMD: visual molec-
computer simulations. Comp. Phys. Commun. 91, 275–282.ular dynamics. J. Mol. Graph. 14, 33–38.
Roux, B., and Mackinnon, R. (1999). The cavity and pore helices inJensen, M.Ø., Tajkhorshid, E., and Schulten, K. (2001). The mecha-
KcsA K channel: electrostatic stabilization of monovalent cations.nism of glycerol conduction in Aquaglyceroporins. Structure 9,
Science 285, 100–102.1083–1093.
Roux, B., Berne`che, S., and Im, W. (2000). Ion channels, permeation,Jensen, M.Ø., Park, S., Tajkhorshid, E., and Schulten, K. (2002).
and electrostatics: insight into the function of KcsA. BiochemistryEnergetics of glycerol conduction through aquaglyceroporin GlpF.
39, 13295–13306.Proc. Natl. Acad. Sci. USA 99, 6731–6736.
Schumaker, M.F., Pome`s, R., and Roux, B. (2000). A combined mo-Jensen, M.Ø., Tajkhorshid, E., and Schulten, K. (2003). Electrostatic
lecular dynamics diffusion model of single proton conductiontuning of permeation and selectivity in aquaporin water channels.
through gramicidin. Biophys. J. 79, 2840–2857.Biophys. J., 85, 2884–2899.
Jorgensen, W.L., Chandrasekhar, J., Madura, J.D., Impey, R.W., Stillinger, F.H. (1979). Dynamics and ensemble averages for the

Structure
74
polarization models of molecular interactions. J. Chem. Phys. 71,
1647–1651.
Stillinger, F.H., and David, C.W. (1978). Polarization model for water
and its ionic dissociation products. J. Chem. Phys. 69, 1473–1484.
Sui, H., Han, B.G., Lee, J.K., Walian, P., and Jap, B.K. (2001). Struc-
tural basis of water-specific transport through the AQP1 water chan-
nel. Nature 414, 872–878.
Tajkhorshid, E., Nollert, P., Jensen, M., Miercke, L.J.W., O’Connell,
J., Stroud, R.M., and Schulten, K. (2002). Control of the selectivity
of the aquaporin water channel family by global orientational tuning.
Science 296, 525–530.
Torrie, G.M., and Valleau, J.P. (1974). Monte-Carlo free energy esti-
mates using non-Boltzmann sampling—application to subcritical
Lennard-Jones fluid. Chem. Phys. Lett. 28, 578–581.
Torrie, G.M., and Valleau, J.P. (1977). Non-physical sampling distri-
butions in Monte-Carlo free energy estimation: umbrella sampling.
J. Comp. Physiol. 23, 187–199.
Tuckerman, M., Laasonen, K., Sprik, M., and Parrinello, M. (1995).
Ab initio molecular dynamics simulation of the solvation and trans-
port of H3O and OH
ions in water. J. Phys. Chem. 99, 5749–5752.
Unger, V.M. (2000). Fraternal twins: AQP and GlpF. Nat. Struct. Biol.
7, 1082–1084.
Weber, T.A., and Stillinger, F.H. (1982). Reactive collisions of H3O
and OH
studied with the polarization model. J. Phys. Chem. 8,
1314–1318.
Yu, C.-H., and Pome`s, R. (2003). Functional dynamics of ion chan-
nels: modulation of proton movement by conformational switches.
J. Am. Chem. Soc., 125, 13890–13894.
Zhu, F., Tajkhorshid, E., and Schulten, K. (2001). Molecular dynamics
study of aquaporin-1 water channel in a lipid bilayer. FEBS Lett.
504, 212–218.
Zhu, F., Tajkhorshid, E., and Schulten, K. (2002). Pressure-induced
water transport in membrane channels studied by molecular dynam-
ics. Biophys. J. 83, 154–160.