New Research In
Physical Sciences
Social Sciences
Featured Portals
Articles by Topic
Biological Sciences
Featured Portals
Articles by Topic
 Agricultural Sciences
 Anthropology
 Applied Biological Sciences
 Biochemistry
 Biophysics and Computational Biology
 Cell Biology
 Developmental Biology
 Ecology
 Environmental Sciences
 Evolution
 Genetics
 Immunology and Inflammation
 Medical Sciences
 Microbiology
 Neuroscience
 Pharmacology
 Physiology
 Plant Biology
 Population Biology
 Psychological and Cognitive Sciences
 Sustainability Science
 Systems Biology
Liquid–vapor oscillations of water in hydrophobic nanopores

Edited by Peter C. Agre, The Johns Hopkins University School of Medicine, Baltimore, MD (received for review November 10, 2002)
Related Article
Abstract
Water plays a key role in biological membrane transport. In ion channels and waterconducting pores (aquaporins), onedimensional confinement in conjunction with strong surface effects changes the physical behavior of water. In molecular dynamics simulations of water in short (0.8 nm) hydrophobic pores the water density in the pore fluctuates on a nanosecond time scale. In long simulations (460 ns in total) at pore radii ranging from 0.35 to 1.0 nm we quantify the kinetics of oscillations between a liquidfilled and a vaporfilled pore. This behavior can be explained as capillary evaporation alternating with capillary condensation, driven by pressure fluctuations in the water outside the pore. The freeenergy difference between the two states depends linearly on the radius. The freeenergy landscape shows how a metastable liquid state gradually develops with increasing radius. For radii > ≈0.55 nm it becomes the globally stable state and the vapor state vanishes. Onedimensional confinement affects the dynamic behavior of the water molecules and increases the self diffusion by a factor of 2–3 compared with bulk water. Permeabilities for the narrow pores are of the same order of magnitude as for biological water pores. Water flow is not continuous but occurs in bursts. Our results suggest that simulations aimed at collective phenomena such as hydrophobic effects may require simulation times >50 ns. For water in confined geometries, it is not possible to extrapolate from bulk or short time behavior to longer time scales.
Channel and transporter proteins control flow of water, ions, and other solutes across cell membranes. In recent years several channel and pore structures have been solved at near atomic resolution (1–6), which together with three decades of physiological data (7) and theoretical and simulation approaches (8) allow us to describe transport of ions, water, or other small molecules at a molecular level. Water plays a special role here: it either solvates the inner surfaces of the pore and the permeators (for example, ions and small molecules like glycerol), or it is the permeant species itself as in the aquaporin family of water pores (9–11) or in the bacterial peptide channel gramicidin A, whose water transport properties are well studied (12–14). Thus, a better characterization of the behavior of water would improve our understanding of the biological function of a wide range of transporters. The remarkable water transport properties of aquaporins [water is conducted through a long (≈2 nm) and narrow (≈0.3 nm diameter) pore at bulk diffusion rates while at the same time protons are strongly selected against] are the topic of recent simulation studies (15, 16).
The shape and dimensions of biological pores and the nature of the porelining atoms are recognized as major determinants of function. How the behavior of water depends on these factors is far from understood (17). Water is not a simple liquid because of its strong hydrogen bond network. When confined to narrow geometries like slits or pores it displays an even more diverse behavior than already shown in its bulk state (18, 19).
A biological channel can be crudely approximated as a “hole” through a membrane. Earlier molecular dynamics (MD) simulations showed pronounced layering effects and a marked decrease in water selfdiffusion in infinite hydrophobic pore models (20, 21). Recently, water in finite narrow hydrophobic pores was observed to exhibit a distinct twostate behavior. Either the cavity is filled with water at approximately bulk density (liquidlike) or it is almost completely empty (vaporlike) (22, 23). Similar behavior was seen in Gibbs ensemble Monte Carlo simulations in spherical (24) and cylindrical pores (25).
In our previous simulations (23) we explored model pores of the dimensions of the gating region of the nicotinic acetylcholine receptor. Hydrophobic pores of radius R ≥ 0.7 nm were filled during the whole simulation time (up to 6 ns) whereas narrow ones (R ≤ 0.4 nm) were permanently empty. Changing the pore lining from a hydrophobic surface to a more hydrophilic (polar) one rendered even narrow pores water, and presumably ion, conducting. At intermediate radii (0.4 nm < R < 0.7 nm) the porewater system was very sensitive to small changes in radius or character of the pore lining. In a biological context, a structure close to the transition radius would confer the highest susceptibility to small conformational rearrangements (i.e., gating) of a channel.
We have extended these simulations to beyond 50 ns to explore the longer timescale behavior of the waterpore system. Starting from the observed oscillations in water density (Fig. 1a) we analyze the kinetics and the free energy of the system. We compare the water transport properties of the pores to experimental and theoretical data.
Methods
Model. The pore model was designed to mimic the dimensions of a biological pore [e.g., the gate region of the nicotinic acetylcholine receptor (26)], while retaining the tractability of a simple model. Cylindrical pores of finite length were constructed from concentric rings of pseudo atoms (Fig. 1b). These hydrophobic pseudo atoms have the characteristics of methane molecules, i.e., they are uncharged Lennard–Jones spheres with a van der Waals radius of 0.195 nm. The pore consists of two mouth regions (radius R_{M} = 1 nm, length L_{M} = 0.4 nm) and an inner pore region (L_{P} = 0.8 nm) of radius 0.35 nm ≤ R ≤ 1 nm, which is the minimum van der Waals radius of the cavity. The centers of the porelining atoms are placed on circles of radius R + 0.195 nm. The model was embedded in a membrane mimetic, a slab of pseudo atoms of thickness ≈1.5 or 1.9 nm. Pseudo atoms were harmonically restrained to their initial position with a force constant of 1,000 kJ·mol^{1}·nm^{2}, resulting in positional fluctuations of ≈0.1 nm, comparable to those of pore lining atoms in membrane protein simulations although this does not allow for global collective motions as in real proteins.
Simulation Details. MD simulations were performed with GROMACS 3.0.5 (27) and the simple point charge (SPC) water model (28). The Lennard–Jones parameters for the interaction between a methane molecule and the water oxygen are ε_{CO} = 0.906493 kJ·mol^{1} and σ_{CO} = 0.342692 nm from the GROMACS force field. The integration time step was 2 fs and coordinates were saved every 2 ps. With periodic boundary conditions, longrange electrostatic interactions were computed with a particle mesh Ewald method [real space cutoff 1 nm, grid spacing 0.15 nm, fourthorder interpolation (29)], whereas the shortrange van der Waals forces were calculated within a radius of 1 nm. The neighbor list (radius 1 nm) was updated every 10 steps.
Weak coupling algorithms (30) were used to simulate at constant temperature (T = 300 K, time constant 0.1 ps) and pressure (P = 1 bar, compressibility 4.5 × 10^{5} bar^{1}, time constant 1 ps) with the x and y dimensions of the simulation cell held fixed at 6 nm (or 3.9 nm for the 80ns simulation of the R = 0.35nm pore). The length in z was 4.6 nm in both cases (ensuring bulklike water behavior far from the membrane mimetic).
The large (small) system contained ≈700 (300) methane pseudo atoms and 4,000 (1,500) SPC water molecules. Simulation times T_{sim} ranged from 52 to 80 ns (altogether 460 ns). Bulk properties of SPC water were obtained from simulations in a cubic cell of length 3 nm (895 molecules) with isotropic pressure coupling at 300 K and 1 bar for 5 ns.
Analysis.Time courses and density. For the density time courses (Fig. 1a) the pore occupancy N(t), i.e., the number of water molecules within the pore cavity (a cylinder of height L_{P} = 0.8 nm containing the porelining atoms) was counted. The density n(t) is given by N(t) divided by the wateraccessible pore volume V = L_{P} πR_{eff}^{2} and normalized to the bulk density of SPC water at 300 K and 1 bar (n_{bulk} = 53.67 ± 0.03 mol·liter^{1}). The effective pore radius for all pores is R_{eff} = R  δR. Choosing δR = 0.03 nm fixes the density 〈N〉/V in the (most bulklike) R = 1.0nm pore at the value calculated from the radial density, , where R_{0} = 1.05 nm is the radius at which n(r) vanishes.
The density n(r) was determined on a grid of cubic cells with spacing 0.05 nm. Two and onedimensional densities were computed by integrating out the appropriate coordinate(s). A probabilistic interpretation of n(r) leads to the definition of the potential of mean force (PMF) of a water molecule βF(r) = ln[n(r)/n_{bulk}] with β^{1} = k_{B}T, via Boltzmann sampling of states.
Freeenergy density and chemical potential. The Helmholtz free energy as a function of the pore occupancy N at constant T = 300 K for a given pore with volume V was calculated from the probability distribution p(N) of the occupancy as βF(T, V, N) = ln p(N) and transformed into a freeenergy density f(T, n) = F/V. A fourthorder polynomial in n was leastsquarefitted to βf(T, n). The chemical potential μ(T, n) = ∂f(T, n)/∂n was calculated as the analytical derivative of the polynomial.
The βf curves obtained for different radii R from the simulations are only determined within an unknown additive constant f_{0}(T; R) but a thermodynamic argument shows that all of these curves coincide at n = 0. For n = 0 no water is in the pore, so the freeenergy differential is simply dF = S dT + 2γ dA with the constant surface tension of the vacuum (inside the pore)–water (outside) interface of area A = πR^{2}. At constant T this implies F(R) = 2γA + const(T), so that the freeenergy density f(T, n = 0; R) = F(R)/V = 2γA/(LA) = 2γ/L of a pore with radius R and length L is independent of R. Hence, all freeenergy density curves necessarily coincide at n = 0 and f_{0} is a function of T only.
Kinetics. The time series n(t)/n_{bulk} of the water density in the pore was analyzed in the spirit of singlechannel recordings (31) by detecting open (high density; in the following denoted by a subscript o) and closed (approximately zero density; subscript c) pore states, using a Schmitt trigger with an upper threshold of 0.65 and a lower threshold of 0.15. A characteristic measure for the behavior of these pores is the openness 〈ω〉 = T_{o}/T_{sim}, i.e., the probability for the pore being in the open state (23) with errors estimated from a blockaveraging procedure (27). The distribution of the lifetimes t_{α} of state α = {o, c} are exponentials (data not shown). The maximumlikelihood estimator for the characteristic times τ_{o} and τ_{c} is the mean τ_{α} = 〈t_{α}〉.
The freeenergy difference between the closed and open states, ΔF = F_{c}  F_{o}, can be calculated in two ways. First, we obtained it from the equilibrium constant K = T_{c}/T_{o} = (T_{sim}  T_{o})/T_{o} = 〈ω〉^{1}  1 of the twostate system as βΔF_{kin} = ln K. Second, ΔF was determined from p(N) as the ratio between the probability that the pore is in the closed state and the probability for the open state: βΔF_{eq} = ln P_{c}/P_{o} = ln Σ_{N}_{≤}_{N}_{‡} p(N)/Σ_{N}_{>}_{N}_{‡} p(N). The definition of state used here is independent of the kinetic analysis. It only depends on N^{‡}, the pore occupancy in the transition state, the state of lowest probability between the two probability maxima that define the closed and open states. The relationship involving K can be inverted to describe the openness in terms of ΔF(R), 〈ω(R)〉 = (1 + exp[βΔF(R)])^{1}.
Dynamics. The three components of the selfdiffusion coefficient were calculated from the Einstein relations 〈(x_{i}(t)  x_{i}(t_{0}))^{2}〉 = 2D_{i} (t  t_{0}). The simulation box was stratified perpendicular to the pore axis with the central layer containing the pore. During T_{sim} the mean square deviation (msd) of water molecules in a given layer was accumulated for 10 ps and after discarding the first 2 ps, a straight line was fit to the msd to obtain D_{i}. These diffusion coefficients were averaged in each layer for the final result.
The current density (flux per area) was calculated as j_{0} =Φ_{0}/A from the equilibrium flux Φ_{0} = M/T_{sim} with the total number of permeant water molecules M and the effective pore cross section A = πR_{eff} ^{2} for pores or A = L_{x}L_{y} for the bulk case, i.e., a simulation box of water with periodic boundary conditions. Permeant water molecules were defined as those whose pore entrance and exit zcoordinate differed. In addition, distributions of permeation times were computed.
Results and Discussion
The water density in the pore cavity oscillates between an almost empty (closed) and filled (open) state (Fig. 1a). We refer to the waterfilled pore state as open because such a pore environment would favorably solvate an ion and conceivably allow its permeation. Conversely, we assume that a pore that cannot sustain water at liquid densities will present a significant energetic barrier to an ion. As shown in Fig. 1c, water molecules can pass each other and often permeate the pore in opposite directions simultaneously.
Even though the oscillating behavior was already suggested by earlier 1ns simulations (23) only at these longer times do clear patterns emerge. The characteristics of the porewater system change substantially with the pore radius. The oscillations (Fig. 1a) depend strongly on the radius. The water density (Fig. 2) shows large pores to be waterfilled and strongly layered at bulk density. With decreasing radius the average density is reduced because of longer closed states even though layer structures remain. The narrowest pores appear almost void of water.
The sudden change in behavior is borne out quantitatively by the openness (Fig. 3a), which indicates a sharp increase with increasing radius around R = 0.55 nm. Although the range of radii over which this transition takes place appears to be small (0.45–0.7 nm) the crosssectional area doubles. The maximum number of water molecules actually found in the cavity in our simulations more than doubles from 21 to 46 in this range of R, so that the average environment that each water molecule experiences changes considerably.
Density. The radial densities in Fig. 2 show destabilization of the liquid phase with decreasing pore radius. Above R = 0.45 nm distinctive layering is visible in the pore, and for the larger pores appears as an extension of the planar layering near the slab. For R < 0.45 nm no such features remain and the density is on average close to 0. The open state can be identified with liquid water and the closed state with water vapor. In the continuously open 1nm pore, the average density 〈n(t)〉/n_{bulk} is 0.82, but 0.032 in the closed 0.35nm pore. Brovchenko et al. (33) carried out Gibbs ensemble Monte Carlo simulations of the coexistence of liquid TIP4P water with its vapor in an infinite cylindrical hydrophobic pore of radius R = 1.075 nm. At T = 300 K they obtained a liquid density of 0.81 and a vapor density close to 0, in agreement with the numbers from our MD simulations.
Analysis of the structure in the radial PMF (Fig. 3b) lends further support to the above interpretation. Water molecules fill the narrow pores (R ≤ ≈0.45 nm) homogeneously as it is expected for vapor. For the wider pores, distinct layering is visible as the liquid state dominates. The number of layers increases from two to three, with the central water column being the preferred position initially. As the radius increases, the central minimum shifts toward the wall. For R = 0.7 nm the center of the pore is clearly disfavored by 0.2 k_{B}T. In the largest pore (R = 1 nm), the influence of curvature on the density already seems to be negligible as it is almost identical to the situation near a planar hydrophobic slab.
Kinetics. Condensation (filling of the pore) and evaporation (emptying) occur in an avalanchelike fashion as shown in Fig. 1c. In our simulations both events take place within ≈30 ps, roughly independent of R. However, the actual evaporation and condensation processes seem to follow different paths, as we can infer from the analysis of the kinetics of the oscillations. The time series of Fig. 1a reveals that the lifetimes of the open and closed state behave differently with increasing pore radius (Fig. 3c). In the range 0.45 nm ≤ R ≤ 0.6 nm, the average time a pore is in the closed state is almost constant, τ_{c} = 1.40 ± 0.37 ns; outside this range no simple functional relationship is apparent. The average open time can be described as an exponential τ_{o}(R) = aexp(R/ζ) with a = 1.3 × 10^{5} ns and ζ = 4.9 × 10^{2} nm for 0.35 nm ≤ R ≤ 0.7 nm.
1/τ_{o} is related to the “survival probability” of the liquid state and 1/τ_{c} to that of the vapor state. These times characterize the underlying physical evaporation and condensation processes. Their very different dependence on R implies that these processes must be quite different. The initial condensation process could resemble the evaporation of water molecules from a liquid surface. Evaporating molecules would not interact appreciably, so that this process would be rather insensitive to the area of the liquid–vapor interface A = πR^{2} and hence R. The disruption of the liquid pore state, on the other hand, displays very strong dependence on the radius. Conceivably, the pore empties once a density fluctuation has created a vapor bubble that can fill the diameter of the pore and expose the wall to vapor, its preferred contact phase. The probability for the formation of a spherical cavity of radius λ with exactly N water molecules inside was determined by Hummer et al. (34). From their study we find that the probability p(λ; n) for the formation of a bubble of radius λ and density below a maximum density n is apparently an exponential. Once a bubble with λ ≈ R develops, the channel rapidly empties but this occurs with a probability that decreases exponentially with increasing R, which corresponds to the observed exponential increase in τ_{o}. In particular, for lowdensity bubbles (n < 0.2 n_{bulk}) we estimate the decay constant in p(λ; n) as 2 × 10^{2} nm, which is of the same order of magnitude as ζ.
From the equilibrium constant K(R) = T_{c}(R)/T_{o}(R) = exp[βΔF(R)] the freeenergy difference between the two states ΔF = F_{c}  F_{o} can be calculated. ΔF increases linearly with the pore radius (Fig. 3a Inset), βΔF(R) = a_{0} + a_{1}R with and . Together with K(R), the gating behavior of the pore is characterized (31). In this sense, the MD calculations have related the input structure to a “physiological” property of the system. (Note, however, that the time scales of ion channel gating and the oscillations observed here differ by 5 orders of magnitude.)
FreeEnergy Density. The Helmholtz freeenergy density f(T, n; R) displays one or two minima (Fig. 4a): one for the empty pore (n = 0) and one in the vicinity of the bulk density. The 0.45nm pore is close to a transition point in the freeenergy landscape: the minimum for the filled pore is very shallow and disappears at smaller radii (R = 0.4 and 0.35 nm). For very large and very small radii, only one thermodynamic stable state exists: liquid or vapor. For intermediate radii, a metastable state appears. Near R = 0.55 nm both states are almost equally probable although they do not coexist spatially because the pore is finite and small. In infinite pores spatially alternating domains of equal length would be expected (35) and were actually observed in MD simulations (36). The oscillating states in short pores, on the other hand, alternate temporally, thus displaying a kind of “timeaveraged” coexistence. For higher densities n/n_{bulk} > 1 the curves start to resemble parabolas, similar to a parabolic f(T, n) seen for cylindrical volumes (data not shown) and spherical cavities (34) in bulk water.
The chemical potential (Fig. 4b) shows the transition from the stable vapor state, μ(T, n) > 0, through the twostate regime to the stable liquid state, μ(T, n) < 0. The features of μ(T, n) indicate that the condensation (and evaporation) processes occur in an avalanchelike fashion: Let the density in the pore be at the transition state, the left zero of μ. If the density is perturbed to increase slightly then μ becomes negative. Every additional molecule added to the pore decreases the free energy further by an amount μ whereas the increase in density lowers the chemical potential even more. This leads to the avalanche of condensation. It only stops when the stable state, the right zero of μ, is reached. Now a further addition of molecules to the pore would actually increase the free energy and drive the system back into the stable state. Similarly, a perturbation that decreases the density in the transition state leads to accelerated evaporation.
From the probability distribution p(N) the freeenergy difference between closed and open state ΔF(R) is calculated, and , consistent with the estimate from the kinetics. ΔF(R) (Fig. 3a Inset) shows the transfer of stability from the vapor state for small R to the liquid state for large R. The coexistence regime is at ΔF(R_{c} = 0.57 nm) = 0.
Dynamics. MD simulations not only allow us to investigate the thermodynamic properties of the system but also the dynamical behavior of individual molecules. A few selected water molecules are depicted in Fig. 1c shortly before the pore empties. They show a diverse range of behaviors and no singlefile like motion of molecules is visible in the liquid state. On evaporation (and condensation) the state changes within ≈30 ps.
The mean permeation time 〈τ_{p}〉 in Table 1 increases with the pore radius, i.e., water molecules permeate narrow hydrophobic pores faster than they diffuse the corresponding distance in bulk water (the reference value). This is consistent with higher diffusion coefficients D_{z} in the narrow pores (up to almost three times the bulk value). The diffusion coefficient perpendicular to the pore axis, D_{xy}, drops to approximately half the bulk value. Marti and Gordillo (37) also observe increased diffusion in simulations on water in carbon nanotubes (D_{z} ≤ 1.65 D_{bulk}) and a corresponding decrease in D_{xy}. Experimental studies on water transport through desformyl gramicidin A (13) can be interpreted in terms of a D_{z} of five times the bulk value. Histograms (data not shown) for τ_{p} show that there is a considerable population of “fast” water molecules (e.g., the black and the dark gray one in Fig. 1c) with τ_{p} between 2 and 10 ps, which is not seen in bulk water. The acceleration of water molecules in the pore can be understood as an effect of onedimensional confinement. The random 3D motion is directed along the pore axis and the particle advances in this direction preferentially. The effect increases with decreasing radius, i.e., increasing confinement. The average equilibrium current density j_{0} follows the trend of the openness closely but more detailed timeresolved analysis shows water translocation to occur in bursts in all pores. In narrow pores, bursts occurring during the “closed” state contribute up to 77% of the total flux (data not shown). For singlefile pores, simulations (14, 22) and theory (38) also point toward concerted motions as the predominant form of transport.
Capillary Condensation. The behavior as described so far bears the hallmarks of capillary condensation and evaporation (19, 39, 40) although it is most often associated with physical systems that are macroscopically extended in at least one dimension, such as slits or long pores. Capillary condensation can be discussed in terms of the Kelvin equation (18), [1] which describes how vapor at pressure p relative to its bulksaturated pressure p_{0} coexists in equilibrium with its liquid. Liquid and vapor are divided by an interfacial meniscus of curvature 1/r (r > 0 if the surface is convex); γ_{lv} is the surface tension between liquid and vapor and ν_{l} the molecular volume of the liquid. Although the Kelvin equation is not expected to be quantitative in systems of dimensions of only a few molecular diameters it is still useful for obtaining a qualitative picture. Curvature 1/r and contact angle θ in a cylindrical pore of radius R are related by R = r cos θ. With Young's equation, γ_{wv} = γ_{wl} + γ_{lv} cos θ, Eq. 1 becomes [2] independent of the interface. For our system, the surface tension between liquid water and the wall, γ_{wl} > 0, and between vapor and the wall, γ_{wv} > 0, are fixed quantities. The hydrophobicity of the wall implies γ_{wv} < γ_{wl}, i.e., the wall is preferentially in contact with vapor; ν_{l} can be considered constant. Hence, for a given pore of radius R there exists one vapor pressure p(R) > p_{0} at which vapor and liquid can exist in equilibrium. Water only condenses in the pore if the actual vapor pressure exceeds p(R). Otherwise, only vapor will exist in the pore. The effect is strongest for very narrow pores. Hence a higher pressure is required to overcome the surface contributions, which stabilize the vapor phase in narrow pores. The pressure fluctuates locally in the liquid bulk “reservoir.” These fluctuations can provide an increase in pressure above the saturation pressure in the pore and thus drive oscillations between vapor and liquid.
Comparison with Experiments, Simulations, and a Theoretical Model. Experiments on aquaporins (9, 10) and gramicidin A (12, 13) yield osmotic permeability coefficients of water, p_{f}, of the order of 10^{12} to 10^{14} cm^{3}·s^{1} (Table 3, which is published as supporting information on the PNAS web site, www.pnas.org). We calculate p_{f} = 1/2Φ_{0}ν_{l} from the equilibrium flux of our MD simulations (14) and find that narrow (R = 0.35 and 0.4 nm), predominantly closed pores have p_{f} ≈ 5 × 10^{14} cm^{3}·s^{1}, that is, the same magnitude as Aqp1, AqpZ, and gramicidin A. As these pores are longer (≈2 nm) and narrower (R < 0.2 nm) than our model pores, strategically placed hydrophilic groups (15) seem to be needed to stabilize the liquid state and facilitate water transport in these cases.
Recently Giaya and Thompson (41) presented an analytical meanfield model for water in infinite cylindrical hydrophobic micropores. They predict the existence of a critical radius R_{c} for the transition from a thermodynamically stable water vapor phase to a liquid phase. The crucial parameter that R_{c} depends on is the water–wall interaction. We choose the effective fluid–wall interaction ε_{eff} = ρ_{w} ε_{fw}, the product of the density of wall atoms with the well depth of the fluid–wall interaction potential, as a parameter to compare different simulations because this seems to be the major component in the analytical fluid–wall interaction. As shown in Table 2, compared with carbon nanotube simulations our pore has a very small ε_{eff} and thus can be considered extremely hydrophobic. This explains why Hummer et al. (22) observe permanently waterfilled nanotubes with a radius of only 0.24 nm although their bare fluid–wall interaction potential is weaker than in our model. The much higher density of wall atoms in the nanotube, however, more than mitigates this. Once they lower their ε_{eff} to double of our value, they also observe strong evaporation. This suggests that the close packing of wall atoms within a nanotube may result in behavior not seen in biological pores. The mean field model agrees qualitatively with the simulations as it also shows a sharp transition and high sensitivity to ε_{eff}.
Conclusions
We have described oscillations between vapor and liquid states in short (L_{P} = 0.8 nm), hydrophobic pores of varying radius (0.35 nm ≤ R ≤ 1 nm). Qualitatively, this behavior is explained as capillary evaporation, driven by pressure/density fluctuations in the water “reservoir” outside the pore. Similar behavior is found in simulations by different authors with different water models [SPC (this work), SPC/E, TIP3P (data not shown), TIP3P (22), TIP4P (25)] in different nanopores [atomistic flexible models (this work), carbon nanotubes (22), spherical cavities (24) and smooth cylinders (25, 42)].
We presented a radically simplified model for a nanopore that is perhaps more hydrophobic than in real proteins [although we note the existence of a hydrophobic pore in the MscS channel (6)]. From comparison with experimental data on permeability we conclude that strategically placed hydrophilic groups are essential for the functioning of protein pores. The comparatively high permeability of our closed pores suggests pulsed water transport as one possible mechanism in biological water pores. Local hydrophobic environments in pores may promote pulsatory collective transport and hence rapid water and solute translocation.
Our results indicate intrinsically collective dynamic behavior not seen on simulation time scales currently considered sufficient in biophysical simulations. These phase oscillations in simple pores, a manifestation of the hydrophobic effect, require >50 ns of trajectory data to yield a coherent picture over a freeenergy range of 6 k_{B}T. We thus cannot safely assume that the behavior of water within complex biological pores may be determined by extrapolation from our knowledge of the bulk state or short simulations alone.
Acknowledgments
We thank all of our colleagues for their interest in this work, especially Joanne Bright, José FaraldoGómez, Andrew Horsfield, and Richard Law. This work was funded by The Wellcome Trust.
Footnotes

