help button home button Biophys. J.
HOME HELP FEEDBACK SUBSCRIPTIONS ARCHIVE SEARCH TABLE OF CONTENTS
This Article
Right arrow Abstract Freely available
Right arrow Full Text (PDF)
Right arrow Alert me when this article is cited
Right arrow Alert me if a correction is posted
Services
Right arrow Similar articles in this journal
Right arrow Similar articles in PubMed
Right arrow Alert me to new issues of the journal
Right arrow Download to citation manager
Right arrow reprints & permissions
Citing Articles
Right arrow Citing Articles via HighWire
Right arrow Citing Articles via Google Scholar
Google Scholar
Right arrow Articles by Jensen, M. O.
Right arrow Articles by Schulten, K.
Right arrow Search for Related Content
PubMed
Right arrow PubMed Citation
Right arrow Articles by Jensen, M. O.
Right arrow Articles by Schulten, K.
Biophysical Journal 85:2884-2899 (2003)
© 2003 The Biophysical Society

Electrostatic Tuning of Permeation and Selectivity in Aquaporin Water Channels

Morten Ø. Jensen, Emad Tajkhorshid and Klaus Schulten

Theoretical and Computational Biophysics Group, Beckman Institute, University of Illinois at Urbana-Champaign, Illinois

Correspondence: Address reprint requests to K. Schulten, E-mail: kschulte{at}ks.uiuc.edu.


   ABSTRACT
TOP
 ABSTRACT
INTRODUCTION
 METHODS
 RESULTS AND DISCUSSION
 CONCLUSION
 ACKNOWLEDGEMENTS
 REFERENCES
 
Water permeation and electrostatic interactions between water and channel are investigated in the Escherichia coli glycerol uptake facilitator GlpF, a member of the aquaporin water channel family, by molecular dynamics simulations. A tetrameric model of the channel embedded in a 16:0/18:1c9-palmitoyloleylphosphatidylethanolamine membrane was used for the simulations. During the simulations, water molecules pass through the channel in single file. The movement of the single file water molecules through the channel is concerted, and we show that it can be described by a continuous-time random-walk model. The integrity of the single file remains intact during the permeation, indicating that a disrupted water chain is unlikely to be the mechanism of proton exclusion in aquaporins. Specific hydrogen bonds between permeating water and protein at the channel center (at two conserved Asp-Pro-Ala "NPA" motifs), together with the protein electrostatic fields enforce a bipolar water configuration inside the channel with dipole inversion at the NPA motifs. At the NPA motifs water-protein electrostatic interactions facilitate this inversion. Furthermore, water-water electrostatic interactions are in all regions inside the channel stronger than water-protein interactions, except near a conserved, positively charged Arg residue. We find that variations of the protein electrostatic field through the channel, owing to preserved structural features, completely explain the bipolar orientation of water. This orientation persists despite water translocation in single file and blocks proton transport. Furthermore, we find that for permeation of a cation, ion-protein electrostatic interactions are more unfavorable at the conserved NPA motifs than at the conserved Arg, suggesting that the major barrier against proton transport in aquaporins is faced at the NPA motifs.


   INTRODUCTION
TOP
 ABSTRACT
 INTRODUCTION
METHODS
 RESULTS AND DISCUSSION
 CONCLUSION
 ACKNOWLEDGEMENTS
 REFERENCES
 
All life-forms need to continuously and precisely regulate water content of their cells and organs through exchange of water with their surroundings. Although all cell membranes are to some extent water-permeable, water transport through lipid bilayers is slow and inadequate for exchange of large volumes of water, e.g., as occurring in kidneys, or for a rapid response of cells to osmotic stress, as in red blood cells. The high water permeability of certain cell membranes was explained by the discovery of aquaporins (AQPs), integral membrane proteins that are specialized for selective transport of water (Preston et al., 1992Go; Agre et al., 1998Go; Borgnia et al., 1999Go; Yasui et al., 1999Go).

AQPs are present in all life-forms, and more than 100 have been characterized (Borgnia et al., 1999Go). Eleven human AQPs have been identified, and their impaired function is implicated in pathological situations, such as nephrogenic diabetes insipidus and congenital cataract (Borgnia et al., 1999Go; Kozono et al., 2002Go). All members of the AQP family provide a passive means of water transport. A subfamily of AQPs, aquaglyceroporins, possesses the extended capability of stereo-selectively conducting small, linear carbohydrates in addition to water (Heller et al., 1980Go; Fu et al., 2000Go; Borgnia et al., 1999Go; Jensen et al., 2001Go, 2002Go). The most extensively studied aquaglyceroporin is the Escherichia coli glycerol uptake facilitator GlpF, which is used by the organism to accelerate absorption of glycerol from the environment, especially at low glycerol concentrations (Richey and Lin, 1972Go).

AQPs form tetramers in the membrane; the biological significance of oligomerization was not clear until recent reports implicated that the central pore allows ion transport (Yool and Weinstein, 2002Go; Hazama et al., 2002Go). Atomic resolution models of human aquaporin-1 (AQP1) based on electron microscopy of two-dimensional crystals provided the first insight into the tetrameric architecture and folding of AQPs (Murata et al., 2000Go; Ren et al., 2001Go). The first molecular dynamics (MD) simulations performed on these models, however, indicated that the models needed further improvement (Zhu et al., 2001Go; de Groot et al., 2001Go). GlpF was the first AQP whose high-resolution structure (2.2 Å) was determined by x-ray crystallography (Fu et al., 2000Go). The first GlpF structure included three glycerol molecules, a natural substrate for GlpF, in each monomer. MD simulations of this glycerol-saturated GlpF (GlpF+G) resulted in a detailed description of protein-substrate interactions and conduction pathway, which for the first time explained why two half-membrane spanning (re-entrant) loops with unusual secondary structures are conserved in the architecture of AQPs (Jensen et al., 2001Go), as sketched for a single monomer in Fig. 1. In each monomer, the transmembrane segment of the channel is formed by six helices and the re-entrant loops that meet each other at the center of the channel. Conservation of the re-entrant loops, which are related by quasi-twofold symmetry, is of structural and functional importance (Murata et al., 2000Go; Ren et al., 2001Go; Fu et al., 2000Go; Zhu et al., 2001Go; Jensen et al., 2001Go; Sui et al., 2001Go). Approximately one-half of each loop is nonhelical and defines a curve-linear conduction pathway, constituted by backbone carbonyl groups that are exposed toward the channel interior forming here hydrogen bonds with the permeant substrate (Jensen et al., 2001Go). The otherwise energetically unfavorable secondary structure of the nonhelical parts is stabilized through multiple hydrogen bonds of backbone N-H groups to two highly conserved glutamate residues (Fu et al., 2000Go; Jensen et al., 2001Go; Nollert et al., 2001Go).



View larger version (93K):
[in this window]
[in a new window]
  		FIGURE 1  Schematic representation of a single GlpF monomer. The six membrane-spanning helices are drawn as green cylinders. The two re-entrant, half-membrane spanning repeats meeting in the center of the channel at the position of the NPA motifs are drawn in red with the M3 and M7 {alpha}-helices represented as cylinders and the conduction pathway represented as red threads.

 
