Article Information

• PDF (536 kb)
• Full Text with Large Figures
• Cited by in Scopus (28)

PubMed

• Articles by Christine Peter
• Articles by Gerhard Hummer

Related Articles

• Nanobiocatalysis and its potential applications

  Nanobiocatalysis and its potential applications Trends in Biotechnology, Volume 26, Issue 11, 1 November 2008, Pages 639-646 Jungbae Kim, Jay W. Grate and Ping Wang Abstract Nanobiocatalysis, in which enzymes are incorporated into nanostructured materials, has emerged as a rapidly growing area. Nanostructures, including nanoporous media, nanofibers, carbon nanotubes and nanoparticles, have manifested great efficiency in the manipulation of the nanoscale environment of the enzyme and thus promise exciting advances in many areas of enzyme technology. This review will describe these recent developments in nanobiocatalysis and their potential applications in various fields, such as trypsin digestion in proteomic analysis, antifouling, and biofuel cells. Abstract | Full Text | PDF (752 kb)

• Nanopore sequencing technology: nanopore preparations

  Nanopore sequencing technology: nanopore preparations Trends in Biotechnology, Volume 25, Issue 4, 1 April 2007, Pages 174-181 Minsoung Rhee and Mark A. Burns Abstract For the past decade, nanometer-scale pores have been developed as a powerful technique for sensing biological macromolecules. Various potential applications using these nanopores have been reported at the proof-of-principle stage, with the eventual aim of using them as an alternative to DNA sequencing. Currently, there have been two general approaches to prepare nanopores for nucleic acid analysis: organic nanopores, such as α-hemolysin pores, are commonly used for DNA analysis, whereas synthetic solid-state nanopores have also been developed using various conventional and non-conventional fabrication techniques. In particular, synthetic nanopores with pore sizes smaller than the α-hemolysin pores have been prepared, primarily by electron-beam-assisted techniques: these are more robust and have better dimensional adjustability. This review will examine current methods of nanopore preparation, ranging from organic pore preparations to recent developments in synthetic nanopore fabrications. Abstract | Full Text | PDF (1679 kb)

• Simulations of Electrophoretic RNA Transport Through Transme...

  Simulations of Electrophoretic RNA Transport Through Transmembrane Carbon Nanotubes Biophysical Journal, Volume 94, Issue 7, 1 April 2008, Pages 2546-2557 Urs Zimmerli and Petros Koumoutsakos Abstract The study of interactions between carbon nanotubes and cellular components, such as membranes and biomolecules, is fundamental for the rational design of nanodevices interfacing with biological systems. In this work, we use molecular dynamics simulations to study the electrophoretic transport of RNA through carbon nanotubes embedded in membranes. Decorated and naked carbon nanotubes are inserted into a dodecane membrane and a dimyristoylphosphatidylcholine lipid bilayer, and the system is subjected to electrostatic potential differences. The transport properties of this artificial pore are determined by the structural modifications of the membrane in the vicinity of the nanotube openings and they are quantified by the nonuniform electrostatic potential maps at the entrance and inside the nanotube. The pore is used to transport electrophoretically a short RNA segment and we find that the speed of translocation exhibits an exponential dependence on the applied potential differences. The RNA is transported while undergoing a repeated stacking and unstacking process, affected by steric interactions with the membrane headgroups and by hydrophobic interaction with the walls of the nanotube. The RNA is structurally reorganized inside the nanotube, with its backbone solvated by water molecules near the axis of the tube and its bases aligned with the nanotube walls. Upon exiting the pore, the RNA interacts with the membrane headgroups and remains attached to the dodecane membrane while it is expelled into the solvent in the case of the lipid bilayer. The results of the simulations detail processes of molecular transport into cellular compartments through manufactured nanopores and they are discussed in the context of applications in biotechnology and nanomedicine. Abstract | Full Text | PDF (1934 kb)

• …more

Copyright © 2005 The Biophysical Society. All rights reserved.
Biophysical Journal, Volume 89, Issue 4, 2222-2234, 1 October 2005

doi:10.1529/biophysj.105.065946

Biophysical Theory and Modeling

Ion Transport through Membrane-Spanning Nanopores Studied by Molecular Dynamics Simulations and Continuum Electrostatics Calculations

Christine Peter¹ and Gerhard Hummer^,  

Laboratory of Chemical Physics, National Institute of Diabetes and Digestive and Kidney Diseases, National Institutes of Health, Bethesda, MD 20892-0520

Address reprint requests to Gerhard Hummer, National Institute of Diabetes and Digestive and Kidney Diseases, National Institutes of Health, Bethesda, MD 20892-0520.

¹ Christine Peter’s present address is Max-Planck Institute for Polymer Research, Ackermannweg 10, D-55128 Mainz, Germany.

Ion transport across membranes, mediated by ion-channel proteins, is central to many biological processes 1,2, such as neural signal transmission and electrical excitability of muscle. Recent progress in the expression, purification, and crystallization of membrane proteins has led to a rapidly growing number of high-resolution structures 3,4, including those of several ion channels 5,6,7,8. Channel structures, together with detailed biochemical and biophysical characterizations using, e.g., amino-acid mutation studies and single-channel conductance measurements, provide key insights into the energetics, selectivity, and gating of ion channels. The availability of atomic structures has also enabled increasingly detailed computational studies that relate ion-channel structure and dynamics to functional mechanisms 9,10,11,12,13,14,15,16,17,18,19,20,21.

