Water in hydrophobic nanopores

The behaviour of water in highly simplified model pores was investigated in molecular dynamics computer simulations. An unexpected two-state behaviour was found. Water either fills the pore lumen at about bulk density and forms a liquid phase, or the pore is almost void of water (vapour). We discuss these findings and their relevance in a biological context.

This page deals exclusively with water but further research on ions in nanopores is described in Phys. Biol. 1 (2004) and J. Am. Chem. Soc. 126 (2004).

Table of contents top site menu

Water

The hydrophobic gating hypothesis

Model system and Molecular Dynamics

Water oscillations and how to quantify them

A hydrophobic gating mechanism for nanopores

Liquid-vapour oscillations: The longer time-scale behaviour

Free energy landscape of the pore-water system

Transport properties

Conclusions (short)

Water

Water is remarkable stuff: It is the most abundant liquid on earth and its physical properties are quite different from many other liquids.

• Ice (the solid phase of water) at 273.15 K has a higher density than water at 273.15 K; most other materials become denser on solidification. The density of water reaches a maximum at 277.14 K (ca. 4°C). This so-called density anomaly is indicated by ice floating on water and lakes freezing from top to bottom).
• Ice can be melted by applying pressure; ordinary liquids can be solidified by pressure (and their solids become even more rigid and not fluid).
• H₂O has a higher melting and boiling point than comparable substances and it has a high heat capacity [75.2 J(mol K)^-1 (liquid), 37.6 J(mol K)^-1 (solid)].
• Its dielectric constant is high (around 80), so it is a very good solvent for ions and its surface tension is high (72 x 10^-3 J m^-2).
• Protons are diffusing very fast through water due to the Grotthus mechanism of hydrogen bond flipping instead of physical diffusion of the proton itself.

All these properties can be related in one way or another to the extensive capability of water to form strong hydrogen bond networks.

Life as we know it is water-based. All organisms rely on water to be available, e.g. as a solvent for ions, as a participant in enzymatic reactions or as a means to enhance the stability of a cell by using pressurised, water-filled cavities. In short, it is essential for biological systems (i.e. cells).

In the context of ion channels water is important as the solvent for the ions. Water molecules surround the ions in a tight hydration shell, which screens the strong electric field of the ions and prevents them from forming salt crystals, thus keeping the ions mobile. The energetic cost to remove the hydration shell is large (of the order of 100 k_BT); this cost would be incurred if the ion was to be transferred into a low-dielectric (i.e. non-polar or hydrophobic) environment. When ions are required to cross the cell's lipid membrane through ion channels they must either take their hydration shell with them or functional groups of the protein must substitute for the hydration shell in part or in full. These functional groups are invariably polar to mimic the electrostatic interaction between water and ion. Some channels like KcsA are known to employ this effect for fine tuned selectivity. In addition, evidence is accumulating that gates in ion channels might also rely on the high cost of ionic dehydration to prevent ions from permeating.

The hydrophobic gating hypothesis

We have the atomic structures of a handful of channels. Amongst these are the bacterial potassium channel KcsA (pdb:1K4C), the mechanosensitive channels of large (1MSL) and small (1MXM) conductance, and the nicotinic acetylcholine receptor nAChR (1OED). These are examples of gated channels. A unifying motive emerges from these structures. The putative gate, one of the constriction sites of the pore, is lined by a ring of hydrophobic residues (see boxed regions in the figure below).

Putative hydrophobic gates in three channels. The proteins are drawn from the Protein Databank coordinates in cartoon representation with alpha-helices represented by cylinders and beta-sheets by arrows, using VMD. The surface of the pore was created with the programme HOLE. Blue indicates a radius which allows more than one water molecule to fit into the diameter, green signifies single-file geometry, and a red pore has to be flexible to accommodate at least one water molecule of diameter ca. 0.28 nm (note that typical bare ionic radii start at 0.095 nm and hence diameters of 0.19 nm).

