INTRODUCTION

TOP ABSTRACT INTRODUCTION THEORY AND METHODS RESULTS DISCUSSION CONCLUSIONS REFERENCES

Synthetic ion channels are emerging as promising nanoscale biosensor devices (Bayley, 1999; Hartgerink et al., 1998), and self-assembled channels of the kind considered in this paper have demonstrated potent antibacterial properties (Fernandez-Lopez et al., 2001). The ease of introducing and evaluating design changes in synthetic systems provides an opportunity to test computational models of ion channels, thus advancing the quantitative understanding of analogous biological systems. The present study is guided by this aim. An earlier work addressed via molecular dynamics (MD) simulations the dynamics of water in a self-assembled, cyclo[-(L-Gln-D-Ala-L-Glu-D-Ala)₂-] octapeptide nanotube (Engels et al., 1995). In the present investigation, we use both MD and continuum electrostatics to study the energetics of transferring an ion from the bulk to specified locations inside a nanotube comprising cyclo[-(L-Ala-D-Ala)₄-] monomers.

Molecular dynamics has been extensively used in biomolecular simulations, and the exemplar of this technique for channel systems are the investigations by Roux and Karplus (1991, 1993), who used an atomic description of the channel, water, and ions to study the translocation of water and alkali metal ions through gramicidin. (See, for example, Hao et al. (1997) and Åqvist and Luzhkov (2000) for more recent applications of MD in calculating ion translocation energetics through biological channels.) Although MD is now routinely used in biomolecular simulation, it still requires a significant investment of time and computational resources. For example, realistic membrane-protein simulations such as those reported by Tieleman and Berendsen (1998) or Bernéche and Roux (2000) are still much too computationally intensive to be used on a routine basis. Furthermore, if one is interested not in the actual dynamics but only the energetics of placing an ion inside an ion channel, and if this estimate is required for a number of different modifications of the channel, then simplified alternatives to MD become desirable.

The use of non-MD approaches to study ion partitioning in channel systems has not been as extensive. Åqvist and Warshel (1989) studied ion translocation through the gramicidin channel using the protein-dipole Langevin-dipole formulation for the channel protein, membrane, and bulk water but used explicit water molecules inside the channel itself. They were able to obtain activation barriers that were in reasonable agreement with experimental estimates. Jordan and coworkers have pursued a more reduced, semimicroscopic description of channels (Partenskii and Jordan, 1992a,b; Partenskii et al., 1994; Sancho et al., 1995): the atomistic protein structure was disregarded, the water external to the channel was treated as a high dielectric continuum, whereas the water inside the channel was represented by freely rotating dipoles. Using this description, they illustrated the nonlinearity of the electric response upon the introduction of an ion into the channel. In a recent study, Roux and MacKinnon (1999) investigated the selectivity of a KcsA (K⁺) channel using continuum electrostatics and were able to show how the ion is stabilized in the low dielectric protein/membrane environment.