Understanding ion selectivity has been a major focus of both structural and computational studies, and enormous progress has been made in particular for potassium channels, which are highly permeable for K⁺ ions but not Na⁺ ions. Structural studies of potassium channels have also provided new insight into ion-channel gating, relating ion permeability and impermeability to differences in the channel structures between the “open” and “closed” states, respectively. The first high-resolution x-ray structures showed the potassium channel in a closed state, but recently the structure of the MthK potassium channel in the open state was reported 22. Remarkably, the open pore of the MthK channel is lined predominantly by nonpolar amino acids. One of the questions we address in this article is a functional role of the nonpolar pore region of potassium channels.

Narrow and relatively hydrophobic pore regions are a repeatedly observed feature of biological channels facilitating the transport of polar species such as water, protons, and ions. They are present in the potassium channels, e.g., KcsA, MthK, or the Kir family 7,22,23. The mechanosensitive channels (e.g., MscS 24), and the nicotinic acetylcholine receptor nAChr 25 have similar hydrophobic regions, which may function as gating regions that control the ion permeability 26,27. Narrow, mostly hydrophobic cavities are also involved in the regulation of water or proton transport in several biological systems, such as aquaporins 28,29, cytochrome c oxidase 30, and bacteriorhodopsin 31.

The presence of nonpolar pores as conduits for ions seems at first sight rather surprising for electrostatic reasons. From a simple dielectric model of ion permeation through a hydrophobic pore, one might expect a substantial barrier for displacing a charged particle from the high dielectric solution into a channel surrounded by a low dielectric medium, uncompensated by strong interactions with polar residues lining the channel wall 32. The high conductivity of potassium channels in the open state thus raises the question whether hydrophobic channels are more permeable for ions than one may think based on this simple dielectric argument. To address this question, we will use a combination of computer simulations and continuum electrostatics calculations.

Biomolecular systems as large and complex as membrane channels are becoming accessible to computer simulations thanks to advances both in the raw computational power and in the simulation methodology 9,10,11,12,13,14,15,16,17,18,33. However, one problem in theoretical ion-channel research is the complexity of the resulting systems of membrane-embedded proteins, making it difficult to identify the principles underlying the functional mechanisms in the wealth of detailed information provided by molecular simulations. The irregularly shaped pores are lined by both hydrophilic and hydrophobic regions of various widths, and single amino-acid residues can play a crucial role at specific sites. Separating the overall free energies for ion transport into individual contributions to quantify the roles of, e.g., “electrostatics”, “hydrogen bonds”, “desolvation”, or “conformational fluctuations” thus proves difficult, with results that can depend on the separation procedure.

Therefore, to answer the general question about a functional relevance of hydrophobic pores for gating we pursue a different strategy. We avoid the complexity of the biological systems by studying simpler model systems designed to capture specific properties. In addition to being more easily interpretable, such systems are also computationally less expensive than the more realistic membrane-protein systems. For ion translocation through channels, this is of particular relevance, because long-range electrostatics plays an important role in the energetics. To avoid surface effects, molecular dynamics simulations of molecular systems are usually performed either with periodic boundary conditions or for solvent droplets, possibly “embedded” in a dielectric continuum. In either case, electrostatic artifacts seem possible because of the long range of the Coulomb interaction, the strong dielectric discontinuity at the membrane, and the weakly screened interactions within the low-dielectric membrane region. By using a comparatively simple system, we can explicitly explore the energetic contributions coming from the simulation boundary conditions. Moreover, using a simple channel system enables us also to perform direct and fully quantitative comparisons of explicit molecular dynamics simulations to continuum electrostatics calculations. Normally, the motions of protein and membrane, and the inhomogeneous dielectric environment permit only qualitative comparisons. Here, we will study a system with a relatively rigid pore and a well-defined dielectric constant of the entire membrane region.

We focus on the basic principles that govern ion transport through narrow hydrophobic pores by simulation of model membranes formed of hexagonally packed carbon nanotubes (CNTs), as displayed in Fig. 1. Central questions of our study are: Which factors determine the energetics of ion translocation through such nonpolar nanopores? How large is the free-energy barrier for ion translocation through a hydrophobic pore with a diameter similar to that of the KcsA channel in the closed-state structure 7? How much does that barrier drop if the pore is widened to a diameter corresponding to that of the open-state structure of the MthK potassium channel 22, and how mobile are ions in wide channels?

Display large version of this figure

Figure 1

Schematic of simulation setup. (A) Snapshot of the channel simulation system with eight copies of the periodic unit cell. Each unit cell contains one open and three capped CNTs, which are periodically replicated in the x, y plane, resulting in a hexagonally packed membrane. In the z direction, the periodic copies of the CNT membrane are separated by a water layer. (B and C) Illustration of the channel and sieve setups, respectively. The content of one unit cell is indicated with gray atoms. (D) Schematic of the simulation system, indicating the translocation of the Na⁺ ion through the open pore in a “channel” system.

Our model system also allows us to explore methodological issues in computational studies of ion translocation 34. By comparing continuum electrostatics calculations to results from molecular dynamics simulations, we will address the question which level of theory/simulation is required to understand various aspects of the translocation process.

The article is organized as follows. We first present different computational approaches to investigate ion transport through membrane channels and describe the simulation methods used in this study. Then, we present our results on the translocation of a Na⁺ ion through model membranes formed by carbon nanotubes. Finally, based on our simulation results, we discuss the role of hydrophobic cavities and the water phase inside them in the process of charge transport through biological membranes.