Of course, the simplest way to close passage to an ion is to physically occlude the pathway. This mechanism is most likely employed in some channels. However, it might require some fairly large-scale domain motions. It is conceivable that nature uses an alternative to complete occlusion. If the pore in the closed state were only narrow enough not to let an ion pass with its complete hydration shell then the cost of (partial) desolvation would be prohibitively high and the pore would appear closed to the ion for all practical purposes. This would require the gate to be lined with a hydrophobic or non-polar surface, so that hydration shell substitution cannot take place.

nAChR is a good example for this effect. Its structure was determined in the open and in the closed state (see N. Unwin (2000) for an overview) and shows the changes that occur on gating. In the open state the receptor presents a pore of at least 0.65 nm radius to ions, which can easily pass with their hydration shell; in addition, the pore is lined by polar residues. The gate in the closed state is a pore of a radius 0.3 nm < R < 0.35 nm---it is not physically occluded and would allow cations to pass (r_cation=0.06...0.17 nm). Unwin et al (1999) suggested that in nAChR the closed gate is formed by a hydrophobic girdle of highly conserved Val and Leu residues which would present a high dehydration barrier to an ion.

In order to test the hypothesis that a hydrophobic pore does indeed present a substantial barrier to ion permeation we carried out molecular dynamics computer simulations on a extremely simplified pore model. First we focused on water without ions because we wanted to understand the behaviour of the "simplest" system before progressing to more complicated ones (ion+water). In addition, we assume that we can use water permeation as a "proxy" for ion permeation through hydrophobic pores because if the environment is unfavourable for water then it is certainly unfavourable for ions, too.

Model system and Molecular Dynamics

A model pore with the typical dimensions of the nAChR gating region (an hour glass shape with pore length L=0.8 nm and radii 0.35 nm ≤ R ≤ 1.0 nm; the dimensions of the mouth regions were kept fixed at L_M=0.8 nm R_M=1.0 nm) was created from pseudo atoms with the characteristics of methane molecules (approximating the surface created by the methyl groups of Leu, Val or Ile). Water was modelled according to the simple point charge scheme (SPC).

Molecular dynamics (MD) simulations were performed with GROMACS at constant temperature T=300 K and pressure P=1 bar. The simulation box was treated within periodic boundary conditions; hence electrostatic interactions were computed with a particle-mesh-Ewald (PME) algorithm. The systems typically contained 4000 SPC water molecules and 700 methane pseudo atoms. Calculations proceeded at a speed of approximately 1.5 ns/day on dual processor workstations. The total accumulated trajectory time was about 0.5 μs.

Water oscillations and how to quantify them

The MD simulations show a very dynamical behaviour of water indeed. The density in pores of radius around 0.5 nm fluctuates between a filled pore at approximately bulk density and a practically empty pore.

This dynamic behaviour is clearly seen in a movie of water in a hydrophobic nanopore (mpeg) with a radius of 0.55 nm.

Oscillations between filled and empty states and definition of the porestate function. Top: time course of the density in a hydrophobic pore (R=0.6  nm). The density is measured in units of the bulk water density 0.996 g cm^-3. Middle: porestate indicator function ω(t). Bottom: z-coordinates of water molecules over time. The porestate indicator function equals one if liquid water fills the pore but equals 0 if no water or not enough water to form a continuous chain of water molecules is in the pore.

The water-pore system oscillates between two clearly distinguishable states, which will be called open (full) and closed (empty) because we assume that an ion could conceivably permeate a pore that can contain liquid water but not the empty one.

In order to quantify the effect of different pore geometries and pore linings on the water behaviour we define a functional measure. A convenient number is the probability that the pore is in the open state, i.e. the expectation value of the porestate indicator function ω(t). We call this the openness and calculate it as the ratio <ω> = T_open/T_sim. If a pore is water-filled during the whole simulation time it will have an openness of 1; if it is open half of the time and closed during the other half it will have openness 0.5, and a permanently empty pore will have <ω>=0. The openness does not contain any information how often the pore switches states, it is only related to the equilibrium constant K of the two state system OPEN↔CLOSED by K = T_closed/T_open = <ω>^-1 - 1. The switching frequency is determined by the "activation energies" for the filling and emptying transitions and only determine how quickly equilibrium will be attained.