In this work we present atomistic MD simulations of ion-channel systems and formulate computationally facile continuum electrostatics models based on insights gleaned from the simulations. The model channel is a self-assembled cyclic octapeptide with an alternating (D,L)-amino acid sequence (Hartgerink et al., 1998; Ghadiri et al., 1994). These cyclic peptides self-assemble into a cylindrical tube with a -sheet (backbone) architecture. The nanotube channel provides a more complicated, yet experimentally accessible, system compared with gramicidin. Like gramicidin, the nanotube is made of alternating (D,L)-amino acids. Further, like gramicidin (Arseniev et al., 1985), the nanotube pore wall is solely lined by carbonyl and amide dipoles. However, there is one critical difference, the nanotube pore diameter of ~6 Å is twice that of gramicidin. (These pore diameters were estimated with Oliver Smart's HOLE program (Smart et al., 1993).) Thus, unlike gramicidin, the nanotube channel contains more than a single column of water molecules (Engels et al., 1995), and the ion can, in principle, coordinate with more than two water molecules, which is the maximum possible in the gramicidin channel. The larger diameter also leads to cation conduction rates about three times greater (Ghadiri et al., 1994) than that in gramicidin. The pore diameter of the nanotube falls in the middle of the range of diameters noted for the nicotinic acetylcholine receptor channel, which has a pore diameter ranging between 3 Å to ~9 Å (Opella et al., 1999) and a greater variety of groups lining the pore.

We study the water-tube transfer energetics of Li⁺, Na⁺, Rb⁺, and Cl, initially using MD-based perturbation free energy methods, thereby exploring the solvation structure and thermodynamics of particles with a range of sizes and charge. To better study the structure and thermodynamics of ion solvation in nanotubes, we focus exclusively on the solvation of a solitary ion inside the channel. An important aspect of our work is the treatment of long-range electrostatic interactions using Ewald summation. To the best of our knowledge, this is the first study to use perturbation free energy simulations including Ewald-electrostatics for ions in nanotubes. Earlier investigations of channel systems have typically used various truncation schemes to study the dynamics and energetics of ion transport (but see Bernéche and Roux (2000)). We find that the structuring of water molecules around the ion is sensitive to the location of the ion inside the channel. We present macroscopic models at various levels of detail and highlight the interrelation of water structuring and continuum modeling. Water structuring should manifest as a nonlinearity of the medium's electrostatic response; however, we find it difficult to unequivocally discern water structuring based solely on the calculated electrostatic response. The possible reasons for this are discussed. The present work demonstrates the subtlety and importance of water structuring in the narrow confines of ion channels and provides ways to describe these effects in continuum electrostatics calculations.

	THEORY AND METHODS

TOP ABSTRACT INTRODUCTION THEORY AND METHODS RESULTS DISCUSSION CONCLUSIONS REFERENCES

Molecular dynamics simulations were performed with the CHARMM (Brooks et al., 1983) molecular mechanics package (version c26b2) using the all-atom PARAM22 forcefield (MacKerell et al., 1998). The Lennard-Jones parameters for Na⁺ and Cl are from Beglov and Roux (1994), and those for Li⁺ and Rb⁺ are from Åqvist (1990). The Li⁺ and Rb⁺ parameters, originally used with the simple point charge water model (Åqvist, 1990), were used with the TIP3P (Neria et al., 1996; Jorgensen et al., 1983) water model in CHARMM. (In CHARMM, the water protons have a small radius as opposed to the original TIP3P (Jorgensen et al., 1983) water model.)

The continuum electrostatics calculations were performed with the MEAD (Macroscopic Electrostatics with Atomic Detail) (Bashford and Gerwert, 1992; Bashford, 1997) suite of programs. The atomic partial charges from the PARAM22 forcefield and the corresponding atomic radii compiled by Nina et al. (1997) were used. The Born radii of the ions used in the MEAD calculations were determined from simulations as discussed below.

Molecular structure

The cyclic octapeptide c(_DAA)₄ used in this study is a simplification of the c(_DLW)₄ octapeptide used experimentally (Sanchez-Quesada and Ghadiri, personal communication; see also Ghadiri et al., 1994; Motesharei and Ghadiri, 1997). The monomer c(_DAA)₄ peptide was made from a D-Ala-L-Ala dimer by applying a fourfold rotation about the z axis. Initial bond lengths and bond angles in the D-Ala-L-Ala dimer were taken from the equilibrium bond lengths and bond angles of the PARAM22 forcefield (MacKerell et al., 1998). Initial dihedral angles were taken from the crystal structure of the blocked octapeptide dimer c(_LW_DA)₄ (Ghadiri et al., 1993). This cyclic peptide was energy minimized by the adopted basis Newton-Raphson (ABNR) method, which was carried to convergence; the fourfold symmetry was maintained during the energy minimization.

To create an octapeptide dimer having a backbone structure matching, as closely as possible, the experimental crystal structure, a duplicate of the octapeptide was generated, rotated by 180° about a line passing through the -carbon atoms of residues 1 and 5, translated by 5.006 Å along the z axis and finally rotated by 9.17° about the z axis. The resulting octapeptide dimer was energy minimized using the ABNR method. The energy minimization was carried out in a symmetrical infinite-tube environment generated by translating dimer images ±9.4 Å along the z axis. The 9.4 Å was arrived at by a series of similar energy minimizations using different z translations from which it was found that 9.4 Å gave the lowest final energy. Crystallographic studies of the blocked dimer suggest a unit cell dimension of ~4.7 Å (Ghadiri et al., 1993) corresponding to a dimer translation of ~9.4 Å, in good agreement with what we find via the energy-minimization protocol.

From this energy-minimized octapeptide dimer, a 4-mer nanotube was created in two stages. First the octapeptide dimer was filled with a canonical count of water (Engels et al., 1995): one water molecule in each -plane, and two water molecules in each midplane (Fig. 1). This dimer was energy minimized in a symmetrical-infinite environment. Second a 4-mer nanotube was built from this energy-minimized, water-filled dimer by appropriate translational symmetry operations. The choice of a 4-mer nanotube is made for computational convenience. In experiments the tube length would depend on the bilayer used, which for a DMPC bilayer would amount to having six to eight monomers, based on the intermonomer spacing of ~5 Å.

View larger version (34K): [in this window] [in a new window]   FIGURE 1   Schematic of the 4-mer. The -plane of each monomer is defined as the plane coming closest (in the least square sense) to passing through the -carbon atoms of the monomer. The midplanes are defined as the planes midway between the -planes. (A) Slightly off-center view along the axis of the channel. Only the innercoordination shell of the ion is shown for clarity. (B) View perpendicular to the channel axis. The front of the channel has been cut away for clarity. The water molecules are arranged in a canonical fashion: two water molecules in the midplane and one in the -plane.

Molecular dynamics

The water-filled nanotube is solvated in a preequilibrated box of TIP3P (Neria et al., 1996; Jorgensen et al., 1983) water molecules. Any water molecules external to the nanotube whose oxygen atoms were closer than 2.8 Å were deleted as overlapping waters. A Na⁺ ion was placed in the midplane region bounded by -planes B and C (Fig. 1). The volume of this system was adjusted based on a partial specific volume of 1 cc/gm for water and ion and 0.73 cc/gm for the nanotube (Engels et al., 1995). Overlaps with the internal water molecules were removed by minimization. The system so generated measures 32.3 Å × 32.3 Å × 47.07 Å and comprises 320 nanotube atoms, 1549 water molecules, and the sole ion. In the initial course of this investigation we used this large simulation volume. However, comparable accuracy is obtained even with a small box that is derived from this system by deleting selected water molecules to give a cubic box with edge length 28.948 Å. This smaller system comprises the nanotube, ion, and 718 water molecules. From this initial configuration, three additional boxes were obtained by replacing the Na⁺ ion with Li⁺, Rb⁺, and Cl ions, respectively, followed by 50 steps of ABNR minimization. The four boxes thus obtained were the starting configurations for further studies.

In the initial phase of equilibration, each of the systems obtained above were equilibrated for more than 120 ps. During the equilibration cycle and for the rest of the simulation, the tube backbone atoms, i.e., the N (imino), -carbon, and carbonyl-carbon atoms were restrained with a harmonic (Hookian) spring of force constant 2m kcal/mol-Ų, in which m is the mass of the respective atom in atomic mass units. The ion was held in place by a harmonic spring of force constant 10m kcal/mol-Ų. Dynamics were propagated using the Leapfrog Verlet (Allen and Tildesley, 1997) algorithm using a time step of 1 fs. The SHAKE (Ryckaert et al., 1977) algorithm was used to rigidly fix bonds to hydrogen atoms. The temperature was maintained at 298 K by coupling to a heat bath with a coupling constant of 0.5 ps (Berendsen et al., 1984). The nonbonded interaction list was truncated at 12 Å on an atom basis; the van der Waals interactions were switched to zero from 10 Å to 12 Å; and Ewald sums with conducting boundary conditions were used for electrostatic interactions with the Ewald damping coefficient  set to 0.2333 Å¹. The reciprocal space cutoff was determined by performing test calculations to determine a suitable value. This lead to a spherical cutoff such that ²  38(2/L)². , the reciprocal space lattice vector, is given by  = 2(l/L_x, m/L_y, n/L_z) (Feller et al., 1996), in which L_x, L_y, and L_z are the box dimensions. Because the Ewald formulation is meaningful only for a neutral system, a uniform background compensating charge density is added to ensure electroneutrality of the simulation cell. This charge density is just the negative of the net charge smeared over the entire simulation volume. (The interested reader should consult Tosi's (1964) clear presentation of Ewald summation in the presence of a neutralizing background.)

In the next phase of equilibration, the ion charge was slowly changed from the full formal charge to zero charge as follows. A series of equilibrated trajectories with ion charges of q (q = ±1 a.u.;  = 0.95, 0.85, ... , 0.05) were generated, in which the starting configuration of each trajectory was taken from the final configuration of the earlier trajectory. For example, the  = 0.85 trajectory began with the ending coordinates of the  = 0.95 trajectory. For each of the trajectories, initial velocities were assigned from a Gaussian distribution such that the initial temperature was 625 K. By scaling the particle velocities, the systems were cooled to 298 K over 1 ps followed by 1 ps of equilibration with temperature maintained at 298 K by velocity rescaling. This "scrambling" of water configuration was primarily to prevent water molecules from getting trapped in some particular configuration, as this is a matter of concern when the water molecules are confined to a narrow channel. These systems were further equilibrated for 8 ps with temperature maintained at 298 K by coupling to a heat bath. The nonbonded interactions were treated as mentioned above. The end point of the  = 0.05 trajectory was taken as the starting point for subsequent perturbation free energy calculations. We adopted this approach to minimize hysteresis in subsequent perturbation calculations, as the present approach would likely ensure that the all-important water configuration around the ion is not trapped in a state corresponding to the ion at a higher charge state. This is validated by the modest error bars on the calculated free energies (Results section). Note that the perturbation free energy calculations were done in two different ways. As discussed below, one method for calculating the free energy incorporated the "scrambling" of water at each charge state ; the other did not, but the simulation extended for much longer, namely 300 ps.

The starting configuration for the simulations with the ion in the -plane was obtained via a somewhat different protocol. We used the small 28.948 Å simulation box from the very start, instead of using a bigger box and deleting water molecules. The ion was inserted in -plane B (Fig. 1) and the ion and tube atoms restrained as before. This system was equilibrated for 200 ps using a 2-fs timestep, followed by a 300-ps perturbation free energy calculation, again using a 2-fs time step. Bonds to hydrogen atoms were fixed using SHAKE (Ryckaert et al., 1977). The charge was turned off in steps of 0.05, to obtain the  = 0.05 state. At each step a 10-ps equilibration was followed by 20 ps of data collection. Thus, this procedure of obtaining the  = 0.05 state also provided an estimate of the discharging free energy (data not shown). All subsequent perturbation calculations for the ion in the -plane were as for the ion in the midplane, i.e., using a 1-fs timestep and included a "scrambling" at each charge state. Likewise a few 300-ps charging perturbation calculations using a 2-fs timestep were performed for the ion in the midplane (data not shown). Adopting different styles of perturbation and different timesteps permitted us to evaluate the effect of these variables on calculated numbers. The average energetic values obtained using the 300-ps perturbation calculations (data not shown) are similar to the values reported in this work indicating that the calculations are reasonably well converged.

Solvation energetics

(1)

(2)

(3)

1

Solvation structure

Analysis of electrostatic response

(4)

(5)

Continuum electrostatics

The macroscopic description of ion solvation is based on treating the peptides and the bulk water as continua with dielectric constants of  = 4 and  = 80, respectively. The ion with  = 1 is represented as a spherical cavity with a charge at its center. The molecular surface (Connolly, 1983), determined by rolling a 1.4 Å radius probe sphere over the atomic radii, delineates regions of different dielectric constants. The potential is described by the Poisson equation, which must be solved numerically because of the complex shape of the dielectric boundaries. This is implemented with the above-mentioned MEAD software.

The electrostatic contribution to the excess free energy of the ion in the nanotube, with the ion in vacuum ( = 1) as the reference state, is given by

(6)

in which _rxn is the reaction field acting on the ion in the multidielectric environment of the channel and the water with the all channel charges set to zero; _channel is the potential at the position of the ion charge due to the peptide channel charges; and q is the charge of the ion itself. _channel is independent of the ion charge, whereas _rxn is proportional to the ion charge. Their sum, _rxn + _channel, is the continuum analog of of Eq. 4. Here the  dependence enters through the ion-charge dependence of _rxn. Eq. 6 above for the ion in bulk water is just A_bulk = (1/2)q_rxn, where now _rxn is q/a · (  1), which is the Born model for ionic solvation. The factor of (1/2) in Eq. 6 is simply a consequence of _rxn depending linearly on charge.

The Poisson equation was solved by a finite difference approach as implemented in the MEAD (Bashford and Gerwert, 1992; Bashford, 1997) suite of programs. The initial coarse grid had 51 grid points with a grid spacing of 1 Å and the fine grid had 101 grid points with a grid spacing of 0.1 Å. Notationally we refer to this grid choice by 51 × 1:101 × 0.1. We have also performed calculations at grid level 51 × 1:101 × 0.25:111 × 0.1 with similar results.

Some of the macroscopic calculations included a limited number of solvation shell water molecules as part of the ion channel environment. In these cases, the assignment of water coordinates was performed as follows. Initial water coordinates were taken from the last frame of the MD simulation done in calculating the solvation structure. Water molecules not part of the first solvation shell were deleted and the nanotube coordinates reset to their original values, and the system was energy minimized (using ABNR) for 250 steps with the nanotube backbone and the ion harmonically restrained. Configurations with various fractional charge states were generated by energy minimization starting with the structure at the preceding charge state. Thus, for example, the structure corresponding to the  = 0.75 state was derived by energy minimization of the structure corresponding to  = 1.0. We used this method as a means of allowing explicitly included solvent shell water to respond to the charging process by explicit reorientation, while keeping the nanotube harmonically restrained, thus confining its response to that modeled by the dielectric constant of 4.

	RESULTS

TOP ABSTRACT INTRODUCTION THEORY AND METHODS RESULTS DISCUSSION CONCLUSIONS REFERENCES

Molecular dynamics

The calculated excess free energies of charging the ions in the bulk solvent are presented in Table 1 along with experimental values and results of similar calculations by others. Because we are also interested in continuum electrostatic calculations, we relate the bulk hydration free energies obtained via microscopic simulations to their macroscopic counterpart, the Born model for ion solvation. This allows us to calculate the only free parameter in the Born model, the ionic radius, by the formula (Born, 1920),

(7)

in which A_bulk is given in kcal/mol and q in atomic units, and a is the Born radius in Å. The radii thus obtained are intimately related to how water coordinates with the ion and are usually much different from the crystal (or Pauling) radii (for example, see Latimer et al., 1939). Molecular interpretations of the Born radius are usually based on the first peak of the ion-water radial distribution function (Roux et al., 1990) or an average of the crystal radius and the first-peak of the radial distribution function (Babu and Lim, 1999). In this study, together with the microscopic charging free energy calculation of A_bulk, we adopt Eq. 7 as the definition of an effective Born radius, a.

The electrostatic component of the transfer free energy into the nanotube (columns labeled MD; Table 3) is obtained as the difference between the charging free energies in nanotube (Table 2) versus bulk (Table 1). It is clear that the cations tend to partition readily into the tube, whereas the anion does not. This trend is in agreement with experimental observations (Sanchez-Quesada and Ghadiri, personal communication). There are two reasons for the favorable partitioning of the cations: enhanced ion-water interaction and favorable interaction of the cation with the nanotube. The issue of ion-water interaction is addressed below, and the role of nanotube-ion interaction is addressed in the Discussion section.

View this table: [in this window] [in a new window]   TABLE 2   Excess free energy (in kcal/mol) for charging the ion in nanotube mid- and -planes

The radial distribution of water around the ions (Fig. 2) indicates that the first shell is somewhat narrower around the ion inside the nanotube than in the bulk. This indicates more tightly bound water molecules in the channel. On the other hand, in comparison with bulk, the first shell in the channel has fewer water molecules, indicating desolvation (Table 4). The coordination number is determined from the integral of the radial distribution function. The angular distribution of water dipoles around the ion is a subtler metric of ion-water interaction (for example, see Jayaram et al., 1989; Hyun et al., 1995). In the channel midplane (Fig. 3), for Li⁺ and Na⁺, the water around the ion in the nanotube is more highly structured than in bulk with angles mostly confined to very favorable dipole orientations (cos  > 0.75), whereas for Rb⁺ in the midplane distribution of dipoles is bulk-like. In the channel -plane (Fig. 3) the distribution is bulk-like for Li⁺, whereas enhanced structuring is observed for Na⁺ and Rb⁺. Similar calculations at a charge state of  = 0.5 indicates that the angular distribution of water dipoles around ions of all types are similar in the bulk and inside the nanotube (data not shown). Thus, any structuring beyond that found in bulk appears when ion charges are above ~50% of their full formal charge.

View larger version (0K): [in this window] [in a new window]   FIGURE 2   Ion-oxygen (water) radial distribution function, g(r), in bulk (dotted line) and inside the tube (solid line). (A) Ion in the midplane. (B) Ion in the -plane. The curves have been arbitrarily shifted along the ordinate to fit on the same graph and hence no values are specified for the ordinate. Inside the tube, the absence of a pronounced second maxima is an artifact of using a radial distribution for waters that are constrained by the cylindrical geometry of the tube: beyond the first shell the g(r) has little meaning inside the tube. The raw distribution functions have been smoothed (Allen and Tildesley, 1997).

View larger version (1K): [in this window] [in a new window]   FIGURE 3   Probability distribution, p(cos ), of the orientation of water dipoles around the ions. Dotted lines, ion in bulk; solid line, ion in nanotube midplane; and broken line, ion in the nanotube -plane.

For Cl (Fig. 3) in the channel midplane, the angular distribution of water is quite different from the bulk. In the bulk the angular distribution peaks around cos   0.6. This corresponds to a linear hydrogen bond between Cl and water (see Solvation structure in theory and methods). This means that in bulk, the bifurcated hydrogen bond to the ion, corresponding to cos  = 1, is less favored than the single H-bond that leaves the other H-atom free to H-bond with other water molecules. These results are consistent with those of Ichiye and coworkers (Hyun et al., 1995). (Their definition of  differs by  from our definition. Also, if we normalize our distributions (data not shown), the peak magnitudes are in the same range as theirs, but the agreement is not exact, as our Cl parameters are somewhat different than theirs.) In contrast, the distribution in the nanotube midplane shows a large fraction of water molecules in the bifurcated hydrogen bond configuration. In the -plane, however, the angular distribution of water dipoles is similar to that for the ion in the bulk albeit with some skew toward the bifurcated H bond (Fig. 3).

The data points in Figs. 4 and 5 show the potential calculated using Eq. 5, i.e., the electrostatic response function both inside the nanotube and in the reference bulk solvent. These data are fit to quadratic functions of  in the plots. However, the quadratic terms are much smaller than the linear terms, and the curves can be said to be nearly linear.

View larger version (29K): [in this window] [in a new window]   FIGURE 4   Average potential, Eq. 5, at the ion in the midplane at charge state . (), Ion in the nanotube; (×), ion in the bulk solvent. The solid line is least squares fit to the bulk data. The dashed line is least squares fit to the nanotube data. The fitted equations are also shown.

View larger version (29K): [in this window] [in a new window]   FIGURE 5   Average potential, Eq. 5, at the ion in the -plane. Rest as in Fig. 4.

At the  = 0 state, the calculations predict a non-zero, positive potential at the site of the (uncharged) ion in bulk, consistent with results found by earlier investigators (Pratt et al., 1994; Hummer et al., 1996; Ashbaugh and Wood, 1997). The potential predicted at the site of Na⁺⁰ is ~9.0 kcal/mol-e, in good agreement with the range of values observed by earlier investigators for particles of similar dimensions. As a test of our simulation protocol, we calculated ₀ at Na⁺⁰ in a box of SPC water molecules. The Na⁺ ion potential parameters were those of Straatsma and Berendsen (1988). We get ₀ = 10.3 kcal/mol-e and A = 94.8 kcal/mol, which agrees well with ₀ = 10.1 kcal/mol-e and A = 95.7 ± 0.6 kcal/mol-e reported by Ashbaugh and Wood (1997) using the same set of parameters (data not tabulated). Likewise, we calculate ₁ = 164.7 kcal/mol-e before and ₁ = 215.3 kcal/mol-e after finite-size correction (data not tabulated), respectively. These values are also in good agreement with those read from Fig. 1 in Ashbaugh and Wood (1997), further reinforcing our confidence in our calculations.

Continuum electrostatics

The electrostatic component of the excess free energy of the ion inside the nanotube, calculated via two different macroscopic models, is shown in Table 2. These models are schematically illustrated in Fig. 6. Model A is the simplest: the water inside the nanotube is assumed to be a uniform medium of dielectric constant 80. The ion radii are the Born radii from Table 1. Model B explicitly considers the first shell water molecules as part of the low-dielectric nanotube environment, i.e. in model B, the first-solvation shell water molecules are modeled atomistically and assigned  = 4, equal to dielectric constant of the nanotube. The water molecules beyond the first-solvation shell are described as a continuum of  = 80. The number of water molecules is set equal to the coordination number observed in the simulations (Table 4), and at each charge state ( = 0.25, 0.5, 0.75, 1.0), the position of the water molecule is reoptimized to mimic the response of a molecular fluid instead of a continuous medium (see Methods).

View larger version (56K): [in this window] [in a new window]   FIGURE 6   Models of ion solvation within the macroscopic formulation. The environment outside the tube has a dielectric constant of 80. Model A, channel lumen has a dielectric constant of 80; model B, explicit water molecules are included as part of the nanotube environment.

Table 3 shows the calculated bulk to nanotube transfer free energies, and Figs. 7 and 8 depict the solvation free calculated via the continuum models against their atomistic counterpart. For Li⁺ and Na⁺ in the midplane and for Na⁺ and Rb⁺ in the -plane, model B describes the simulation results better than model A. (The difference between MD results and continuum results are smaller for model B.) For Li⁺ in the -plane, models A and B have errors of approximately equal magnitude but opposite sign. For Cl both models A and B are inadequate for the ion in the midplane, and in the -plane model A better describes simulation results.

View larger version (23K): [in this window] [in a new window]   FIGURE 7   Comparison of free energies (kcal/mol) calculated via macroscopic models versus MD simulations for ion the mid-plane. (), Model A; (), model B. In each panel, the charge states are successively  = 0.25, 0.5, 0.75, and 1.00 from bottom left-hand corner to top right-hand corner.

View larger version (23K): [in this window] [in a new window]   FIGURE 8   Comparison of free energies (kcal/mol) calculated via macroscopic models versus MD simulations for ion the -plane. Rest as in Fig. 7.

	DISCUSSION

TOP ABSTRACT INTRODUCTION THEORY AND METHODS RESULTS DISCUSSION CONCLUSIONS REFERENCES

We briefly discuss the bulk free energy calculations before addressing the main concerns of this work, namely the formulation of continuum models based on insights from simulations. As Table 1 illustrates, although the trends in the solvation free energies are similar to experimental estimates, there is disagreement in magnitude between our calculations using Ewald summation and those performed with ions in a water droplet surrounded by a vacuum (Åqvist, 1990; Beglov and Roux, 1994; Roux, 1996). Using Ewald summation, we find that the charging free energy for cations is ~10 kcal/mol less negative than the corresponding value obtained in simulations using water droplets. For Cl the Ewald summation results are more negative by ~15 kcal/mole than the corresponding values from water-droplet simulations. Darden et al. (1998) find a similar discrepancy and have suggested water polarization at the droplet-vacuum interface as a possible cause. Thus, although we use the same parameters, using Ewald electrostatics leads to numbers (predictably) different than those with droplet simulations. However, the relative trends are consistent between different methods. Likewise, for charging the ion in the nanotube, this difference between droplet simulations and our Ewald approach can be expected to hold. Thus, in calculating the free energy change for bulk to tube partitioning, both Ewald and droplet approaches should lead to similar results. This is heartening, as to the best of our knowledge the Ewald-electrostatics approach for charging ions in nanotubes is novel and the trends obtained using a different method would be similar to those presented in this paper.

As Table 3 indicates, cation partitioning is favored into the nanotube. As cation size is increased, the partitioning free energy is reduced. Experimentally, K_m, the affinity constant, for ion-nanotube interactions suggests that Li⁺ readily partitions into the nanotube, Na⁺ does to a lesser extent followed by Rb⁺ (Sanchez-Quesada and Ghadiri, personal communication). This trend is validated by our MD simulations (Table 3).

The role of the nanotube in determining the changes in free energies for transition from the mid- to -plane is two-fold. One aspect is the nanotube's influence on ion water interactions and the second is the interaction of the ion with nanotube partial charges. Engels et al. (1995) have shown that the water occupancy of an ion-free nanotube is ~two to three water molecules in the midplane and ~one to two water molecules in the -plane. These numbers reflect how much space there is to accommodate water molecules in each zone, and the observed coordination numbers and the attendant changes in free energies can be considered in this light.

Li⁺ in the midplane coordinates with ~five water molecules (Table 4), three within the midplane region and one each from the -plane above and below. In the -plane, one might guess that Li⁺ would coordinate with six waters, three each from the adjacent midplanes. However, the calculated value is four, two each from the adjacent midplanes. It is likely that tighter coordination with fewer water molecules is preferred over the converse. Clearly a loss of one water molecule influences the less favorable free energy of solvation in the -plane (Table 3). For Na⁺ there is no change in the coordination number in going from mid- to -plane (Table 4), and the change in free energy in going from mid- to -plane is not as great (Table 3). For Na⁺, the midplane to -plane transition energy is consistent with earlier studies by Demchuk and Bashford using a completely different umbrella-sampling technique (unpublished observations). The smaller N_coor, in comparison with Li⁺, in the midplane is likely due to the larger size of particle. For Na⁺ in the -plane, as for Li⁺, it appears that a smaller, tightly coordinated first shell is preferred. For Rb⁺, N_coor increases in the -plane (Table 4) and the free energy change in going from mid- to -planes is also favorable (Table 3). A similar trend is seen for Cl (Tables 3 and 4), and the sharp change in free energy in going from midplane to -plane is again consistent with the calculations by Demchuk and Bashford (unpublished observations).

The role of the direct nanotube-ion interaction can be seen by decomposing _channel from model B into contributions from the first solvation-shell water molecules, _water, and the nanotube, _tube (Table 5). The changes in free energy suggested by changes in coordination number are validated by the changes in _water: increased coordination numbers are reflected in more favorable ion water interactions. Further, _tube is uniformly negative throughout the channel. This clearly influences the preferential solvation of cations inside the tube. Also _tube decreases in magnitude from the mid- to the -plane. This is due to the diminished effect of the negatively charged carbonyl oxygens lining the midplane and an enhanced effect of the positively charged C atoms.

Table 5 also reveals expected trends for _rxn, the energetic response of the dielectric. For Na⁺ and Cl, the Born radii are nearly equal, and correspondingly _rxn is also nearly equal, but opposite in sign. Further, the Born equation, Eq. 7, suggests that _rxn should be larger for smaller radii, as is seen in Table 5. Finally, on account of desolvation and because of the low dielectric nanotube, the reaction field component of the free energy would be less negative inside the nanotube; correspondingly (1/2)_rxnq is less negative than the bulk value (Table 1) as expected.

We emphasize that in our calculations, the ion does not have a contact interaction with the peptide channel (Fig. 1), as it is always restrained near the axis of the channel. That such a calculation is credible is seen from the calculations by Crozier et al. (2001). In their simulations of a model channel of diameter 8.1 Å, it is seen that the ion hardly moves farther than 1 Å from the channel axis. This study is relevant to the present discussion, as the choice of dipoles lining the channel wall and the chosen diameter are similar to that in the nanotube.

Macroscopic models and simulations

Model A is a purely linear response model because 1) the water inside the tube is described as a continuum dielectric and 2) the nanotube has a fixed structure. The simulations in the bulk and in the nanotube suggest some deviation from linearity (Figs. 4 and 5), but the nonlinearity is quite modest. Hence, our intuitive expectation would be that a linear response model would hold for the ion in the bulk and in the nanotube. However, calculations with model A belie this expectation. In some instances the model A results are closer to MD and in others the model B results are closer. Below we reconcile this observation with information obtained from simulations. To facilitate the discussion we refer to the  = 0.75, 1.0 states as the high-charge limit, and the  = 0.25, 0.5 states as the low-charge limit.