↵* To whom correspondence should be addressed. Email: mark{at}biop.ox.ac.uk.

This paper was submitted directly (Track II) to the PNAS office.

Abbreviations: MD, molecular dynamics; SPC, simple point charge; PMF, potential of mean force.
 Received November 10, 2002.
 Accepted April 4, 2003.
 Copyright © 2003, The National Academy of Sciences
References
 ↵
Doyle, D. A., MoraisCabral, J., Pfützner, R. A., Kuo, A., Gulbis, J. M., Cohen, S. L., Chait, B. T. & MacKinnon, R. (1998) Science 280, 6977.pmid:9525859

Chang, G., Spencer, R. H., Lee, A. T., Barclay, M. T. & Rees, D. C. (1998) Science 282, 22202226.pmid:9856938

Fu, D., Libson, A., Miercke, L. J., Weitzman, C., Nollert, P., Krucinski, J. & Stroud, R. M. (2000) Science 290, 481486.pmid:11039922
 ↵
Bass, R. B., Strop, P., Barclay, M. & Rees, D. C. (2002) Science 298, 15821587.pmid:12446901
 ↵
Hille, B. (2001) Ion Channels of Excitable Membranes (Sinauer, Sunderland, MA), 3rd Ed.
 ↵
 ↵
 ↵