Theoretical ion-channel research includes a large variety of simulation techniques, ranging from implicit-solvent simulations that rely on continuum electrostatics to computationally much more expensive classical (explicit-solvent) or ab initio molecular dynamics (MD) simulations. Continuum descriptions using Poisson-Nernst-Planck theory as well as explicit descriptions of ions in Brownian dynamics or MD simulations are used to study ion transport 20,35,36,37,38. In (all-atom) MD simulations, the inclusion of solvent greatly increases the computational cost, without even considering electronic degrees of freedom. Typically, large system sizes and long simulations are required to obtain reliable results in simulations of membrane proteins. In continuum models, one faces several well-known limitations 39,40,41,42: 1), bulk properties of water such as density and dielectric constant are applied uniformly to the water phase in the system, even though the properties of water might locally deviate strongly from those of bulk water; 2), the assumption of a uniformly distributed polarizable medium is problematic for cavities that contain few and inhomogeneously distributed water molecules; and 3), surface polarization effects cannot be accounted for in simple continuum models since they are a result of surface ordering of water molecules without an external electric field.

The carbon-nanotube model membranes used in our study result in systems that are small compared to biological membrane channels and thus computationally relatively inexpensive. The small size and relative simplicity of our system allow us to investigate the thermodynamics and kinetics of ion translocation through a hydrophobic pore with explicit-solvent MD simulations. We used an umbrella sampling technique 43,44 to obtain the free-energy for moving an ion from bulk solvent into a nanopore embedded in a membrane. For reference, we also computed the free energy of ion translocation using a dielectric-continuum model. By comparing the free energies obtained from all-atom MD simulations and continuum electrostatics, we identify those aspects of the ion translocation energetics determined specifically by the local properties of solvent molecules. This comparison is supported by the analysis of several structural properties of the ion and the solvent in the MD simulations, among which the local electrostatic potential is of particular interest for charge translocation.

MD simulations were performed with the sander module of AMBER 6.0 45. The periodically replicated rectangular unit cell contains four CNTs that are hexagonally packed in the x, y plane and between 1183 and 1188 TIP3P water molecules 46. Two types of CNTs of (10,10) armchair type were used: water-permeable tubes with open ends, and impermeable tubes with capped ends. The tubes consist of 360 (open) and 340 (capped) sp² carbon atoms 47,48 and are 20.7 (open) and 20.4Å (capped) long with a radius of 6.7Å (tube axis to carbon nuclei). This setup results in system sizes of ∼5000 atoms. We used two types of membranes, one consisting exclusively of open tubes (“sieve” geometry), and the other combining three capped and one open tube within a unit cell in the x, y plane. The resulting system mimics a channel in an otherwise sealed membrane and is denoted as “channel” geometry. The two setups are illustrated in Fig. 1. Each unit cell contains a single Na⁺ ion 49. For comparison and to complement the results of the simulations with a single ion, three additional types of systems were simulated: 1), the above (10,10) CNT membranes without ions; 2), the (10,10) CNT membranes with five pairs of Na⁺ and Cl⁻ ions (initially placed randomly in the bulk water phase); and 3), sieve and channel systems consisting of narrower (6,6) CNTs (4Å radius) with a single Na⁺ ion and 516 TIP3P water molecules. Simulations of systems with five ion pairs corresponding to an ionic concentration of ∼0.2 mol/L will allow us on one hand to compare free energies from an umbrella sampling procedure to unbiased simulations. On the other hand, the comparison of “infinite dilution” and “finite concentration” simulations will also allow us to assess possible ion-concentration dependences in the low-concentration regime. Finally, simulations at finite concentration will give us insights into the rate of ion translocation across the membrane.

In all simulations reported here, the preassembled CNT membranes remained stable with only sub-Ångstrom fluctuations out of plane. All simulations were carried out with a time step of 2 fs; the simulation temperature was kept constant at 300K by weak coupling to a temperature bath with a relaxation time of 0.2ps 50. The cutoff distance for Lennard-Jones interactions was set to 10Å. Particle-mesh Ewald summation was used for electrostatic interactions 51. For an equilibration period of 4ns, the systems were simulated at a constant pressure of 1atm using the weak coupling algorithm with a relaxation time of 0.5ps. Anisotropic pressure scaling was applied, i.e., all three box lengths were allowed to fluctuate independently. After equilibration, production simulations and free-energy calculations were carried out at constant volume. Simulation lengths for analysis of equilibrium properties of the free simulations (without umbrella potentials) were 10ns (with ions) and 12ns (without ions), respectively.

The free energy of an ion as a function of the axial distance z from the membrane center and the radial distance r from the tube axis was determined by umbrella sampling 43. A harmonic biasing potential K_r,i(r – r_i)²/2+K_z,i(z – z_i)²/2 was applied to the ion, with counter forces distributed uniformly over all carbon atoms of the CNT membranes. r and z are the axial distance from the membrane center and the radial distance from the tube axis respectively; r_i and z_i are the target positions; and K_r,i and K_z,i are the corresponding force constants in the ith umbrella-sampling window. In simulations with channel geometry, the ion was pulled into the single open tube. In the sieve geometry, one specific tube was selected. The ion-membrane target position z_i was varied from 20 to 0Å in 0.5-Å increments with a force constant K_z,i of 2.0kcal mol⁻¹Å⁻². In the radial direction, three types of restraints were used. Restraints with a target position r_i=0Å and two different force constants K_r,i of 0.1 and 0.07kcalmol⁻¹Å⁻² were applied to all z_i windows. Additionally, restraints with a target position r_i of 5Å and a force constant of 0.1kcalmol⁻¹Å⁻² were applied to windows with 12Å≤z_i≤20Å (where z=0 corresponds to the membrane center). Altogether, this results in 99 overlapping windows with 390-ps sampling time in each window (38.6ns total). The weighted histogram analysis method 44,52,53 was used for a collective analysis of the data.

