INTRODUCTION

TOP ABSTRACT INTRODUCTION METHODS RESULTS AND DISCUSSION CONCLUSIONS REFERENCES

The ability of living cells to transport water, ions, and water-soluble molecules across their cell membrane is mediated by proteins that function as transporters and channels. Water channels conduct water across the membrane and play important roles in cell osmotic regulation. Aquaporins (AQPs) are a family of water channel proteins present in all life forms (Borgnia et al., 1999). All members of the AQP family permeate water, but block the transport of protons (Tajkhorshid et al., 2002). A large variety of AQPs have been identified in plants, which are very dependent on water in their local environment (Johansson et al., 2000) and utilize osmotic pressure for many functions. AQPs are also widely distributed in various organs of the human body, such as kidney, eye, and the brain, and their malfunction has been connected to diseases such as diabetes insipidus and congenital cataracts (Borgnia et al., 1999). Several AQPs have also been characterized in bacteria (Hohmann et al., 2000).

The ability of a membrane to conduct water is characterized by the ratio of net water flow to the hydrostatic or osmotic pressure difference across the membrane. A comprehensive introduction to osmotic permeation can be found in Sperelakis (1998). In the presence of a hydrostatic pressure difference, P, across the membrane, water flows from the high-pressure side to the low-pressure side of the membrane. Under physiological conditions, the respective volume flux J_v (cm/s), defined as the net flow of water (cm³ /s) per unit area of the membrane (cm²) is proportional to P

(1)

where L_P (cm³/N·s) is referred to as the hydraulic permeability.

When the two sides of a membrane have the same hydrostatic pressure but different concentrations of an impermeable solute, an osmotic pressure difference will be established, and water will flow from the side with lower solute concentration to the other side. In dilute solutions, the flux is linearly proportional to the solute concentration difference C

(2)

where J_W (mol/s·cm²) is the molar water flux presented as the number of moles of water passing through the unit area of the membrane per second; C (mol/cm³) is the solute concentration difference; and P_f (cm/s) is defined as the osmotic permeability of the membrane.

In dilute solutions, the water flux produced by a solute concentration difference C is identical to that produced by a hydrostatic pressure difference P = RTC, where R is the gas constant and T is the temperature.

Therefore, L_P and P_f are related by a constant factor

(3)

where V_W is the molar volume of water (18 cm³ /mol).

Water is known to diffuse through lipid bilayers and, hence, all cellular membranes are at least somewhat water-permeable. However, specialized water channels are mainly responsible for the intrinsically high water permeability of certain cellular membranes (Borgnia et al., 1999). Because each channel conducts water independently of other channels, one can define the hydraulic permeability l_P (cm⁵/N · s) and osmotic permeability p_f (cm³/s) for a single water channel as in Eqs. 1 and 2

(4)

(5)

where j_V (cm³/s) and j_W (mol/s) are the volume and molar flux through a single channel, respectively. p_f and l_p obey a similar relation as in Eq. 3, namely,

(6)

In this paper conduction properties of the membrane are denoted by capital letters (e.g., J, L_P, P_f), while the properties of a single channel are denoted by small letters (e.g., j, l_P, p_f). When permeation through other parts of the membrane is negligible, L_P and P_f are equal to the density of a channel (number of channels per unit area of the membrane) times l_P and p_f, respectively. Obviously, p_f (or l_P) is the main physical quantity characterizing a water channel. P_f of a membrane can be measured experimentally, but only if the density of the channel in the membrane is known can p_f be determined.

The channel properties l_P and p_f are determined by the size and architecture of the interior of the channel. Since some of the water channel proteins, e.g., AQPs (Fu et al., 2000; Murata et al., 2000; Sui et al., 2001), are structurally known, it is desired to relate their l_P and p_f values to their structures. Molecular dynamics (MD) simulations are potentially able to provide this structure-function relationship because they can reveal dynamic processes at atomic resolution. The accuracy of the simulations could be tested by comparing the calculated l_P and p_f to observed values in experiments.