The other half of each loop, known as M3 and M7 in GlpF, and HB and HE in AQP1 (Murata et al., 2000Go; Ren et al., 2001Go; Sui et al., 2001Go), is {alpha}-helical. These helices terminate at the channel center, where two conserved (Asn-Pro-Ala) NPA motifs (Fu et al., 2000Go; Murata et al., 2000Go; Ren et al., 2001Go; Sui et al., 2001Go) meet each other. Multiple hydrogen bonds between the NPA motifs (Zhu et al., 2001Go; Jensen et al., 2001Go) result in a configuration in which one of the amide hydrogen atoms of each asparagine points toward the pore. Together they provide two very close hydrogen-bonding sites for the permeating substrate at the center of the channel.

These structural details of GlpF (Fu et al., 2000Go) were used to refine AQP1 models (de Groot et al., 2001Go; Zhu et al., 2001Go), before a 2.7 Å resolution structure of AQP1 was determined by x-ray crystallography (Sui et al., 2001Go). In addition to a better description of the channel interior, the high-resolution structure of AQP1 revealed the position of four water molecules inside the channel. Three-dimensional structures of GlpF have also been solved under low glycerol concentration conditions, resulting in water-filled, glycerol-free GlpF (GlpF-G) and revealing the presence of a single file of water molecules inside the channel (Tajkhorshid et al., 2002Go). Comparison of the GlpF+G and GlpF-G structures shows a remarkable similarity between the conformations of side chains in the two channels, indicating that the channel can accommodate glycerol with only minor conformational changes occurring mainly at its narrowest part, denoted as the selectivity filter (SF) (Fu et al., 2000Go; Jensen et al., 2001Go; Tajkhorshid et al., 2002Go; Grayson et al., 2003Go).

Substrate selectivity is crucial in the function of AQPs; these proteins should provide a facilitatory mechanism for transmembrane water flow without affecting the electro-chemical properties of the membrane. Water pores in AQPs are impermeable to charged species. The narrow size of the channel combined with its hydrophobic interior seems to be an effective means of preventing ions from entering the channel; hydrated ions are too large to enter the pore, and, besides an unfavorable increase in entropy, the solvation energy provided by the channel interior is too weak to dehydrate ions. Selectivity against protons, essential in conservation of the transmembrane proton gradient and in the bioenergetics of the cell, may, however, not be explained as easily, since protons can, be transported very fast through a chain of water molecules according to the Grotthuss mechanism (de Grotthuss, 1806Go; Schulten and Schulten, 1986Go; Agmon, 1995Go; Pomès and Roux, 1996Go, 2002Go; Brewer et al., 1999Go; Smondyrev and Voth, 2002Go).

Proton exclusion is indeed effectively accomplished by AQPs through tuning of the orientation of water molecules inside the channel (Tajkhorshid et al., 2002Go) that prevents them from adopting a uniform arrangement required for a proton wire, e.g., as observed in gramicidin A (Pomès and Roux, 1996Go, 2002Go). This mechanism of proton exclusion has been investigated most recently for water conducting nanotubes (Zhu and Schulten, 2003Go). An alternative mechanism of proton exclusion, based on the disruption of the single file water in the SF region, i.e., the narrowest part of the channel, or in the NPA region was suggested by others (de Groot et al., 2001Go; Ren et al., 2001Go).

Since substrates of AQPs, e.g., water and linear sugar molecules, are neutral, electrostatic forces might not seem to be critical for the function of the channel. It turns out, however, that such forces, which have proven crucial in the function of ion channels, e.g., the KcsA potassium channel (Roux and MacKinnon, 1999Go; Bernèche and Roux, 2001Go), also play important roles in the selectivity mechanisms employed by AQPs (Jensen et al., 2001Go; Tajkhorshid et al., 2002Go; Grayson et al., 2003Go). We report here a detailed analysis of water dynamics and electrostatic interaction of permeating water molecules with the GlpF channel using trajectories from multinanosecond MD simulations. MD has been applied to study water and glycerol permeation in AQPs (Zhu et al., 2001Go, 2002Go; Jensen et al., 2001Go, 2002Go; de Groot et al., 2001Go; Tajkhorshid et al., 2002Go), as well as water transport in nanotubes (Hummer et al., 2001Go; Zhu and Schulten, 2003Go). To investigate in further detail how water conduction and proton exclusion are established by conserved features of the protein architecture of AQPs, we analyze here the dynamics of water inside the channel, and investigate and decompose the electrostatic forces acting on water in various regions of the channel during its permeation.


   METHODS
TOP
 ABSTRACT
 INTRODUCTION
 METHODS
RESULTS AND DISCUSSION
 CONCLUSION
 ACKNOWLEDGEMENTS
 REFERENCES
 
In this section we describe the modeling of tetrameric GlpF in a hydrated lipid membrane, the details of simulation protocols, and the procedures applied for trajectory analysis.

Modeling
The crystal structure of a GlpF monomer in its glycerol-bound form (GlpF+G) (Fu et al., 2000Go) was obtained from the Research Collaboratory for Structural Bioinformatics Data Bank (PDB) (Bernstein et al., 1977Go), PDB entry code 1FX8 (Fu et al., 2000Go). Missing atoms of Arg257 and hydrogen atoms were added using Insight II (Insight, 1998Go) and XPLOR (Brünger, 1992Go), respectively. The tetrameric structure of GlpF was generated with VMD (Humphrey et al., 1996Go) using the coordinate transformation matrices provided in the PDB file.

A 16:0/18:1c9-palmitoyloleylphosphatidylethanolamine (POPE) lipid bilayer was used to mimic the E. coli membrane. Initially 480 POPE lipids, all in cis-conformation were placed on a regular hexagonal lattice with a lattice vector of 8.1 Å. The membrane plane was in the xy plane, and the membrane normal was along the z-direction. The GlpF tetramer including channel-bound water and glycerol molecules resolved crystallographically was inserted into the membrane using the location of tyrosine residues and the proposed mid-membrane residues (Thr18, Gly49, Asn68, Pro69, Ala70, Gly96, Thr156, Gly184, and Gly243) (Fu et al., 2000Go) to adjust the normal position of the membrane relative to the protein. An optimal position, where equal distances from the lipid headgroups to the above residues resulted, was obtained by separately adjusting the two monolayers relative to the protein. As a result the hydrophobic surface of the protein was embedded by the hydrophobic core of the membrane, with a vertical phosphorus-phosphorus separation between the two monolayers initially measuring 39.6 Å. Upon visual inspection in VMD (Humphrey et al., 1996Go), 83 lipids overlapping with the protein were removed.

The program SOLVATE (Grubmüller, 1996Go) was used to add a 20 Å ellipsoidal solvation shell around the membrane and the solvent-exposed parts of the protein. From this system a rectangular box, based on the positions of the nitrogen atoms of the ethanolamine headgroups of the lipids, was cut out. In the resulting box, all water molecules located inside the hydrophobic core of the membrane were discarded. The normal dimension of the rectangular system ensured ~10 Å hydration of each POPE monolayer. Four water molecules were replaced by chloride ions, to neutralize the net charge of the system.