The (unbiased) probabilities of the ion position P(z, r) that are directly obtained from the weighted histogram analysis method are radial distribution functions in the r direction and thus weighted by the size of the volume element 2πr. We compute the “local” free energy G(x, y, z)=–k_BT ln(P(z, r)/(2πr)), where the distances x and y from the tube center in the direction of the corresponding axes satisfy r=(x²+y²)^1/2. k_B is Boltzmann’s constant, and T is the absolute temperature.

To obtain a one-dimensional potential of mean force (PMF), i.e., the free-energy profile G(z) along the pore axis, we integrate P(z, r) in the x, y plane. In the Results section, we compare the resulting PMFs, where P(z, r) is integrated from r=0 to R=5Å, i.e.,

To analyze the hydration structure of ions in the bulk phase and in the center of the nanotubes, additional simulations were carried out with ions restrained only in the z direction (K_z,i=2.0kcalmol⁻¹Å⁻², z_i=0 and 20Å).

Local electrostatic potentials play a central role in the translocation of ions through molecular pores. In explicit-solvent simulations, those potentials are determined by the fluctuating molecular structure of water in the bulk phase, the membrane interface, and the pore, as well as the position of ions and other charged or polarizable groups. In our simulations, we determine an average local charge-density distribution ρ(r) by binning atomic charges of the solvent and ions into a three-dimensional grid. The local electrostatic potential (and field) generated by ρ(r) is determined by the Poisson equation 54

(1)

(2)

A continuum electrostatics description approximating the system in the channel-type simulations consists of a rectangular periodic box filled with high-dielectric (ϵ_s=78) medium and a low-dielectric (ϵ_i=1) slab with a cylindrical (high-dielectric) pore. The ion is treated as a sphere with an internal permittivity of ϵ_i=1 and a single charge at its center. We use the iterative, FFT-based method to solve the Poisson equation for such a system 56,57. The region of low dielectric permittivity (ϵ_i=1) consists of the interior of the ion and the low-dielectric slab. Geometric parameters are the size of the periodic box (32×28.8×51.2Å, corresponding to the simulated system), the thickness of the low dielectric (ϵ_i=1) slab (20Å), the pore radius (5.3Å), and the ionic radius (1.5Å). Algorithmic parameters are applied as described before 56,57. A grid spacing of 0.2Å was used.

The solvation free energy of the ion was computed for several ion positions inside the box. As indicated in Fig. 4, we systematically varied the distance from the membrane center in the z direction and the distance from the channel axis in the x, y direction. The absolute solvation free energies obtained by this grid-based algorithm depend somewhat on the position of the ion center relative to the grid. However, since the aim of the study was not to compute absolute solvation free energies of the ion, these dependences could be largely eliminated. Grid effects were suppressed by subtracting from the solvation free energy of the ion in the system with membrane the solvation free energy of the ion at the same coordinate in an otherwise identical system without membrane.

In this section, we first determine the free energy for pore translocation of ions by umbrella sampling simulations. Then, we compare the energetics from all-atom simulations to those obtained from continuum electrostatics calculations. An analysis of the local water structure, electrostatic potentials, and internal energies will provide insights into the observed thermodynamics of ion translocation. Finally, the mobility of the ion in the pore and in bulk solution will be compared.

Fig. 2 shows the two-dimensional free-energy surface G(z, r) of a sodium ion in a carbon nanotube as a function of the axial distance z of the ion from the membrane center and the radial distance r from the tube axis, as defined in the Methods section. We find that the free-energy difference between the bulk water phase and a narrow cylindrical region in the channel interior is small overall. Inside the tube, a radial distance of ∼2Å from the tube axis is energetically favored for both geometries (channel and sieve). For the channel geometry, this region is slightly higher in free energy than the bulk water phase; for the sieve geometry, this region is slightly lower (and thus favored) compared to the bulk. In contrast, ion positions near the tube axis (r∼0) are highly unfavorable. (Note that the data in Fig. 2 have been corrected for the “Jacobian” 2πr.)

Display large version of this figure

Figure 2

Free-energy surfaces for ion and water penetration into membrane pores. (Upper panels) Two-dimensional free-energy surface G(z, r) of a sodium ion in a CNT as a function of the axial distance of the ion from the membrane center z and the radial distance from the tube axis r as obtained from umbrella sampling (left, sieve geometry; right, channel geometry). (Lower left panel) Free-energy surface for water determined from the average density as –k_BT ln(ρ/ρ₀)+4kJ/mol (where ρ₀=0.033Å⁻³; same scale as the ion free energies). The white bar indicates the positions of the carbon atom centers. (Lower right panel) Radial ion, water, and carbon distribution in the tube (averages z=−5 to 5Å; ion and water densities presented relative to the corresponding density in the bulk phase; carbon density in arbitrary units). As discussed in the Methods section, the “Jacobian” 2πr has been removed, so that uniform ion and water distributions would produce flat free energy profiles.

Fig. 3 shows the one-dimensional PMF of the ion projected onto the z coordinate by integrating the probability in the radial dimension from r=0 to r=5Å. We see that the resulting apparent free-energy barrier for the ion to pass through a nanotube with an effective diameter of ∼10Å is remarkably low (∼3–6kJ/mol).

Display large version of this figure

Figure 3

Free energy of ion translocation through membrane pore. The two-dimensional PMF of the ion was projected onto the z coordinate by integrating the probability in the radial dimension from r=0 to r=5Å. Results are shown for MD simulations of the wide nanotubes (black solid line, umbrella sampling/channel system; red solid line, umbrella sampling/sieve system; red dashed line, free simulation/sieve system) and continuum electrostatics (cyan dashed line). (Purple line and inset) MD simulations of narrow tubes in channel system.