In equilibrium MD simulations, only slight net water flow (due to thermal fluctuation) through the channel can be observed, and the results cannot be directly used to calculate the osmotic permeability p_f. Using the steered molecular dynamics (SMD) method (Isralewitz et al., 2001; Izrailev et al., 1998), however, one can apply external forces to water molecules inside the channel to accelerate motion in the desired direction through the channel. To determine the experimentally measurable properties (e.g., p_f) from SMD trajectories, one needs to apply suitable forces that can be related to the hydrostatic pressure difference between the two sides of a membrane. Direct generation of a pressure difference cannot be easily implemented technically in a typical system under periodic boundary conditions, which contains a layer of membrane and a layer of water, the water layers on both sides of the membrane being actually connected into a single layer.

In this paper we present a method to produce in MD simulations a hydrostatic pressure difference across the membrane for systems under periodic boundary conditions, making it possible to computationally observe a net water flow through channels, and to calculate l_P or p_f from simulations.

We demonstrate the applicability of the suggested method through a series of simulations on the Escherichia coli glycerol uptake facilitator (GlpF). GlpF (Heller et al., 1980), a tetrameric membrane protein, is a member of the AQP family (Borgnia et al., 1999) and is known to permeate water (Borgnia and Agre, 2001). We choose this protein because of the availability of high-resolution structures (Fu et al., 2000; Tajkhorshid et al., 2002) that proved very stable in equilibrium MD simulations (Jensen et al., 2001).

	METHODS

TOP ABSTRACT INTRODUCTION METHODS RESULTS AND DISCUSSION CONCLUSIONS REFERENCES

For typical MD simulations under periodic boundary conditions, the unit cell is rectangular. We assume here that the membrane lies in the xy plane of the orthogonal unit cell, and, therefore the z axis is normal to the membrane. The area of the membrane in the unit cell, A, is equal to the product of the x and y dimensions of the cell. We use P₁ and P₂ to denote the hydrostatic pressure at the upper and lower surfaces of the membrane, respectively. In equilibrium MD, P₁ and P₂ are obviously equal to each other.

A water pressure gradient can be produced in MD simulations by applying a constant force f in the z-direction on all water molecules, as illustrated in Fig. 1. Under this constant force field the pressure in the water is no longer uniform, but dependent on the z-position. Here we assume the membrane is fixed in its position, which could be achieved in specific simulations by constraints, or by applying counter forces, as will be described later. Under periodic boundary conditions, the water molecules between two membranes in adjacent unit cells feel three external forces on them: the forces exerted by the two membranes, P₁A and P₂A, and the applied forces, the sum of which is nf (n is the total number of water molecules in a unit cell). In a stationary state, these three forces are balanced, i.e., P₁A + P₂A + nf = 0. Therefore, the pressure difference across the membrane can be related to the applied force by the formula:

(7)

The water flux through the channel in the membrane can be easily measured by counting the water molecules passing through the channel during the simulation. Thus, this method allows one to quantitatively calculate l_P or p_f, which can be compared with experimental values.

View larger version (45K): [in this window] [in a new window]   FIGURE 1   Illustration of the method to produce a pressure gradient under periodic boundary conditions. The membrane and water molecules outside the unit cell are the "images" of those inside. A constant force f, shown by small arrows, is exerted on all water molecules.

Unlike the case of conventional SMD simulations, in which only a small number of atoms are pulled (Isralewitz et al., 2001), the hydrostatic pressure established in the present method is generated by pulling a large number of water molecules and, therefore, the new simulations mimic macroscopic hydrostatic pressure in experiments. Furthermore, one does not need to know or assume the mechanism of water passage to set up the simulations, and the system will determine by itself which water molecules enter or move through the channel. Thus one can observe at the atomic level a permeation event that is similar to what happens in reality.