For the study of water conduction, we modified the GlpF+G/POPE/water system described above by replacing three glycerol molecules (per monomer), as present in GlpF+G (1FX8) (Fu et al., 2000Go), by water. Water insertion using DOWSER (Hermans et al., 1998Go) placed in the constriction region of each monomer eight water molecules in a single file in close accord with the later reported crystal structures of GlpF-G (PDB entry codes 1LDI and 1LDA) (Tajkhorshid et al., 2002Go) and with the resulting equilibrium distribution of water inside the channel as obtained from MD simulations (Tajkhorshid et al., 2002Go).

The resulting water-saturated tetramer (GlpF-G) was substituted into a hydrated POPE membrane, pre-equilibrated for 100 ps with fixed protein at constant temperature and pressure (310 K, 1 atm) (Jensen et al., 2001Go). The total GlpF-G/POPE/water system counted 106,245 atoms including the GlpF tetramer (15,356 atoms), 317 POPE lipid molecules (39,625 atoms), 17,068 water molecules (51,204 atoms), and four glycerol molecules (56 atoms; harmonic constraints with a force constant of 5 kcal/mol/Å2 were applied to all glycerol carbon atoms to ensure that glycerol molecules remained in the bulk water region) and four chloride ions in the bulk region, with initial dimensions of 122.0 x 112.7 x 77.0 Å3. The corresponding dimensions of the pure GlpF tetramer are 73.4 x 73.4 x 57.0 Å3. The system was subsequently minimized and subjected to MD for another 100 ps with fixed protein. The protein was then released, the full system was minimized, and the simulation was conducted for 5 ns in the NPT ensemble with the first nanosecond considered equilibration, and the last 4 ns used for analysis.

The program NAMD (Kalé et al., 1999Go) and the CHARMM27 parameter set (Schlenkrich et al., 1996Go; MacKerell Jr. et al., 1998Go; Feller and MacKerell, 2000Go) were used for all simulations. Full periodic boundary conditions were imposed. The particle mesh Ewald method (Darden et al., 1993Go) was used for computation of electrostatic forces without truncation. Langevin dynamics was used to maintain a constant pressure of 1 atm and a constant temperature of 310 K (Feller et al., 1995Go). The applied temperature represents the natural environment of E. coli, and is above the gel-liquid crystalline phase transition temperature of POPE. A W48F/F200T double mutant of GlpF with increased water permeability (Tajkhorshid et al., 2002Go) was prepared from the GlpF-G/POPE system (equilibrated for 1 ns) by mutating the appropriate residues. The mutant was simulated for 5 ns in the NPT ensemble using the same protocol as for the wild-type (WT), while again considering the first 1-ns equilibration (Tajkhorshid et al., 2002Go; the mutant structure is accessible with PDB entry code 1LDF).

To investigate the electrostatic influence of the protein on the structure and dynamics of water, in addition to the native WT system above denoted (I), we carried out simulations on four charge-modified systems denoted (II)–(V): in (II) the charges of the NH2 groups of Asn68 and Asn203 of the NPA motifs were turned off; in (III) charge contributions from the backbone atoms of the M3 (residues 70–79) and M7 (residues 205–217) helices were turned off; (II) and (III) were combined in (IV); and in (V) a "hydrophobic" channel was created by combining (IV) with turning off the charges of the carbonyl groups of residues 64–67 and 196–201, which constitute the conduction pathway (Jensen et al., 2001Go). Due to the (charge) group concept of the CHARMM27 parameter set (Schlenkrich et al., 1996Go; MacKerell Jr. et al., 1998Go; Feller and MacKerell, 2000Go), we also set the charges of the Gly209 side chain and of C{delta} and H{delta} of Pro210 to zero, thereby avoiding introduction of net charges. In systems (III)–(V), all atoms of the M3 and M7 helices with charges turned off were fixed, as were the carbonyl groups in (V). The initial configuration used in systems (II)–(V) was the one of system (I) already equilibrated for 1 ns. We simulated systems (II)–(V) for 0.5 ns under the same conditions as used for simulating system (I).

Analysis
All trajectories were analyzed with the molecular graphics and analysis program VMD (Humphrey et al., 1996Go). Numeric values presented for a protein monomer, or atom subsets thereof, as emerging from the analysis procedures described below, are all averaged over the four monomers of the tetramer. Furthermore, various properties are calculated for water molecules occupying different regions of the channel. In assigning the calculated values to different regions of the channel, the position of the oxygen atom of water is used. When a pair of water molecules is considered, e.g., in calculation of pair correlations, the midpoint of the two oxygen atoms is adopted.

Water permeation
Water diffusion inside the channel occurs in single file and is highly correlated. Such a concerted translocation of water in single file, as also encountered in carbon nanotubes (Hummer et al., 2001Go; Berezhkovskii and Hummer, 2002Go; Zhu and Schulten, 2003Go), can be described by a continuous-time random-walk model (Berezhkovskii and Hummer, 2002Go). In that model, the hopping rate k is the frequency of water file translocation by a distance equal to that of the average water-water spacing . is the mean hopping time, i.e., the average time between successive single file translocations (Berezhkovskii and Hummer, 2002Go). From 2Dt {approx} <{Delta}z2(t)> (for t -> {infty}) we obtain the one-dimensional diffusion coefficient D with <{Delta}z2(t)> being the center of mass mean square displacement along the channel axis z:

(1)
N is here the number of gliding time origins and nCR is the number of water molecules within the constriction region (CR {equiv} -6 Å < z < 14 Å) between t' and t' + t. Any water molecule diffusing outside the constriction region within a given time window starting at t' is discarded in the remaining part of the calculation within that time window. With being the average number of single-file water molecules in the constriction region we obtain the bi-directional conduction rate k±, i.e., the number of permeating water molecules per channel per time, as (Berezhkovskii and Hummer, 2002Go; Zhu and Schulten, 2003Go; Zhu et al., 2003Go). We note in passing that the uni-directional conduction rate equals k±/2).

Water-water correlation
To quantify the correlated motion of water within the constriction region of the channel, we, as others (de Groot et al., 2001Go) define, for each position z along the channel axis, a nearest-neighbor cross correlation coefficient c(z), as an average over all pairs of adjacent water molecules, i and j visiting z (meaning that the midpoint of the two correlated oxygen atoms is z):

(2)
where


and

Single file disruption
For each position z along the channel axis, a single file disruption ratio, dr(z), based on the distance between any two adjacent water molecules i and j visiting z, was calculated as

(3)
Again, and NDWP (z) and NWP (z) are the total number of disrupted water pairs and the total number of water pairs, respectively, counted at position z within the constriction region. If the distance between the oxygen atoms of the i,jth nearest-neighbor water pair exceeded 3.75 Å, the single file was considered disrupted at z, i.e., at the midpoint of that water pair.