Since the difference in the two-dimensional free energy between points in the bulk water phase and some regions inside the nanotube is small, a major contribution to the PMF is the restriction of the accessible space in the x, y direction, once the ion enters the tube. Not accounted for in Fig. 3 is an additional entropic contribution of k_BT ln 4=3.5kJ/mol that makes entering a tube more likely in a simulation with sieve geometry (four open tubes) than in one with channel geometry (one open tube).

Beyond this trivial entropy contribution, we find that it is slightly more favorable for the ion to enter a tube in a sieve system than in a channel system. As Fig. 3 shows, the free-energy barrier is ∼2kJ/mol lower. This barrier reduction can be explained by the fact that the neighboring nanotubes in the sieve geometry are filled with water. The presence of this high-dielectric medium inside the membrane decreases the dielectric penalty for the ion to go into the pore by increasing the effective dielectric constant of the surrounding membrane.

For comparison we also analyzed 10ns of unrestrained simulations (i.e., with a freely moving ion without umbrella potential). As predicted by the umbrella sampling simulations, the barrier for the ion to pass though the nanotube membrane is low and several passages can be observed for the sieve systems. In a simulation with a single Na⁺ ion, the ion enters one of the tubes once; in a simulation with five Na⁺/Cl⁻ ion pairs, one of the Na⁺ ions enters a tube once and stays inside for ∼1.5ns, and later in the course of the simulation one Na⁺ ion passes through a tube. We find at most one ion at a time in the membrane pores. A crude estimate for the PMF based on the simulation with five ion pairs in a sieve geometry is included in Fig. 3. It reproduces the magnitude of the free-energy barrier obtained from umbrella sampling. No passages are observed for the free simulations of the channel system, confirming that the free-energy barrier is considerably higher because of the lower number of open pores and the higher free-energy barrier of the single pore, as discussed. We note further that we did not observe any Cl⁻ ions penetrating into the nanotube pores.

Also included in Fig. 3 is the PMF for the corresponding systems with narrower pores formed by (6,6) nanotubes (pore radius 4Å with respect to carbon nuclei, corresponding to a ∼2.5Å “effective” radius accessible to water, as obtained by subtracting the “radius” of the carbon atoms lining the pore). Here, the free-energy barrier estimated by umbrella sampling is considerably higher (∼120kJ/mol instead of 5kJ/mol). As in the wider pores, the barrier for the sieve geometry is only marginally lower (by ∼3kJ/mol). In summary, a hydrophobic pore with an effective radius of ∼5–5.5Å does not pose a high energetic barrier for a Na⁺ ion. In contrast, a narrow pore with an effective radius of ∼2.5Å poses a very high barrier.

This dramatic change in ion permeability with a relatively small change in the pore diameter provides an explanation for potassium-channel gating by pore-width modulation 58. The hydrophobic gating region of the KcsA-channel structure in the closed form has a diameter of ∼4Å at its narrowest point, whereas the gating region of the MthK structure in the open form has widened considerably—with a diameter of 12Å at its narrowest point 22. Based on our estimates for a much simpler pore, such a change lowers the barrier for ion translocation from >100kJ/mol (40k_BT) to near zero.

In this section, we explore to what extent a continuum electrostatics description of the system captures the essence of the all-atom free-energy calculations. A comparison of continuum calculations and explicit-solvent simulations will also allow us to explore the role of the solvent molecular structure for channel electrostatics. Fig. 4 shows the continuum electrostatic contribution to the solvation free energy of an ion (single charge; radius 1.5Å) for several positions of the ion inside the box. As sketched in Figure 4A, the distance from the membrane center in the z direction and the radial distance from the channel axis in the x, y plane is systematically varied. The free-energy profiles in the z direction are shown in Figure 4B for several radial distances. The free energies are given relative to the solvation free energy of the ion at the center of the bulk water phase, where no dependence from the position relative to the channel axis (x and y coordinates) is observed. General observations are as expected: 1), the free-energy increases as the ion enters the channel and gets closer to the membrane center, thus getting deeper into the low-dielectric slab; 2), the free-energy barrier increases as the ion approaches the channel wall, i.e., as the layer of high dielectric, polarizable medium between the ion and the pore wall gets thinner. In Figure 4C, we compare the free energies from umbrella sampling and continuum electrostatics for a two-dimensional cut through the membrane (as a function of z and r). Overall, the free energies from the continuum electrostatics calculations are in a range comparable to those determined from explicit-solvent MD simulations. However, as a function of the distance from the channel axis, the free energy of the ion in the explicit-solvent simulations differs significantly from the continuum-electrostatics predictions. As expected by symmetry, the continuum electrostatics calculation finds the pore axis to be the most favorable location for the ion inside the pore. In explicit water, the ion almost exclusively stays in a cylindrical layer ∼2Å from the channel axis. To “coarse grain” such molecular effects, we also compare radially integrated probabilities from continuum electrostatics and MD simulations (from r=0 to 5Å). The resulting PMFs are shown in Figure 4B and Fig. 3. The overall barrier is approximately twice as high in a continuum electrostatics description as in the simulations with explicit solvent, but with an absolute difference of only ∼4–5kJ/mol (∼2k_BT). The free-energy profiles obtained from continuum electrostatics also have a somewhat different shape than those from explicit-solvent simulations. The latter level off after the ion has entered the pore, corresponding to a flat free-energy barrier, whereas the profiles from continuum electrostatics have a clear maximum at the membrane center. It should be noted, however, that a quantitative comparison of the two free-energy profiles requires some caution, since the continuum electrostatics results depend on model parameters such as tube and ion radius (i.e., the definition of the dielectric boundary), and the dielectric permittivity of the high dielectric medium. For example, the PMF obtained from the continuum model assuming an ionic radius of 2Å is ∼1.5kJ/mol lower, and is thus in better agreement with the MD data than the barrier height obtained for an ionic radius of 1.5Å. Similarly, the barrier height obtained assuming a dielectric constant of 94, which corresponds to the TIP3P water model 59, is 0.4kJ/mol lower than with a dielectric constant of 78 (experimental dielectric constant of water). For comparison, lowering the dielectric constant to 50 increases the barrier by ∼1.5kJ/mol.