Pohl, P., Saparov, S. M., Borgnia, M. J. & Agre, P. (2001) Proc. Natl. Acad. Sci. USA 98, 96249629.pmid:11493683
 ↵
 ↵
 ↵
 ↵
 ↵
Tajkhorshid, E., Nollert, P., Jensen, M. O., Miercke, L. J., O'Connell, J., Stroud, R. M. & Schulten, K. (2002) Science 296, 525530.pmid:11964478
 ↵
de Groot, B. L. & Grubmüller, H. (2001) Science 294, 23532357.pmid:11743202
 ↵
Finkelstein, A. (1987) Water Movement Through Lipid Bilayers, Pores, and Plasma Membranes: Theory and Reality (Wiley, New York).
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
Unwin, N. (2000) Philos. Trans. R. Soc. London B 355, 18131829.pmid:11205343
 ↵
Lindahl, E., Hess, B. & van der Spoel, D. (2001) J. Mol. Mod. 7, 306317.
 ↵
 ↵
 ↵
 ↵
Sakmann, B. & Neher, E., eds. (1983) SingleChannel Recordings (Plenum, New York).
 ↵
 ↵
 ↵
Hummer, G., Garde, S., García, A. E., Pohorille, A. & Pratt, L. R. (1996) Proc. Natl. Acad. Sci. USA 93, 89518955.pmid:11607700
 ↵
 ↵
 ↵
Martí, J. & Gordillo, M. C. (2001) Phys. Rev. E Phys. Plasmas Fluids Relat. Interdiscip. Top. 64, 02150410215046.
 ↵
Berezhkovskii, A. & Hummer, G. (2002) Phys. Rev. Lett. 89, 06540310654034.
 ↵
Rowlinson, J. S. & Widom, B. (1982) Molecular Theory of Capillarity (Clarendon, Oxford).
 ↵
 ↵
 ↵