A hydrophobic gating mechanism for nanopores

Short MD simulations (simulation time T_sim < 6 ns) already tentatively confirmed the proposed hydrophobic gating mechanism (Beckstein et al, J. Phys. Chem. B 105 (2001), 12902--12905).

Here the openness stands for the probability that an ion could permeate the pore. As shown in the graph at the right, it is close to 0 at small radii (0.35 nm to 0.5 nm). Then it increases sharply to attain its full value of 1 for R > 0.6 nm. The graph demonstrates that

• a hydrophobic pore of the dimensions of the closed nAChR gate is practically impermeable to liquid water although at least three water molecules could fit into the pore diameter, and that
• a pore of the dimensions of the open nAChR gate contains water at bulk density even though it has a hydrophobic pore lining.

When the pore lining is made more hydrophilic by adding small dipoles parallel to the pore axis (the dipoles have the same dipole moment as the dipole in the peptide bond) a closed hydrophobic pore can be opened. Adding two dipoles (blue square) to the closed R=0.4 nm or 0.5 nm pore switches it to completely open. (One dipole is not strong enough and the results are inconsistent due to the short simulation time and perhaps some geometrical stabilisation effects in the 0.35 nm pore.) This again resembles the nAChR open gate: by replacing hydrophobic residues by polar ones the equilibrium between open and closed states is firmly shifted towards the open state.

In conclusion, the MD simulations indicated that the hydrophobic gating hypothesis might be true and protein channels could indeed employ the hydrophobic barrier presented by a cylindrical pore with a non-polar surface. By combining (1) pore radius reduction with (2) hydrophobic (hydrophilic) pore lining in the closed (open) state a protein can obtain a sufficient gating effect without having to resort to extreme conformational changes.

Liquid-vapour oscillations: The longer time-scale behaviour

The previous calculations lacked a better resolution of the hydrophobic gating effect which could be seen in fairly large error bars near the transition radius of R_c=0.55 nm; only a few switching events could be recorded resulting in poor statistics. Thus, the simulations were extended from 6 ns to beyond 50 ns. The results of these simulations are published in Beckstein and Sansom, PNAS 100 (2003), 7063-7068 doi: 10.1073/pnas.113684410.

The density in the pore oscillates between a liquid state with approximately 0.81 bulk density and a vapour state with density close to 0. The time scale for these oscillations is greater than nano seconds and it is now also possible to get good estimates for the full range of radii.

Oscillating water density in model pores of increasing pore radius R. The water density n(t) (in units of the bulk water density) over the simulation time shows strong fluctuations on a greater than ns time scale (50 ps moving average smoothing). Two distinctive states are visible: open at approximately n_bulk (liquid water), and closed with very few or no water in the pore (water vapour).

Openness. Wide pores are permanently water filled whereas narrow ones cannot sustain liquid water. The broken sigmoidal curve is derived from the free energy difference between open and closed state as determined by the equilibrium probability distribution.

The openness curve does not differ significantly from the one from the short simulations but because of the much more thorough sampling of states it is now possible to obtain meaningful local densities (shown in the next graph) and collect statistics to obtain free energies for the open and closed pore states.

Density of water in hydrophobic nanopores. Molecular dynamics simulations of water in hydrophobic model pores show strong one-dimensional confinement effects on water behaviour. The water density exhibits layering for liquid water (red/orange). Pores of radii smaller than 0.45 nm predominantly contain water vapour (dark blue) although they are still much wider than single water molecules. The density is measured relative to the density of bulk water (1.00 corresponds to n_bulk=53.7 mol L^-1, SPC water at 300 K and 1 bar). Graphic created with XFarbe