Tables 6 and 7 shows the deviations of models A and B from the simulation results. In the low-charge limit it is seen that model A almost always describes the simulated trends better than model B. In this limit, as we noted earlier, the solvent structuring around the ion is bulk-like. Further, the solvent is fairly disordered when ion charges are low (data not shown; see Jayaram et al., 1989). Hence, in the low-charge limit either model B is inappropriate or using a low dielectric constant for the first shell water is inappropriate. This conclusion is similar to that reached by Partenskii and Jordan (1992a,b). They find that in the low-field limit the dipolar chain in their model of an ion-channel with a single file of water has a high dielectric constant. In their model the low field limit is for ion charges less than 20 to 30% of their full formal value, a reflection of the single file of water; a wider channel should shift this limit higher, as it does in the nanotube.

View this table: [in this window] [in a new window]   TABLE 7   Calculated root mean square error for free energies of charging calculated via model A and model B in the -plane

For the cations in the midplane, Fig. 3 indicates that in the high-charge limit, the first solvation shell around Li⁺ and Na⁺ is structured, whereas it is bulk-like for Rb⁺. Correspondingly, in this limit, for Li⁺ and Na⁺ model B describes the data better (Table 6 and Fig. 7) and model A does so for Rb⁺. In the -plane, the solvent around Li⁺ is bulk-like, whereas it is more structured around Na⁺ and Rb⁺. Correspondingly, for Li⁺ model A best describes the data in the high-charge limit, and model B does so for Na⁺ and Rb⁺ (Table 7 and Fig. 8).