We tested the method on the GlpF channel with two sets of calculations (referred to as set 1 and set 2), each including four independent simulations. The starting configuration for all simulations was an equilibrated system (shown in Fig. 2) that included the GlpF tetramer, a patch of palmitoyloleoyl phosphatidylethanolamine (POPE) lipid bilayer, and ~17,000 water molecules. The whole system contains 106,189 atoms. For a detailed description of system build-up and equilibration we refer the reader to Jensen et al. (2001). We note that in the present paper, the +z direction is defined as going from the periplasmic side to the cytoplasmic side of the membrane, or pointing downward in Fig. 2.

View larger version (112K): [in this window] [in a new window]   FIGURE 2   Side view of the unit cell including the GlpF tetramer, POPE lipid molecules, and water molecules. The protein is shown in tube representation, lipids in line representation (hydrogen atoms not shown), phosphorus atoms of lipids are drawn as VDW spheres, and water molecules are shown in line representation (top: periplasmic side; bottom: cytoplasmic side).

In all simulations, because the water molecules had external forces on them, the membrane needed to be kept at its position to avoid translation of the system along the direction of the external forces. This was done in our simulations by applying constant forces in the opposite direction on all heavy atoms of the lipid molecules and on the C atoms of the protein, so that the total external force on the whole system was zero. For each simulation the total counter force on the membrane was equal to the total external force on water, and was distributed between the protein and the lipids according to the ratio of their estimated areas in the membrane plane. Finally, the counter forces on the protein and on the lipids were distributed evenly among C atoms and among lipid molecules, respectively.

In the first set of simulations (set 1), external forces were applied to the oxygen atoms of all water molecules, including those inside the channels. Application of a different force in each simulation resulted in different hydrostatic pressure gradients. In two of the simulations, a pressure difference of ~200 MPa was induced in +z and z directions, respectively; in the other two simulations, the induced pressure difference was ~400 MPa.

Forces on water molecules inside the channel might artificially increase the measured conduction rate. Because we want to describe the conductivity of water induced only by the pressure difference across the membrane, and not directly by forces on the water molecules inside the channel, we performed another set of four simulations (set 2), in which we used the same external forces as in set 1, but water molecules inside the channels were excluded from force application. In set 2 we defined a large enough rectangular "exclusion region," which completely excluded all vestibules and constriction regions of the four channels from external forces. In every simulation step, external forces were applied only on water molecules located outside the mentioned region. This ensured that water molecules inside and near the channels were not subject to the external force. For example, all of the water molecules shown in Fig. 3 are exempt from external forces at that simulation step (the water molecules can become subject to external forces only when they leave the channel region).

View larger version (55K): [in this window] [in a new window]   FIGURE 3   A GlpF monomer with channel water. The protein is shown in tube representation. Water molecules located in the channel and the vestibules are drawn as VDW spheres. These water molecules are all inside the exclusion region defined in set 2 of the simulations.

The first 50 ps of each simulation was discarded, and the rest of the trajectory was used for analysis. To calculate the water flux through a channel, a plane normal to z in the channel was defined, and the net water flow was evaluated by counting the number of water molecules crossing the plane (+1 if water crossed the plane downward, and 1 if upward). In our analysis we selected two such planes in the central part of the channel, and used the average count to determine the water flux.

The simulations were performed under periodic boundary conditions, with fixed volume and constant temperature (310 K) achieved by Langevin dynamics. Full electrostatics calculation was done using the particle mesh Ewald (PME) method (Essmann et al., 1995). All simulations were performed using the CHARMM27 force field (MacKerell, Jr., et al., 1998; Schlenkrich et al., 1996), the TIP3P (Jorgensen et al., 1983) water model, and the MD program NAMD2 (Kalé et al., 1999). The overall simulation time was ~5 ns, and different simulations were run on the supercomputing centers at NCSA and Pittsburgh, and on a local Linux cluster; 1 ns of simulation took ~13 days on 64 Cray T3E processors at Pittsburgh, 6.3 days on 80 Origin2000 processors at NCSA, or 8.2 days on the 32 processor Linux cluster.

	RESULTS AND DISCUSSION