The density plots show that close to a hydrophobic surface water is strongly layered and at least three layers (spanning roughly 1 nm) are visible. Interestingly, the density is strongly enhanced near protruding edges. In larger pore, the surface layering simply "wraps" into the pore. Because these pictures are averages over the whole simulation time pores that are closed an appreciable proportion of time (R<0.6 nm) show a decreased average density in the pore. Nevertheless, the same layering structure is visible down to below 0.55 nm. For even smaller pores, no water structure ever appears, and the smallest pores appear almost void of water.

Free energy landscape of the pore-water system

A molecular dynamics trajectory samples configurations or states of the system according to the Boltzmann distribution. We are interested in the occupancy N of the pore, i.e. the number of water molecules in the pore (from the occupancy the density in the pore can be calculated with the pore volume V as n = <N>/V). Thus we label states by N and consider all states equivalent with same N. The probability P(N) that the pore contains exactly N water molecules can be easily obtained from the MD trajectory by simply counting states. If the simulation time is long enough then all states of interest of the system are sampled sufficiently frequently to yield meaningful statistics. The graphic on the left shows how the pore occupancy distribution reflects the two-state behaviour. One peak of high probability corresponds to an empty pore (N=0), the other one to a filled pore. These two states are separated by a "transition state", which is taken to be the state of lowest probability between the two maxima. (Note that this definition is only useful in the cases when there are actually two maxima observed. Otherwise, the definition in terms of the porestate indicator function seems to capture the essentials quite well, and it coincides with the definition in the two-state (or two-maxima) regime)

The free energy of a given pore state is then simply

F (T,V,N) = -k_BT ln P(N) + C,

f (T,n) = V^-1 F (T,V^-1 V, V^-1 N).

This has the advantage that the free energy density landscapes of different pores can be more easily compared than the absolute free energies themselves because the trivial volume (or radius) dependency is removed from the densities.

Top: Free energy density f (T,n) at constant T=300 K.
Bottom: Chemical potential μ(T,n). n is the water density in the pore, normalised to n_bulk=53.7 mol L^-1. f is given in units of k_BT and the inverse of the liquid molecular volume of bulk water (v_l^-1=n_bulk). Two minima correspond to the observed two-state behaviour. The vapour state becomes metastable with increasing radius and for R>0.55 nm the liquid state is globally stable. f (T,n;R=1.0 nm) is drawn with an arbitrary offset.

The Helmholtz free energy density f (T,n;R) displays one or two minima: one for the empty pore (n=0) and one in the vicinity of the bulk density. The 0.45 nm pore is close to a transition point in the free energy landscape: the minimum for the filled pore is very shallow and disappears at smaller radii (R=0.4 nm and 0.35 nm). For very large and very small radii, only one thermodynamic stable state exists: liquid or vapour. 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. The oscillating states in short pores, on the other hand, alternate temporally, thus displaying a kind of "time-averaged" coexistence.

From the free energy density the chemical potential μ can be calculated,

μ(T,n) = ∂f (T,n) ⁄ ∂n

μ is the free energy to add or remove a water molecule to a pore when there is already water at density n in there.

The chemical potential shows the transition from the stable vapour state, μ(T,n)>0, through the two-state regime to the stable liquid state, μ(T,n)<0. The features of μ(T,n) indicate that the condensation (=filling) and evaporation (=emptying) processes occur in an avalanche-like 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 μ while 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.

The graph of the chemical potential also shows clearly that for the very narrow pores (R=0.4 nm and R=0.35 nm) no stable liquid state exists. The R=0.45 nm pore is very close to a transition point in the free energy landscape which divides the pores that only have a stable vapour state from the ones which can have a vapour and a liquid stable state. There is probably another transition radius 0.7 nm<R<∞, beyond which there is no stable vapour state and only a stable liquid state remains.

Transport properties

MD simulations offer the unique opportunity to watch the behaviour of single molecules. This is primarily of value to develop a "feeling" for the system of interest; hard numbers are won by thermodynamic averaging over many particles and/or long time scales. Nevertheless it is instructive to look at the permeant water molecules in a pore as depicted in the graph on the left below.