Cl solvation provides an interesting counterpoint to the studies with cations. For the cations, the hydrogen bonding pattern with water is similar for the ion in the bulk and in the nanotube, only the magnitude of ordering changes. For Cl in the midplane, the orientational ordering of water is much different than that seen in bulk with a trend towards bifurcated hydrogen bonding with water (Fig. 3). In the high-charge limit, model B is clearly a poor descriptor of solvation and so is model A (Fig. 7 and Table 6). This is understandable in light of the fact that the Born radius was calculated with linear hydrogen bonding dominant. It may indeed be appropriate to recalculate a Born radius for solvation via bifurcated hydrogen bonds, but we have not undertaken this task as our aim was to formulate a consistent model with a limited number of adjustable parameters. The above point about redefining Born radii becomes clearer when one considers the solvation of Cl in the -plane. Here the orientational ordering is fairly similar to that seen in the bulk solvent (Fig. 3), and because no pronounced structuring is seen, model A might be expected to describe the MD trends better. Indeed this is what is seen (Fig. 8 and Table 7).

In the above analysis, we have tried to correlate the trends from continuum models with solvent structuring. Since saturation (nonlinearity) in electrostatic response must go hand-in-hand with increase in structural ordering, the lack of a significant nonlinearity in the calculated response (Figs. 4 and 5) is somewhat unanticipated. One plausible explanation is that inaccuracies in our simulation mask any nonlinearity. Thus, when the linearity or nonlinearity of the electrostatic response is unable to provide a clear direction in the choice of continuum models, the structural calculations appear to be better guides in choosing an appropriate continuum model.