Water orientation
The equilibrium water orientation inside the channel was characterized by the Legendre polynomials: P1(z) = <cos({theta})z> and P2(z) = 0.5 <3cos2({theta})z - 1>. {theta} is the angle between the membrane normal, directed along the z-axis (Fig. 2), and the dipole moment of water. z is the position of the oxygen atom of water along the channel axis. The orientation of water was further analyzed by calculating the probability distribution p[cos({theta}),z] yielding the likelihood for observing a given water orientation cos({theta}) at position z along the channel axis.



View larger version (136K):
[in this window]
[in a new window]
  		FIGURE 2  The simulated system. Direction of periplasm is upwards and cytoplasm is downwards. GlpF monomers are rendered in cartoon representation in different colors. Aliphatic tails of lipid molecules are rendered in dark yellow with some lipid molecules omitted for clarity; headgroups are shown in sphere representation. Water molecules are shown in red on the two sides of the membrane. The single file of water molecules in one monomer is shown as vdW spheres in red and white.

 
Electrostatic interaction energy
The electrostatic interaction energy U(z) between a water molecule W located at position z along the channel axis and the environment (E; the membrane, other water molecules, the protein, or subsets of the latter) was evaluated as a time-averaged sum of the nontruncated Coulomb interaction between the three atoms of the water molecule and all atomic charges of the environment E in question:

(4)
Averaging involved 400 configurations separated by 10 ps and collected over the last 4 ns of the simulation of system (I), i.e., WT GlpF. zO is the z-coordinate of the oxygen of the water molecule. The channel axis was divided into slabs of ~0.8 Å width.

Electric field
The electric field E(z) generated at position z along the channel axis by the partial charges of the environment E was computed using, as trial field points, the position of the oxygen atoms of water molecules inside the channel:

(5)
where rO and ri are the position of the oxygen atom of the water and the position of the partial charge from the environment, respectively. Averaging and binning of E(z) were performed as in the calculation of U(z) (Eq. 4).

To examine the electrostatic interaction between the environment (or parts thereof) and a trial point dipole moment µ of magnitude (Jorgensen et al., 1983Go), representing a water molecule inside the channel, normal (Uz(z)) and in-plane (Uxy(z)) interaction energies were calculated as

(6)

(7)
where Ex(z), Ey(z), and Ez(z) are the x, y, and z components of the electric field E(z) as obtained from Eq. 5, respectively, and Exy(z)=|Exy(z)| is the in-plane field amplitude. Binning and averaging were again performed as in Eqs. 4 and 5, and the positions of the oxygen atoms of water molecules were used for positioning the trial point dipole along the channel axis z. In calculating Uz and Uxy, the dipole moment of water µ is assumed to be aligned favorably with Ez and Exy, respectively.


   RESULTS AND DISCUSSION
TOP
 ABSTRACT
 INTRODUCTION
 METHODS
 RESULTS AND DISCUSSION
CONCLUSION
 ACKNOWLEDGEMENTS
 REFERENCES
 
In this section, we first briefly describe the conformational fluctuation of the channel during the simulations. Then we discuss water permeation, water translocation mechanism, and water equilibrium configurational behavior. A detailed analysis of the electric field generated by the channel and water-protein electrostatic interaction, which constitutes the most important basis of selectivity in AQPs, is finally provided along with a suggestion of barriers raised against ion conduction through the water pore.

Conformational fluctuations of the channel
Fig. 3 a presents the radius of the water pore along with a snapshot (Fig. 3 b) from the simulation of WT GlpF-G. The latter, drawn at a scale matching the pore radius profile, illustrates the structural details of the channel interior, including Asn68 and Asn203 of the conserved NPA motifs, Arg206 at the SF, and backbone carbonyl groups of the re-entrant loops, which participate in substrate conduction.



View larger version (31K):
[in this window]
[in a new window]
  		FIGURE 3  (a) Channel radius calculated along the channel axis by the program HOLE (Smart et al., 1996Go) for crystal structures of GlpF+G (PDB code 1FX8) and GlpF-G (PDB code 1LDI) and for the simulated WT GlpF-G. For the simulated channel, the fluctuations (standard deviations) are shown as dashed lines. In all structures the narrowest part of the channel is located at the SF region (-5 Å < z < -2 Å). Removal of glycerol from the channel results in a narrowing of the SF and, in the case of the simulated system, a narrowing of the region near the NPA motifs. (b) Snapshot from the MD simulation of WT GlpF-G drawn at a scale matching the pore-radius profile in a. The bipolar orientation of the hydrogen bonded water file precluding proton conduction in GlpF is discernible in b, i.e., opposite orientation of water molecules in the two half channels is observed. This pattern propagates from the central water molecule hydrogen bonded to Asn68/Asn203 at the two NPA motifs (residues 68–70/201–203). These motifs are shown in green with the side chains of Asn68/Asn203 displayed in ball and stick representation. Recognizable are also hydrogen bonds between the water molecules and between water molecules and the carbonyl groups in the nonhelical parts of the two re-entrant loops, which are shown in gray. The exposed carbonyl oxygen atoms of residues 65–66 and 198–200 are displayed as red spheres. The carbonyl groups of Ala201 and His66 are those closest to the NPA motifs. The carbonyl groups of Met202 and Leu67, which participate in hydrogen bonding to the H{delta}s of Asn68 and Asn203, are not shown. The {alpha}-helical parts of the loops are shown in blue. Residues Trp48, Phe200, and Arg206, constituting the SF, appear in light blue; for clarity, nonpolar hydrogen atoms in these residues are omitted. The polarity of the electrostatic dipoles generated by the M3/M7 half-helices is indicated by + (blue) and - (red). Water molecules in the constriction region of the channel are marked by red stars.

 
Comparison of MD simulations of GlpF+G (Jensen et al., 2001Go; Tajkhorshid et al., 2002Go) with WT GlpF-G (Tajkhorshid et al., 2002Go, and present study) indicates that removal of glycerol from the channel results in a slight narrowing of the pore at the SF region. Calculating the radius of the channel in different regions by the program HOLE (Smart et al., 1996Go) indicates that the narrowest region of the channel in WT GlpF-G is 10–15% smaller than that in GlpF+G and is shifted ~2.0 Å along the channel axis toward the channel center (see Fig. 3 a), a finding confirmed when comparing the crystal structure of GlpF+G (PDB-code 1FX8) (Fu et al., 2000Go) with those of WT GlpF-G (PDB-codes 1LDI and 1LDA) (Tajkhorshid et al., 2002Go). After removing glycerol, the narrowest part of the channel is still located in the SF region, and the length of the constriction region (20 Å) where water molecules move concertedly, is preserved. Although the crystal structures of GlpF+G and WT GlpF-G are very similar in the remainder of the channel, we notice in simulated WT GlpF-G, i.e. after removal of glycerol another narrowing of the channel right above the NPA motifs ( Å).