Permeant water molecules in a R=0.55 nm pore as it switches from the open to the closed state. z-coordinates of the water oxygen atoms are drawn every 2 ps. The mouth and pore region are indicated by horizontal broken and solid lines. Five trajectories are shown explicitly. The white water molecule permeates the pore within 54 ps whereas the black one only requires about 10 ps.

Transport properties of water in hydrophobic pores. <ω> is the openness. The mean permeation time τ_p is measured relative to the bulk value, τ_p,bulk = 29.9±0.1 ps. The equilibrium current density j₀ is the total number of permeant water molecules per unit time and unit area ( j_0,bulk = 320±3 ns^-1 nm^-2). The diffusion coefficient along the pore axis D_z is normalised to the bulk value of SPC water at 300 K and 1 bar ( D_bulk = 4.34±0.01 nm² ns^-1).

The most obvious observation is the great variety of behaviours. Some water molecules permeate the pore rather slowly (like the white one; others almost jump through the pore (e.g. the black or the orange one). Waters can pass each other and can even permeate the pore in opposite directions at the same time. This indicates that these pores are to be considered as multi-pass pore (as opposed to very narrow single-pass pores where water transport occurs in single-file). The pore empties fairly rapidly; on average filling and emptying takes about 30 ps, rather independent of the pore radius.

Analysis of the water in the pore shows that the 1D confinement exerted by the pore on the water changes its transport properties quite profoundly. This is most clearly brought out by normalising the observed quantities by the corresponding bulk value. As shown in the right graph, the apparent diffusion coefficient in z direction (along the pore axis) increases to almost three times its bulk value for very narrow pores. Similarly, water molecules only take about half the time to permeate narrow hydrophobic pores when compared with the average permeation time in wide (or infinite = bulk) pores. The number of molecules per ns and area, however, drops in line with the openness. This indicates that the major contribution to water transport through these pores occurs in the liquid phase. However, this is not true for the very narrow pores (0.35 nm and 0.4 nm). Additional analysis shows that 77% of permeant water molecules cross the 0.35 nm pore while the state is considered close (and ca. 50% for the 0.4 nm pore). These transport events occur in bursts when a few water molecules travel in small clusters or chains through the otherwise empty pore lumen. When the equilibrium flux rates Φ₀ of the narrow pores, which are of the order of 1...3.2 molecules ns^-1, are compared to experimental flux rates of biological water transporters, one finds that they are of similar magnitude (see Table S1 on Osmotic permeability coefficients and flux [ pdf ])

(Click on the table for a larger version including references or view it in pdf format. Also published as a supplementary table to the article in PNAS.)

That means, that already the "closed" pores mimic the transport capabilities of real proteins. Hence, for water transport the description "closed" is a misnomer for these pore states. It rather suggests that in aquaporins transport does not occur by ordinary diffusion but also by a burst-like mechanism. A string of hydrogen-bonded water molecules may slide rather effortlessly through the "greasy" hydrophobic pore.

Conclusions (short)

A hydrophobic environment can act as a additional barrier to water and ion permeation as in the hydrophobic gating mechanism. 1D hydrophobic confinement also influences the transport of water profoundly. It can accelerate water transport and (probably) induce cooperative transport effects which can result in a burst-like transport mechanism.

Thermodynamic analysis of the pore-water system shows how two-state behaviour can emerge from seemingly simple components (a hole in a slab, surrounded by water). As explained in Beckstein et al. PNAS (2003), the oscillations are a consequence of pressure/density fluctuations in the water bath outside the pore. A narrow hydrophobic pore is more difficult to fill with water than a large one. As the pressure fluctuations can only correspond to changes in free energy of about k_BT a narrow pore can be less often forced to become filled than a larger one. This explains qualitatively the dependence of the density fluctuations on the pore radius.

The hydrogen bonding capabilities of water are critical for this behaviour. A simple non-associative liquid simply fills the pore and does not exhibit these seemingly peculiar oscillations between empty and filled pore states. Water, after all, is somehow special....

SBCB   |  O.B.  |  Research  |  Publications  |  Non-work