TOP ABSTRACT INTRODUCTION METHODS RESULTS AND DISCUSSION CONCLUSIONS REFERENCES

During the simulations with the higher pressure difference (400 MPa), the surfaces of the lipid bilayer became more irregular and less flat than in equilibrium simulations. Such a perturbation was not obvious in the simulations applying the lower (200 MPa) pressure difference. Thus, for simulations in which a very high pressure difference is induced, one may want to apply to the lipid molecules harmonic constraints in the z dimension, instead of constant counter forces, to ensure the stability of the lipid bilayer.

The overall structure of the protein appeared very stable in all the simulations. However, in the simulations with the higher pressure difference (400 MPa), in some monomers, a few water molecules moved into the internal cavities of the protein. This behavior is in agreement with the proposed mechanism of pressure-induced denaturation of proteins (Hummer et al., 1998). In particular, in one of the simulations in which the periplasmic side had higher pressure, a water molecule entered the space between the side chain of Arg₂₀₆ and the main chain of Phe₁₃₅ in one of the monomers and broke the two H-bonds between the guanidinium group of Arg₂₀₆ and the backbone oxygens of these two residues. Not being stabilized by the H-bonds, and due to the insertion of the water molecule, the long side chain of Arg₂₀₆ was then displaced from its original position and blocked the channel. To preserve the correct position of Arg₂₀₆, we repeated this simulation in the presence of harmonic constraints preserving the mentioned two H-bonds in each monomer. Such side chain displacement caused by water was not observed in the simulations applying the lower pressure difference (200 MPa).

During the simulations, as expected, the applied forces produced a water density gradient in bulk water, which is shown in Fig. 4. When the forces on water are in the +z direction, the water density increases with z in bulk water; when the forces are along z, the density decreases with z.

View larger version (16K): [in this window] [in a new window]   FIGURE 4   Water density distribution along the z-direction in the bulk water region. Due to the periodicity of the system, the bulk water is sandwiched between two membranes in the adjacent periodic cells (refer to Fig. 1). The cytoplasmic side of the membrane is on the left, and the periplasmic side is on the right side of the figure. Data points marked by circles and stars are taken from two simulations, where an external force of 1.54 pN along +z and z was applied on all water molecules, respectively. The density is measured by averaging the number of water molecules within a slab of 2 Å thickness over 100 frames taken from the last 100 ps of the trajectory. The standard deviation of the frames is shown as error bars.

The water density difference was discernible in the vestibules of the channel, as shown in Fig. 5. When the external forces on water were directed in the +z direction (Fig. 5, middle panel), the hydrostatic pressure above the membrane was higher than below; consequently, more water molecules were found in the periplasmic vestibule of the channel, and fewer appeared in the cytoplasmic side. This density gradient resulted in a net water flow through the channel along the +z direction. Similarly, when the forces on water were directed in the z direction (Fig. 5, right panel), more water molecules gathered in the cytoplasmic vestibule of the channel, and a net water flow along the z direction was observed.

View larger version (44K): [in this window] [in a new window]   FIGURE 5   Snapshots from three simulations of GlpF. The +z direction is pointing downward in this figure. The channel is represented by two half-membrane spanning repeats, each including a short helix followed by a loop, drawn in tube representation. Water molecules inside and near the channel are shown in CPK representation. The major H-bond partners of water molecules, namely the side chains of Asn68 and Asn203 of the NPA motifs, the side chain of Arg206 in the selectivity filter, and the backbone oxygens of the loops, are drawn in licorice representation. A detailed description of protein-substrate interactions in GlpF is given in Jensen et al. (2001). The direction of the applied forces on water molecules is shown by arrows. Left: equilibrated system without any external force; middle: a downward force is applied on every water molecule, making the hydrostatic pressure (and thus the water density) of the periplasmic side of the channel higher than that of the cytoplasmic side; right: an upward force is exerted on every water molecule, with the same magnitude as in the case above, producing a higher pressure (and thus a higher water density) in the cytoplasmic side of the channel.