Thus, of our simulation approaches, one gives us structural information regarding ion solvation, and these simulations can be performed with relative ease. Based on these simulations, if one finds enhanced structuring of water around the ion relative to bulk values, then it is fair to expect a model like our model B to be a better model for describing the solvation thermodynamics. This is a principal result of our paper. The other simulations are computationally more expensive perturbation free energy calculations. These simulations can be analyzed to understand the electrostatic response, and based on that one can construct appropriate continuum models (see Jayaram et al., 1989). However, in the present case, we find that due to likely inaccuracies in our simulations, it is difficult to use the electrostatic response (a derivative of a free energy) as a sole guide in the formulation of continuum models.

Choice of dielectric constants

In our continuum electrostatic calculations, the only free parameter is the assumed dielectric constant of the nanotube and the first shell water molecules, the peptide charges and radius parameters having been determined from PARAM22 and the ion Born radii from the bulk simulations, respectively. The calculations discussed above were all performed with this dielectric constant set to 4. Calculations with dielectric constants of 3 and 5 have also been done (data not shown). In all these cases, the qualitative trends predicted are as discussed earlier. Quantitative differences, however, exist: With a dielectric constant of 3 model B uniformly overpredicts solvation, but the differences are all within 5 kcal/mol. Similarly with a dielectric constant of 5, the values are uniformly underpredicted.

The assignment of a low dielectric constant to the first solvation shell water molecules is based on their being partially immobilized, as evidenced by the dipole distributions. In particular, model B, in which the first-shell water and the peptide are assigned the same low dielectric constant, represents an assumption that the first-shell water dipoles are immobilized to approximately the same extent as the peptide dipoles. Interestingly, in the high-field/high-charge limit Jordan and coworkers (Partenskii and Jordan, 1992a,b) find that the dipolar chain in their model of the ion channel is best described with a dielectric constant between 3 and 5, a value that also depends on the location of the ion inside the tube. This is so because in the high field limit the dipolar chain is nearly immobilized.

Proteins often have been modeled as media of dielectric constant 4. This choice has found support in theoretical estimates of dielectric constants (Gilson and Honig, 1986; Simonson et al., 1991; Smith et al., 1993