The RMSD values between the crystal structure of GlpF+G (PDB-code 1FX8) and those of WT GlpF-G are 0.26 Å (1FX8 vs. 1LDI) and 0.51 Å (1FX8 vs. 1LDA) (Tajkhorshid et al., 2002Go). The RMSD between the backbone atoms of the simulated WT GlpF-G and the crystal structure converges in each of the four monomers to 0.9–1.5 Å.

Water translocation
In all WT simulations, i.e., systems (I)–(V), and in simulations of the double mutant system, water molecules in the constriction region of the channel moved due to thermal fluctuations. Water molecules form a single file inside the channel and translocate concertedly. This is discernible from Fig. 3 b where a representative snapshot illustrates the highly constrained water configuration in WT GlpF-G. Trajectories of all water molecules visiting the constriction region of each of the four monomers during the 5 ns of the simulation of WT GlpF-G are presented in Fig. 4. Although several features of water motion in the channel are analyzed and quantified in this article, Fig. 4 is the most transparent way of presenting water dynamics in AQPs. The movements of neighboring water molecules in the single file along the channel axis are highly correlated. Water permeation through the channel, therefore, can be viewed as the translocation of the whole single file through the channel, which occurs in hops (shifts) of ~2–3 Å (Fig. 4). The movement of the single file involves simultaneous exchange of protein-water hydrogen bonds, which is a consequence of a narrow pore and has also been observed during glycerol-water co-translocation in GlpF+G (Jensen et al., 2001Go). Defects (water-water interchanges or crossings) in the single file occasionally occur, whereas no persistent disruption of the single file in any part of the channel is evident from the trajectories (Fig. 4).



View larger version (132K):
[in this window]
[in a new window]
  		FIGURE 4  Water translocation in the constriction region of WT GlpF-G during the 5-ns time course of the simulation. The direction of cytoplasm and periplasm is along increasing and decreasing z, respectively. A running average of 50 ps is given to improve clarity. Individual water molecules are assigned different colors each time they enter or re-enter the channel.

 
No preferential direction of water movement inside the channel was observed, i.e., water moved in either direction and no net transmembrane flux occurred, although local density fluctuations together with thermal fluctuations caused 18 water molecules to completely permeate the channel during the 4 ns of analysis. Considering that AQPs are passive water channels, the absence of net diffusion is expected in equilibrium MD simulations with constant isotropic pressure and periodic boundary conditions. For a study of pressure-induced water transport in AQPs, which readily leads to a net water flow across the membrane allowing for appropriate comparison with water flux induced by osmotic or hydrostatic pressure differences in experiments, we refer the reader to the work of Zhu and co-workers (Zhu et al., 2002Go, 2003Go; Zhu and Schulten, 2003Go).

Water-water interchanges, i.e., single file defects, observed in the simulation of WT GlpF are quantified in Table 1. Although water dynamics inside the channel is largely coupled to the motion of groups lining the pore, the radius of the pore seems to be the main structural determinant of such interchanges. Interchanges in the SF (z {cong} -2 Å), i.e., the narrowest region of the channel, are essentially rare. Not surprisingly, interchanges are more likely to occur in wider and more fluctuating regions of the channel ( Å; compare to Fig. 3). Thus one expects that in AQP1, which has a narrower pore (Sui et al., 2001Go), interchanges would be fewer than in GlpF. The large standard deviation found at Å is due to the fact that monomers 1 and 4 exhibit a large number of interchanges in this region (184 and 281, respectively), whereas monomers 2 and 3 do not (46 and 35, respectively).


View this table:
[in this window]
[in a new window]
  		TABLE 1  Water-water interchanges per nanosecond Xww along the channel axis observed over the last 4 ns of the simulation of WT GlpF

 
Interchanges of water molecules inside the channel are associated with a slight steric repulsion between the interchanging water molecules ( kcal/mol) and between the water molecules and their surroundings ( kcal/mol). It has to be mentioned, however, that interchanging water molecules have attractive electrostatic interactions with each other (~-3.9 ± 3.7 kcal/mol) and with their surroundings (~-22.8 ± 6.3 kcal/mol), which can easily overcome steric barriers introduced by the repulsive van der Waals (vdW) interactions. Despite occasions of water interchange, single file arrangement of water is favored by the channel, and in most parts of the simulation time, water molecules avoid crossing each other inside the channel of GlpF.

Channel occupancy
The water occupancy in WT GlpF-G over the 5-ns period of the simulation is shown for each monomer in Fig. 5. The average occupancy n in each monomer is fluctuating closely around the average in monomers 1 and 2, whereas monomers 3 and 4 display more dramatic deviations from the average at certain time instances. However, in all four monomers the average occupancy number n is the same as the overall average, water molecules.



View larger version (48K):
[in this window]
[in a new window]
  		FIGURE 5  Water occupancy within the constriction region of each of the four monomers (M1–M4) in WT GlpF-G during the 5-ns time course of the simulation. For each monomer, the average occupancy n given in each panel leads to the overall average occupancy of water molecules.

 
In the simulation of system (V), in which the channel interior is almost completely hydrophobic, at times the channel is almost instantaneously emptied (but later refilled). This observation illustrates the profound role of carbonyl groups in continuously maintaining the water file inside a narrow channel with an otherwise hydrophobic interior. For a larger pore radius, even a purely hydrophobic channel exhibits rapid water conduction, as demonstrated by simulations of water conduction through carbon nanotubes (Hummer et al., 2001Go; Zhu and Schulten, 2003Go). The sensitivity of a channel's water occupancy to vdW interaction between water and the channel in nanotubes (Hummer et al., 2001Go; Beckstein et al., 2001Go) and the channel emptying observed in system (V) of the present study suggest that a balance between channel radius and hydrophobicity is required for accommodation of water (Beckstein and Sansom, 2003Go).

A largely hydrophobic channel interior in AQPs keeps direct hydrogen bonding with the permeant substrate at a minimum, and ensures rapid water diffusion through the channel, as in hydrophobic carbon nanotubes (Hummer et al., 2001Go; Zhu and Schulten, 2003Go). A slightly more polar interior is found in AQP1 (Sui et al., 2001Go) than in GlpF (Fu et al., 2000Go; Tajkhorshid et al., 2002Go; de Groot et al., 2003Go). Presumably, the more polar interior in AQP1 is necessary since the narrower channel in AQP1 relative to that of GlpF has to be compensated by more polar residues to attract water into the channel. An amphipathic channel lining in GlpF, on the other hand, is a prerequisite for efficient conduction of amphipathic substrates like glycerol (Fu et al., 2000Go; Jensen et al., 2001Go, 2002Go; Grayson et al., 2003Go).