Display large version of this figure

Figure 4

(A) Setup for continuum-electrostatics calculations: rectangular periodic box with high-dielectric (ϵ_s=78) bulk medium and low-dielectric (ϵ_i=1) slab with a cylindrical (high-dielectric) pore (10.6-Å diameter). Ion positions are varied from z=0 to z=25Å at distances from the tube axis of r=0Å (black), 2Å (red), 3Å (green), 3.8Å (blue), and 5.3Å (magenta). (B) Continuum-electrostatics free-energy profiles along the z direction for several radial distances from the tube axis (color code as in A) and PMF obtained from radially integrated probabilities (dashed cyan line). Note that the free-energy curves were shifted to zero at the center of the bulk phase. (C) Comparison of free energies from umbrella sampling (thin lines) and continuum electrostatics (thick lines).

The shape of the free-energy profile obtained from explicit-solvent simulations of the membranes with narrow pores agrees well with that obtained from continuum electrostatics (not shown). However, the height of the barrier is substantially lower in the continuum case. This difference reflects mostly difficulties in defining a “correct” continuum system for a narrow molecular pore. Naively taking the parameters of the narrow-pore system (effective tube radius 2.5Å, ion radius 1.5Å) and solving the continuum electrostatics problem results in a much too low free-energy barrier of ∼40kJ/mol compared to ∼120kJ/mol from MD simulations. The reason is that in this simple approach the thin space between ion and membrane is filled by high-dielectric medium, which leads to dielectric screening that is not present in the “real” system. There, the space between ion and carbon wall is empty (except for a single water molecule). We realize that such small, empty volumes are more accurately accounted for in some continuum electrostatics simulations, for example through the rolling probe algorithm 60. This discussion highlights possible problems of continuum models for narrow cavities.

In summary, continuum electrostatics correctly predicts that a hydrophobic pore of ∼10Å in diameter poses a comparatively low free-energy barrier for an ion, although not as low as was found in the explicit-solvent simulations. Continuum electrostatics also predicts that ions face a much higher barrier for the narrow pores. However, without optimizing the boundaries between high and low dielectric, that barrier is substantially lower than that of the explicit-solvent simulations. Moreover, for the wide pores continuum electrostatics does not explain the preferred path of ions through the channel—not at the center but in a narrow, cylindrical shell close to the pore wall. To explain this behavior, we next study the local water structure inside the pore and analyze the hydration shell of the ion.

In the lower right panel of Fig. 2, we compare the ion and water densities inside a nanotube as a function of the distance r from the tube axis. The ion densities are obtained from the umbrella sampling simulations (after unbiasing). Water (and nanotube carbon) densities are obtained from free simulations of the nanotube systems without ions. We find that water forms cylindrical shells, with a core of high water density in the tube center and a high-density layer along the tube wall (there is no difference in the water density profiles inside the nanotubes between sieve and channel systems). Similarly layered water structure inside narrow hydrophobic pores 61,62,63,64 or small cavities 65 has been observed in previous studies. We find that the ion is located predominantly in the cylindrical space between the two regions of high water density. This finding reproduces results by Lynden-Bell and Rasaiah 61, who investigated ion diffusion in hydrophobic tubes of different radii 61 using a rather different simulation setup (an infinite tube generated by periodically replicating a tube segment in one dimension; a repulsive cylindrical wall; minimum-image convention; and a different water model).

Fig. 5 shows the average local electrostatic potential on a cut through two nanotube centers of a sieve system parallel to the xz plane. (All observations made in the following section are qualitatively identical for the channel system.) In the simulation, the ion is restrained in the z direction so that it remains close to the membrane center. Note that by adding the restraint, only the z position of the ion deviates from equilibrium. From the restraint simulations, we can thus evaluate the equilibrium averages of the potential (and later ion hydration structure) conditional on the ion being at the membrane center. The left panel of Fig. 5 shows the total electrostatic potential, in contrast to the electrostatic potential generated by the solvent charges, shown on the right. In Fig. 5, the left tube contains the ion, which generates a region of high positive potential (on average) at the preferred location of the ion, a ring around the tube axis of radius ∼2Å. The right tube is filled solely with water, and the local potential inside is a good representation for the local potential inside a nanotube in a simulation without any ions (not shown). Similar to the structure in the water density (Fig. 2, lower panels), one also observes a layer structure in the electrostatic potential. A one-dimensional potential profile along the x and z axes respectively, is displayed in Fig. 6 (calculated from simulations without ions), together with the distribution of oxygen and hydrogen atoms. We see that the layers of low (negative) electrostatic potential coincide with the regions of high oxygen and hydrogen density, and the layers of higher electrostatic potential coincide with the regions of low water density and with the regions of the low dielectric membrane interior. A small excess of hydrogen atoms is pointing toward the carbon walls and toward the low water density region, leading to a net polarization and the observed structure in the local electrostatic potential. The structure of the hydrogen bond network at inert or hydrophobic surfaces and the fraction of hydrogen atoms that point toward this surface has been analyzed extensively 66,67,68,69,70, leading to an electrostatic surface potential of ∼+0.5V 70,71, consistent with the results shown in Figure 5 and Figure 6. This surface potential results from the ordering of water molecules without any external field and cannot be accounted for by continuum electrostatics.