Under equilibrium conditions, the water molecules in the channel usually formed a single file, in agreement with previous simulations of GlpF (Jensen et al., 2001; Tajkhorshid et al., 2002) and aquaporin-1 (AQP1) (Zhu et al., 2001). After the application of the pressure difference, however, the channel appeared to adopt more water molecules in the part connected to the high-pressure side. This was especially noticeable in the part of the channel located between the NPA (Asn-Pro-Ala) motifs and the cytoplasmic side of GlpF (the lower parts in Fig. 5): when the cytoplasmic side had a higher pressure, the increased water density on that side increased the number of water molecules entering the channel, and they were disordered, i.e., no longer in single file. This behavior probably arose because this part of the channel is close to bulk water in the cytoplasmic side and has a relatively large capacity to accommodate more water molecules when water density on that side increases. However, water molecules mainly formed single file in the part of the channel located between the NPA motifs and the periplasmic side (the upper parts in Fig. 5), even in the presence of an increased water density on that side. This may be due to the fact that this part of the channel is more constricted, so that it cannot hold more water molecules than a single file under such a density. The so called "selectivity filter" of GlpF, located in this part, including the highly conserved Arg₂₀₆, may contribute to the constriction of the channel.

Tables 1 and 2 show the data obtained from simulation sets 1 and 2, respectively. In each simulation, a different force (in either magnitude or direction) was applied to water molecules, inducing different pressure gradients. The average of counts from the four monomers was used to calculate the water flux in the simulation, and their standard deviation gave an estimate of the fluctuation of the water flux. These values have been plotted versus the applied pressure difference in Fig. 6. For each set a line with the best-fit slope is drawn through the data points. From the resulting slope the hydraulic permeability of a single channel can be determined. The respective values are l_p = (2.0 ± 0.3) × 10¹⁷ cm⁵/N·s and l_p = (9.5 ± 0.7) × 10¹⁸ cm⁵/N·s for sets 1 and 2, respectively. Applying Eq. 6 (T = 310 K) one also obtains the osmotic permeability of a single channel, p_f = (2.9 ± 0.4) × 10¹³ cm³/s for set 1 and p_f = (1.4 ± 0.1) × 10¹³ cm³/s for set 2. The results reveal that when forces were not applied on the water molecules inside the channels, the induced flux decreased by a factor of one-half. Because in cellular membranes water transport is induced by the macroscopic effect of external pressure, simulations of set 2 provide a more faithful description of water conductivity in the channels.

View larger version (12K): [in this window] [in a new window]   FIGURE 6   The dependence of water flux on the applied pressure difference. The horizontal and vertical axes represent the pressure difference P (in MPa) and the water flux j (number of water molecules per monomer per ns), respectively. Data points from sets 1 and 2 described in Tables 1 and 2 are shown with circles and stars, respectively. Error bars are the standard deviations of the water flux. Also shown are two lines with the best-fit slope for the two sets of simulations. According to Eq. 4, the line for set 1 (with force on water molecules inside the channel) corresponds to a hydraulic permeability lp of 2.0 × 1017 cm5/N·s with an estimated standard deviation of 0.3 × 1017 cm5/N·s; the line for set 2 (without force on water molecules inside the channel) corresponds to a hydraulic permeability lp of 9.5 × 1018 cm5/N·s with an estimated standard deviation of 0.7 × 1018 cm5/N·s.