Conduction rates
For correlated water transport in single file, the single parameter determining the bi-directional conduction rate, k±, is the hopping rate k or, equivalently, the mean hopping time describing, respectively, frequency and average time of single file translocations by an amount equal to that of the average water-water spacing (Berezhkovskii and Hummer, 2002Go; Zhu and Schulten, 2003Go). From knowledge of the diffusion coefficient D and of the average water-water separation can be calculated (Berezhkovskii and Hummer, 2002Go; Zhu and Schulten, 2003Go; see also Methods). The average diffusion coefficient of water for WT GlpF-G, as calculated from linear regression of a plot of the mean square displacement vs. time for water molecules within the constriction region of the channel (compare to Eq. 1), is 0.42 x 10-5 cm2 s-1, with lower and upper bounds of 0.23 x 10-5 cm2 s-1 and 0.55 x 10-5 cm2 s-1, respectively. These values for the diffusion constants and an average water spacing Å lead to an average hopping time with lower and upper bounds of and respectively. Using the average channel occupancy reported in the previous section, i.e., , the corresponding bidirectional conduction rate (per monomer) is k± = 1.1 ns-1, with lower and upper bounds of k± = 0.5 ns-1 and k± = 1.6 ns-1, respectively.

For WT GlpF-G, the 18 complete permeation events we observed in 4 ns correspond to k± = 1.1 ns-1 which is identical to the average value for k± derived from the continuous-time random-walk model over the same time interval. This demonstrates the validity of applying the latter model, which was developed for analysis of water permeation through carbon nanotubes (Berezhkovskii and Hummer, 2002Go; Zhu and Schulten, 2003Go), to single file transport through AQPs despite defects (water-water interchanges) that arise occasionally in the single file arrangement (compare to Fig. 4 and Table 1). Our average value for the bi-directional conduction rate and its lower and upper bounds agree well with the experimentally deduced values (Borgnia and Agre, 2001Go; Haymann and Engel, 1999Go; Borgnia et al., 1999Go). However, the average hopping time in WT GlpF-G is 9x longer than (Berezhkovskii and Hummer, 2002Go) and 3x longer than (Zhu and Schulten, 2003Go) as found from similar analysis of single file water conduction through carbon nanotubes (Hummer et al., 2001Go; Zhu and Schulten, 2003Go; both nanotube simulations were at 300 K, but used AMBER and CHARMM force fields, respectively). The slower conduction in WT GlpF-G is well-explained by the fact that the channel lining of this protein is one-half hydrophilic and one-half hydrophobic, whereas purely hydrophobic carbon nanotubes are less strongly interacting with water.

Water distribution
The water file inside the channel extends along the conduction pathway, which is formed by the exposed backbone carbonyl groups of re-entrant loops (residues 64–67 and 198–202) and the side chains of Arg206, Asn68, and Asn203 (Jensen et al., 2001Go; Tajkhorshid et al., 2002Go), as illustrated by the snapshot in Fig. 3 b. Only a few polar hydrogen atoms are involved in stabilizing water during permeation; besides Asn68:H{delta} and Asn203:H{delta} of the NPA motifs, Arg206:H{eta} and Arg206:H{varepsilon} of the SF contribute. Arg206:H{eta} marks the beginning of the one-dimensional constriction region of the channel at Å (approaching here from the periplasmic side; see Fig. 3).

The equilibrium distribution of water inside the channel is essentially of overlapping Gaussians-type and correlates closely with the position of exposed carbonyl groups and hydrogen donor groups lining the channel interior (Tajkhorshid et al., 2002Go). A very similar arrangement of water molecules inside the channel was reported in the crystal structure of WT GlpF-G (Tajkhorshid et al., 2002Go).

The average protein-water hydrogen bonding interactions inside the channel and the water-water hydrogen bonds together defining the integrity of the single file are evident from the atomic radial distribution functions between the channel and water, r(GlpF:O–H2O:H) and r(GlpF:O–H2O:O), and between water, r(H2O:O–H2O:H) and r(H2O:O–H2O:O), respectively (Fig. 6). For systems (I)–(IV), the peaks of r(GlpF:O–H2O:H) and r(H2O:O–H2O:H) are consistently located at 1.9 Å (data for systems (II)–(IV) are very similar to system (I) and are not shown). This demonstrates nearly ideal hydrogen bond distances between the channel lining carbonyls and water and between the nearest-neighbor water molecules. The average hydrogen bonding characteristics are not particularly sensitive to the electrostatic modifications invoked in systems (II)–(IV) and r(GlpF:O–H2O:H) and r(GlpF:O–H2O:O) obtained for these systems are virtually identical to those of system (I). However, after turning off the charges of the channel lining carbonyls in system (V), no structure is featured in r(GlpF:O–H2O:H) and the peak in r (GlpF:O–H2O:O) shifts to a larger value (Fig. 6, bottom panels). This demonstrates that the water-channel interactions define the conduction pathway in GlpF and, at the same time, confirms that the electrostatic interactions between water and the channel lining carbonyls help structuring water inside the channel.



View larger version (21K):
[in this window]
[in a new window]
  		FIGURE 6  Atom-atom radial distribution functions between the carbonyl oxygen atoms of the channel and water (channel-water) and between nearest-neighboring water molecules (water-water) located within the constriction region (compare to Fig. 3). The distribution functions are shown for WT GlpF-G (I) in the top panels, and for system (V) in the bottom panels. Because the spatial arrangement of the single file is linear, no radial normalization was invoked. The distribution function of systems (II)–(IV) are very similar to those of system (I), and are not shown.

 
The interaction between the channel lining and water molecules perturbs to some extent the hydrogen bonding between neighboring water molecules (Tajkhorshid et al., 2002Go; de Groot and Grubmüller, 2001Go). To examine this effect, we calculated for each position along the channel axis a disruption ratio and the corresponding average distance between the oxygen atoms of the neighboring water molecules for systems (I)–(V) and for the W48F/F200T double mutant (see Methods section). The results for WT GlpF-G, system (I), and for the double mutant are shown in Fig. 7. To also provide a dynamic description of the relation of the neighboring water molecules inside the channel, the correlation of their motion was calculated and is also shown in Fig. 7. The results for systems (II)–(V) do not differ significantly from the results obtained for system (I), and are consequently not shown.



View larger version (18K):
[in this window]
[in a new window]
  		FIGURE 7  Average configuration of the single file through the channel for WT GlpF-G (I) (top panels) and for the W48F/F200T double mutant (Tajkhorshid et al., 2002Go) (bottom panels): (left) distance between oxygen atoms of nearest-neighbor water molecules dO–O(z), and disruption ratio dr(z) (dashed lines) resolved along the channel axis; (right) correlation coefficient c(z) between nearest-neighbor water molecules.

 
In general, disruption is very rare, indicating that water molecules predominantly remain joined in single file. In system (I) the disruption reaches a minimum in the SF region coinciding with c(z) reaching a maximum in this region (Fig. 7). In the SF, the side chain of the conserved Arg206 interacts strongly with water molecules. However, this has a minimal effect on the interaction of the water with neighboring water molecules. Other side chains in this region, Phe200 and Trp48, are nonpolar and cannot interact with water strongly. Therefore, neighboring water molecules will not be constrained by hydrogen bonds to multiple sites that could keep them apart from each other. At the exit of the SF toward the center of the channel, however, water can form hydrogen bonds with Arg206:H{varepsilon} (at z = -2.5 ± 0.5 Å) and with Phe200:O (at z = -2.5 ± 0.3 Å), where the first peak in dr(z) and dO–O(z) is observed for WT GlpF-G. The next major disruption peak in system (I) is observed at the NPA region where water molecules hydrogen-bond through their O atom, rather than their H atoms, to the polar side chains of Asn68/Asn203. In the W48F/F200T double mutant both the size and the polarity of the SF are increased. Together, these effects result in the presence of multiple water molecules in the SF region of this mutant, which, in turn, changes dr(z) and c(z) (Fig. 7).