Display large version of this figure

Figure 5

Average local electrostatic potential on a cut through two nanotube centers of a sieve system parallel to the xz plane with the ion inside the left tube (restrained to the membrane center in the z direction). (Left panel) Total electrostatic potential. (Right panel) Electrostatic potential generated by the solvent charges.

Display large version of this figure

Figure 6

Average total local electrostatic potential (A and B) and water atom distribution (C and D; solid lines, oxygen; dashed lines, hydrogen) of free MD simulations of the sieve system without ions. (A and C) Profile along the x direction through the centers of two CNTs. (B and D) Profile along the z direction through the low-dielectric/carbon-rich region of the membrane between two nanotubes. The coordinates of the cuts and the width over which the potential was averaged are indicated in A and B.

Does the electrostatic potential in the absence of the ion determine the location of the ion in the pore? With the given layer structure, a positively charged Na⁺ ion should occupy the tube region with negative electrostatic potential. Instead, as shown in Figure 2 and Figure 5, the ion is preferentially found in a region with initially positive potential but low water density. This result suggests that it is thermodynamically more favorable for the ion to move to a region of low water density but unfavorable positive potential, thus forcing a reorientation of the surrounding water molecules, than to move to the region of high water density, where it displaces water and disturbs its structure. This observation illustrates that conclusions based on average electrostatic potentials, without consideration of fluctuations, need to be drawn with caution, especially in the case of potentials due to solvent polarization without external field or ionic charge. Introducing a charged particle amounts to a large perturbation that may alter the solvent polarization qualitatively. This polarization effect can be seen on the right side of Fig. 5, where the electrostatic potential of the polarized solvent is displayed. Water reorientation leads to a highly negative potential at the site of the ion and the water structure is considerably perturbed compared to the nanotube without ion. The solvent polarization leads to an almost perfect screening of the ionic charge, so that in contrast to the electrostatic potential generated by the solvent, the (average) total electrostatic potential remains almost unperturbed by the ion except for its immediate vicinity (compare panels in Fig. 5).

A closer inspection of the water structure around the ion and the energetics of the ion in the tube and in the bulk solution may help to explain why the ion so readily moves into the tube. Fig. 7 presents the radial distribution function of the water oxygen atoms around the ion in a channel geometry. In the simulation, the ion is restrained in the z direction, with reference distances z₀ to the membrane center of 0 (inside the tube) and 20Å (outside the tube/in the bulk water phase). Fig. 7 shows that the extent of the first solvation shell—defined as the number of oxygen atoms that have a distance to the ion of ≤3Å—is practically identical in the tube and in the bulk phase. The cumulative distribution function in the inset shows that the first solvation shell contains between five and six water molecules, which is slightly higher than the average coordination number for Na⁺ of 4.6 predicted by ab initio calculations 72. The average coordination numbers are given in Table 1, together with the fraction of structures that have five and six water molecules, respectively, in the first solvation shell. (In all cases the fraction of structures where the ion is coordinated by <5 or >6 water molecules is <2%.) Inside the nanotube, the fraction of ions with six water molecules in the first solvation shell is lower than in the bulk water phase. Fig. 8 shows the radial distribution of ions inside the nanotube with five- and six-membered solvation shells. We find that ions with five-membered solvation shells are preferentially found closer to the tube wall. Ions with six-membered solvation shells are only found at distances from the wall that are large enough to fit a water molecule between the ion and carbon wall. In Fig. 8, snapshots of the ion and its first solvation shell inside the nanotube are presented for both five- and six-coordinated ions. The spatial separation of the two states of the ion with six- and five-membered solvation shells is reflected in the correlation time of the fluctuation between the two states. As can be seen in Fig. 9 and Table 1, the interconversion between five- and six-coordinated states is much slower in the pore than in the bulk water phase (∼5 vs. 1ps relaxation time). With a two-state analysis, we find that the dominant effect is a ∼5-fold slowdown of the transition from five to six water molecules in the pore compared to the bulk. One possible explanation, which combines this result and the radial distribution of the two states shown in Fig. 8, is that the ion with a five-membered solvation shell “sticks” to the tube wall, and needs to “dissociate” to bind to a sixth water molecule.

Display large version of this figure

Figure 7

Radial distribution function of the water oxygen atoms around the ion inside (solid line) and outside (dashed line) the tube (normalized by the density of TIP3P water of 0.033Å⁻³) in the simulation of the channel system. (Inset) Cumulative distribution of oxygen atoms around ion.

Display large version of this figure

Figure 8

Distribution of ions inside the nanotube as a function of the distance r from the tube axis for ions with five (solid black bars) and six (striped bars) water molecules in their first solvation shell. In the channel-geometry simulation, the ion was harmonically restrained near the membrane center (along the z direction, but not radially). Insets show representative structures.

Display large version of this figure

Figure 9

Normalized correlation function of fluctuations between five- and six-coordinated ions. The ions were harmonically restrained near the membrane center (z₀=0Å) and outside the tube (z₀=20Å). Solid and dashed lines are for channel and sieve simulations, respectively.

To investigate the energetics of ion permeation into the pores, we calculate their average internal energy 〈U〉 (electrostatic and Lennard-Jones interactions of the ion with the rest of the system). Table 2 shows the means and the variances of the internal energy of the ion inside the tube and in the bulk water phase for the channel and sieve systems of wide (10,10) and narrow (6,6) CNTs. The average internal energy of the ion is lower in the bulk phase than inside the pore, the difference being ∼10kJ/mol for the wide tubes and ∼300kJ/mol for the narrow tubes. Following perturbation theory 73,74, higher-order terms of the energy distribution can be used to obtain a rough estimate of the free energy:

(3)

We find that in the wide pores, the slightly higher average energies of the ion inside the pore are compensated in part by a narrowing of the energy distribution. A similar compensating effect was observed for water confined in nanotubes 47. For the narrow pores, the average energy lost by an ion entering the pore is large, ∼300kJ/mol, and that loss cannot be compensated by the narrower distribution. However, the simple second-order perturbation theory is not quantitatively reliable 75 for such broad energy distributions.

We also estimated local diffusion coefficients of the ion inside the pore along z direction and in the bulk fluid. From the positional correlation of ions in the umbrella sampling simulations we can estimate the diffusion coefficient as

(4)

The main results of our simulation study are that narrow nonpolar pores (less than ∼5Å effective diameter at a length of ∼21Å) block ions, but only slightly wider pores (∼10-Å effective diameter) are highly permeable. In the narrow pore, ions get stripped off their solvation shell, resulting in a large energetic penalty; in the wider pore, the first shell remains intact, and the energetic penalty for translocation reduces essentially to the entropic loss for confining an ion in a pore. Remarkably, in the wide pore, we find not only a low barrier, but also a high ionic mobility along the pore axis, comparable to that of an ion in bulk water.

This strong dependence of the pore permeability on the pore diameter can provide a powerful gating mechanism in biological ion channels. Based on our simulations, a relatively modest change in the pore diameter is sufficient to convert a nonpolar pore between conducting and nonconducting states. Indeed, recent crystallographic studies of potassium channels showed hydrophobic pores of less than ∼5-Å diameter in the closed state 7 and ∼12-Å diameter in the open state 22.

With our simple model system, we could also explore to what degree continuum electrostatics calculations capture the energetics of ion translocation through nonpolar membrane pores. We find that the continuum calculations correctly predict that a hydrophobic pore of ∼10-Å effective diameter poses only a comparatively low free-energy barrier for an ion. Whereas the absolute barrier is somewhat higher than that found in the explicit-solvent simulations, relatively modest changes in the atom radii can improve the agreement.

However, probed in more detail, continuum electrostatics does not explain the preferred path of the ions through the channel—not at the center but in a narrow, cylindrical shell close to the pore wall. To explain this behavior, one needs to take the water structure inside the hydrophobic pore into account. As has been suggested elsewhere, the confined water phase inside narrow hydrophobic cavities has unique structural 61,62,65 and dynamic 27,47,64,81,82,83 properties that seem to affect water, proton, and ion transport through membrane channels.

Two studies on ion permeation through nanopores by Beckstein et al. 84,85, which came to our attention only after completion of this work, come to the similar conclusion that the molecular structure of the water phase inside the nanopore is crucial to the understanding of the ion permeation process. At a qualitative level, the results are similar, but Beckstein et al. 84 appear somewhat more pessimistic about the continuum results. They report only one-dimensional free-energy profiles, similar to the G(z) obtained here by integration in the radial direction, and find that the continuum barrier is dramatically higher for narrow pores and slightly lower for wide pores. These qualitative differences from our results (of a continuum barrier that is higher than the MD barrier for the wide pore, but lower for the narrow pore) can probably be explained by differences in the continuum calculations. In particular, we use periodic boundary conditions in both MD and continuum calculations; we do not use a “rolling probe” correction to define the high-dielectric region in the narrow tube; and we do not use the Poisson-Boltzmann equation with a finite ion concentration (neither for the ion in the pore nor in the bulk water phase).

We find that water inside the nanotube forms cylindrical shells of high and low water density, and that the Na⁺ ion is predominantly found in an interstitial cylindrical layer of low water density at a distance of ∼2Å from the pore axis. This can be seen, for instance, in the radial distribution profiles shown in Fig. 2. Between the two water shells, the ion maintains an intact first hydration shell. The ion needs to replace fewer water molecules in the low-density region inside the tube than in the bulk water phase, which is an energetic advantage that partly compensates for the loss of the solvation contributions of the second hydration shell and beyond. Other studies of ions at the water/vapor interface have shown that ions can be found at interfaces, as long as they can maintain their first hydration shell 86.

As reported in a study by Dzubiella et al. 87, the transport properties of an ion through a nanopore are strongly linked to the water structure and the water density in the pore. Qualitative differences between the results of Dzubiella et al. 87 and those of our study in terms of ion permeability of a nanopore can be explained by the smaller effective pore radius in 87, which leads to a different water structure inside the nanopore 61,64,87. This observation further illustrates the sensitivity of ion permeability to modifications in the pore width and the resulting changes in the water structure. (However, in this context it is worth noting that comparisons of pore radii can be difficult because of differences in the definition, and are best performed by using the radial water density profiles.)

The results from continuum-electrostatics calculations, which correctly predict the magnitude of the free-energy barrier for ion translocation through a hydrophobic pore, confirm that a major contribution to this free-energy barrier is due to electrostatic interactions and the different dielectric environment outside and inside the pore. However, atomistic simulations show that the relationship between pore radius, structure of the water phase inside the pore, and ion permeability is rather complex. The local water structure—which cannot be easily accounted for in a simple dielectric model—potentially affects gating mechanisms in biological ion channels.

In conclusion, we want to emphasize that we have explored various methodologies to simulate ion translocation through membrane channels, in particular the treatment of long-range electrostatic forces (lattice-sum or reaction-field techniques) and the influence of system size and periodic boundary conditions. We found that different methodologies do not change the results qualitatively. A detailed study of these aspects will be presented elsewhere (C. Peter and G. Hummer, unpublished data).