No experimental data of single channel permeability, p_f, have been published for GlpF, so the calculated value of p_f cannot be directly compared to experiments. However, p_f measurements are available for AQP1, another member of the AQP family. Different experiments have reported different p_f values for AQP1 (Walz et al., 1994; Zeidel et al., 1992, 1994), with p_f = 5.43 × 10¹⁴ cm³/s considered to be the most accurate one (Walz et al., 1994). It is also reported that AQP1 has a higher permeability than GlpF (Calamita, 2000). Therefore, our calculations have apparently overestimated the value of p_f by a factor of three or more.

It is noteworthy in this regard that the induced pressure differences in our simulations are significantly larger than the osmotic pressure used in experiments. The osmotic pressure of physiological solutions is usually below 10 MPa (e.g., a 200 mM solution of sucrose has an osmotic pressure of ~0.5 MPa), whereas the pressure differences applied in our simulations were 200 and 400 MPa, i.e., more than an order of magnitude higher. The reason for applying such high pressures is that at a simulation time of <1 ns, the low pressures would yield a very low count of conducted water molecules that would not rise above the statistical error.

It has not been experimentally tested for water channels whether at the large pressures applied here Eq. 4 holds, i.e., whether the water flux is still linearly proportional to pressure differences in this range. Under small pressure differences, the number and configuration of water molecules inside the channel should remain the same as in equilibrium, but due to the high pressures applied in our simulations we observed notably more water molecules in the channel (Fig. 5, right panel). This accumulation of channel water in the cytoplasmic vestibule may influence the linear relation between water flux and hydrostatic pressure difference. The data points at only four pressure differences do not permit a convincing test of the linearity of the flux-pressure relationship. Since it is known that some systems respond linearly to very high perturbations, the permeability observed at high pressures may be extrapolated to physiological pressure difference. Simulations at smaller pressure differences are needed to test this, but will require a 10-fold longer simulation time due to lower conduction water counts. Presently, it is practically unfeasible to achieve much longer simulation times for this relatively large system. In fact, the simulations reported here consumed more than an equivalent of 60 days of computation time on 64 processors of a Cray T3E.

The present method actually induces a pressure gradient in the bulk water rather than a pressure step across the membrane, and assumes that the water flux is determined by the difference between the hydrostatic pressures adjacent to the two membrane surfaces. Because in experiments the bulk waters on the two sides of the membrane have different, but uniform pressures, the validity of a comparison of calculated and observed permeabilities is not guaranteed.

	CONCLUSIONS

TOP ABSTRACT INTRODUCTION METHODS RESULTS AND DISCUSSION CONCLUSIONS REFERENCES

A method for inducing a hydrostatic pressure gradient across a membrane in MD simulations has been described. Due to the relatively small water counts during the simulation time, the accuracy of the results is limited; the application of large pressure differences, which are needed to allow statistically significant results, may lead to deviation from the permeability determined under physiological conditions. However, the method appears to be promising for future studies of water permeation in membrane channels that overcomes the present limitations in computational power. Longer simulations with lower pressure differences may be realized in such studies and give more accurate permeabilities. In view of the increasing power of massively parallel computers that are becoming available, longer simulations can be performed in the near future, where a pressure difference as low as those used in experiments can be applied, thus eliminating the possible deviation of the water flux from linearity. Furthermore, water channels with higher permeabilities (e.g., AQP1) will serve better for this purpose because they require a smaller pressure difference to obtain the same water flux. Simulations with both higher and lower pressure differences could also reveal the linear response region of the channel.

The present method of calculating p_f may be used in the future as an alternative to experimental measurements, since the single channel permeability p_f of many membrane proteins has not been experimentally determined, partly due to the difficulty of estimating the number of channels in the membrane. Furthermore, the simulations may be used to compare the permeability of various water channels, such as different members of the AQP family, or to predict the effect of genetic mutation on permeability. The method can also be used to study the influence of hydrostatic pressure differences on water conduction in artificial water channels (Hummer et al., 2001).