The average nearest-neighbor cross-correlation coefficient c(z) along the channel axis (Fig. 7) quantifies the concertedness of water translocation shown in Fig. 4. The correlation curves in Fig. 7 show that, with no exception, water translocation inside the channel is correlated and that the correlation is highest in the periplasmic half of the channel. Interestingly, for WT GlpF-G, the maximum correlation is exhibited in the SF region. In contrast, for the W48F/F200T double mutant the correlation in this region is reduced due to stronger hydrogen bonding interactions owing to the presence of multiple water molecules at the entrance of the channel, and, possibly, due to the higher polarity of residues introduced by the mutations.

Water orientation
The equilibrium water orientation inside the channel of WT GlpF-G (I), characterized by the order parameters P1(z) = <cos({theta})z> and P2(z) = 0.5<3cos2({theta})z -1>, is presented in Fig. 8, b and a, respectively. In Fig. 8 b, P1(z) is furthermore shown for systems (II)–(IV). The probability p[cos({theta}),z] of a given water orientation cos({theta}) as a function of position z inside the channel of WT GlpF-G is shown in Fig. 8 c. Fig. 8 a illustrates that a high degree of water ordering exists inside the channel of WT GlpF-G, whereas disorder of water is approached at the channel outlets as revealed by P2(z). Inside the channel, P2(z) reaches a minimum in the channel center, i.e., at the NPA region around which P2(z) is fairly symmetric. Together, this is indicating a similar degree of ordering in both the periplasmic and the cytoplasmic regions of the channel. Within the channel, the order parameter profile P1(z) in Fig. 8 b is roughly sigmoidal and fairly similar in systems (I)–(III). P1(z) exhibits an inflection point at Å, i.e., at the midpoint between the two H{delta} atoms of Asn68 and Asn203 [z(Asn68:H{delta}) = 6.1 Å ± 0.3 Å; z(Asn203:H{delta}) = 4.7 Å ± 0.4 Å]. The shape of P1(z) reflects the well-known (Tajkhorshid et al. 2002Go) equilibrium, bipolar water orientation in the two halves of the channel with a hydrogen donor/hydrogen acceptor pattern and, hence, the water dipole inversion around the NPA motifs (compare to Fig. 3 b).



View larger version (35K):
[in this window]
[in a new window]
  		FIGURE 8  Water orientation in the channel. (a) Order parameter P2(z) = 0.5 <3 cos2({theta}) - 1> for WT GlpF-G (I). (b) Order parameter P1(z) = <cos({theta})> for systems (I)–(V). {theta} is the angle between the unit vector approximately aligned along the channel axis/membrane normal, and the water dipole vector. Fluctuations in P2(z) and P1(z) are shown in a and b for WT GlpF-G (I). (c) Histogram p[cos({theta}),z] displaying the probability p for observing a given water orientation cos({theta}) at position z along the channel axis of WT GlpF-G (I).

 
A polarized water configuration, in which the dipole moment of all water molecules in the channel are aligned along the channel axis in the same direction, allows proton conduction through the water wire by the Grotthuss mechanism (Agmon, 1995Go; Tajkhorshid et al., 2002Go; Zhu and Schulten, 2003Go). As shown in Fig. 8 b(IV), to induce such a configuration of water, one would need to simultaneously eliminate both the electric fields generated by helices M3 and M7 and the hydrogen bonding of the conserved Asn68/Asn203 side chains (Tajkhorshid et al., 2002Go). This illustrates a profound role of these conserved residues and of the architecture of the protein in the selectivity of the channel against protons. Additional reduction of the electric influence of the channel on permeating water, realized by neutralizing the channel lining carbonyl groups, results in a faster formation of the polarized water chain, which in the case of present simulations shows an opposite polarity; see Fig. 8 b(V). This finely tuned water orientation therefore effectively blocks proton transport in AQPs (Tajkhorshid et al., 2002Go; Zhu and Schulten, 2003Go), preserving the electrochemical gradient across the cell membrane while allowing rapid water permeation through the pore.

The configuration of water is illustrated in more detail by the probability distribution in Fig. 8 c, indicating that there is essentially vanishing probability for locating the proton conducting orientation in either half of the channel, and only the proton blocking orientation is allowed. The observed bipolar orientation of water inside the channel might be used by a hydroxyl ion defect to propagate to the center of the channel. This defect, however, cannot reach the other side of the channel either, and is destined to return and exit the channel from the same side.

The polar hydrogen atoms in the SF, Arg206:H{eta} and Arg206:H{varepsilon}, interact strongly with the permeating water molecules, thereby increasing P1(z) at -7.5 Å < z < -2.5 Å (Fig. 8 b), and in Fig. 8 c shifting the center of p[cos({theta}), z {cong} -5 Å] toward cos({theta}) = 0. Hence Arg206 makes P1(z) deviate from the order parameter profile that would result if the water molecules formed hydrogen bonds only with the channel lining carbonyl groups. The exact control of the water orientation is a result of precisely optimized electrostatic forces between water and the channel interior. For instance, when eliminating the electrostatic dipoles of the M3/M7 {alpha}-helices (compare to Fig. 3 b) in systems (III)–(V), the deviation of P1(z) at z {cong} 5 Å becomes even more pronounced, since, Arg206:H{eta} and Arg206:H{varepsilon}, in a similar way to the H{delta}s of Asn68/Asn203, tend to orient water molecules with their dipole moment perpendicular to the channel axis, whereas the M7 dipole tends to align the dipole moment along z. However, the net electrostatics of the Arg206-M3/M7 competition implies that only in the channel center can water be oriented perpendicularly to the channel axis. We demonstrate in the following sections why this is the case.

Electrostatic interaction energies
Protein-water interactions are critical in controlling the dynamics and configurational behavior of water in AQPs (Tajkhorshid et al., 2002Go; de Groot et al., 2001Go). In this and the following two sections, protein electric field effects are analyzed in terms of 1), electrostatic stabilization energy resulting from interaction between the protein and the permeating water; 2), orientational constraints imposed onto the water file due to the vectorial nature of the protein electric field; and 3), Coulomb blockade of ion permeation through the water pore.

The electrostatic interaction energy U(z) between water molecules moving along the channel axis and the environment, decomposed into contributions from the protein tetramer, the protein monomer, bulk water within a 20 Å distance of the protein (exclusive of other water molecules in the single file), and the lipid bilayer, are shown in Fig. 9. The dominant part of the interaction energy calculated in Fig. 9 a is naturally due to direct hydrogen bonds between water molecules and the nearby groups, such as the channel lining carbonyl groups, as we discuss further below.



View larger version (33K):
[in this window]
[in a new window]
  		FIGURE 9  Electrostatic interaction energy between a water molecule in the constriction region and its surroundings. (a) Interaction with protein tetramer and embedding monomer. (b) Interaction with bulk water and membrane. (c) Interaction with embedding monomer (same as in a), with other water molecules of the water file in the same monomer, and their sum.

 
The electrostatic interaction of permeating water and the full protein tetramer is essentially identical to that between water and the embedding monomer only (Fig. 9 a). Therefore, the tetrameric architecture of AQPs does not reinforce electrostatic stabilization of the water file. Since other monomers do not influence the permeating water, the interaction with the embedding monomer suffices for our further discussion. A maximal (20 kcal/mol) stabilization of water is encountered at the SF in the periplasmic half of the channel (Fig. 9 a). The location of the energy minimum (z {cong} -3 Å) corresponds to the position of the charged guanidinium group of the conserved Arg206 (z(Arg206:H{varepsilon}) = -2.5 ± 0.4 Å).

Despite the polarity of lipid headgroups and the interfacial water layers, bulk water and the lipid bilayer hardly interact with the water file and, consequently, do not contribute to its stabilization. For instance, the stabilizing effect of the membrane amounts to only 1 kcal/mol (Fig. 9 b). Cross-interactions between water files in different monomers are also negligible, and amount to <0.1 kcal/mol (not shown), which is somewhat surprising since water dipoles in different monomers are favorably aligned relative to each other. In contrast, self-stabilizing electrostatic interactions between the water molecules within the same single file are large, -5 to -17 kcal/mol and, therefore, very important (Fig. 9 c). Interestingly, Fig. 9 c reveals that the interaction between a single water molecule and other water molecules of the same single file is the weakest at the SF region. This observation can be explained by the fact that water interacts most strongly with the protein at this region; the positively charged side chain of the conserved Arg206 has a very strong interaction with a water molecule located in this region, as indicated by the location of the maximal protein stabilization effect in the same region (Fig. 9 c). Apart from the SF region, water-water interaction is stronger than water-protein interaction in all parts of the channel, as clearly demonstrated in Fig. 9 c, and accordingly permits fast water permeation in single file.

The observed strong interaction of water with the protein and its weaker interaction with neighboring water molecules at the SF region are somewhat in line with the proposed disruption of the water-hydrogen-bonded network at this position (de Groot et al., 2001Go). Monitoring the intermolecular distance of adjacent water molecules, however, we find very little disruption at this position (compare to Fig. 7). The observed differences might in part be related to the different force fields applied, particularly with respect to the water models used (the GROMACS force field and the SPC water model were used by the authors in de Groot et al., 2001Go; here we use the CHARMM27 parameter set with the TIP3 water model). However, it should be noted that even transient (ps) association of two fragments of a disrupted water file would permit significant proton conduction, since charge translocation occurs on a subpicosecond timescale (Schulten and Schulten, 1986Go; Pomès and Roux, 1996Go; Roux, 2002Go; Zhu and Schulten, 2003Go). In contrast, propagation of an orientational defect in the proton conduction cycle is orders-of-magnitudes slower (Schulten and Schulten, 1986Go; Pomès and Roux, 1996Go; Roux, 2002Go; Zhu and Schulten, 2003Go), and clearly not possible with the bipolar water configuration found in AQPs.

To examine the importance of the protein architecture and the role of conserved residues in the electrostatic tuning of AQPs, we studied the electrostatic interaction of permeating water with various subsets of the protein, namely the side chains of Arg206, Asn68/Asn203, the channel lining carbonyl groups of the re-entrant loops, and the half-helices M3 and M7 (Fig. 10).



View larger version (27K):
[in this window]
[in a new window]
  		FIGURE 10  (a) Electrostatic stabilization energy of the water molecules along the channel due to the side chains of Asn68 and Asn203, the channel lining carbonyl groups of residues 64–67 and 196–202, and the Arg206 side chain. (b) Backbone contribution of the M3/M7 helices and full residue contribution of these helices, exclusive of the contribution of the side chain of Arg206. (c) Channel with the key groups displayed in colors matching the results for the electrostatic interaction energy in a and b. Hydrogen atoms are omitted for clarity.

 
Direct hydrogen bonds to the carbonyl groups of the re-entrant loops, which constitute the conduction pathway (Jensen et al., 2001Go; Tajkhorshid et al., 2002Go; compare to Figs. 3 b and 6), stabilize individual permeating water molecules by 4–7 kcal/mol in the two halves of the channel (Fig. 10 a). Near the NPA motifs at the center of the channel, the contribution of these carbonyl groups becomes minimal, since no direct hydrogen bond between them and water can arise in this region. Instead of interacting with the substrate, carbonyl groups of residues in this region are engaged in important hydrogen bonds that stabilize the structure of the channel interior (Zhu et al., 2001Go; Jensen et al., 2001Go). At the NPA motifs, water gains a smaller, yet significant, stabilization (2 kcal/mol) through hydrogen bonds with the side chains of Asn68/Asn203 (Fig. 10 a). The most pronounced effect is, however, due to the charged side chain of Arg206 reaching 15 kcal/mol in the periplasmic half of the channel (Fig. 9 c). Near the channel center, Arg206 causes a slight destabilizing effect, which is compensated by the side chains of Asn68/Asn203 stabilizing water at this location (Fig. 10 a).

The subtle balance of electrostatic contributions inside the channel ensures a stable bipolar water configuration in the pore. Furthermore, a minimal interaction between the centrally located water molecule and all of the protein but the Asn68/Asn203 side chains in the channel center ensures unhindered and fast dipole inversion of water at this position. Consequently, simultaneous proton blockage and fast water permeation, not limited by the crucial dipole inversion, is efficiently accomplished. Along these lines, as shown in Fig. 10 b, the resultant electric field due to the backbone atoms (C{alpha}–H{alpha}, C=O, and N–H) of the M3/M7 helices also stabilize the centrally (NPA) located water molecule only marginally, namely by <1 kcal/mol.

The electrostatic interaction energy between water and the helices (mainly M7) in the periplasmic region is positive, yet small. This is due to the unfavorable interaction of water dipoles, which are fixed by strong hydrogen bonding to Arg206, and the field generated by backbone atoms of M7 (compare to Fig. 8, b and c). After inclusion of side chains, except that of Arg206 whose very strong interaction can mask the effect of the others, the interaction of water with M3/M7 in this region becomes attractive (Fig. 10 b). Still though, it is <1 kcal/mol per water molecule since no direct hydrogen bonding interactions are present.

Electrostatic field components
To further analyze the electrostatic tuning of permeating water molecules by the channel, normal (Ez) and in-plane (Exy) components of the electric field E of the protein along the channel axis were calculated. The results are summarized in Fig. 11, a and b. Similar to electrostatic interaction energies (Fig. 10), the electric field was further decomposed in terms of group contributions from the protein. We note in passing tha