INTRODUCTION

TOP ABSTRACT INTRODUCTION METHODOLOGY RESULTS AND DISCUSSION CONCLUSIONS APPENDIX REFERENCES

Aqueous proton transport (PT) is implicated in a broad variety of important biological processes, including the passage of protons to or from the inside of a cell via ion channels. Ion channels are usually formed by integral membrane proteins that contain water molecules in their interior, through which protons and other ions may diffuse to enter or exit the inside of a cell (Hille, 1992).

Most information about ion channels is obtained from electrophysiology experiments, where the current-voltage relations of channels can be measured as a function of the type and concentration of the ion that carries the current (see, for example, Phillips et al., 1999; Cukierman, 2000). However, it is often difficult to interpret the data generated from these experiments at more than a qualitative level. Molecular dynamics (MD) simulation, however, is able to provide quantitative information on the likely structure-function relationship of channels (Levitt, 1999). For example, classical MD simulations have been used to study synthetic leucine-serine channels embedded in phospholipid membranes, including a narrow system that is proton-selective (Randa et al., 1999). Although no attempt was made to explicitly incorporate hydrated protons into this study, water molecules were observed to enter the channels during the simulations and to form hydrogen-bonded networks through which PT could in theory occur by a Grotthuss-type hopping mechanism (Agmon, 1995). In fact, with the exception of an equilibrium quantum mechanical study of PT in the gramicidin channel (Pomès and Roux, 1996a), there have been few accounts of MD simulations of actual aqueous PT in ion channels. Aqueous PT in channels is therefore a subject that is worthy of further investigation at the level of MD simulation.

Part of the difficulty associated with the simulation of aqueous PT is the need for a potential energy surface that can accurately describe the reactive process as proton passes (or "hops") through the three-dimensional hydrogen bond networks that are formed by water molecules. At present, this behavior is beyond the capabilities of the empirical potential surfaces that are commonly used in biomolecular MD simulations. Also, quantum mechanical effects may be evident in properties of aqueous PT, mostly due to the light mass of the transferring particle. However, Schmitt and Voth (1998, 1999a, b) have recently addressed these issues with the development of their multi-state empirical valence bond (MS-EVB) potential energy surface.

Several other potentials and methods have also been developed and applied to the problem of PT in water. In particular, Vuilleumier and Borgis have developed an alternative empirical valence bond surface (Vuilleumier and Borgis, 1999), while Marx et al. have applied a combined path integral molecular dynamics/Car-Parrinello simulation technique to study equilibrium properties of an excess proton and 32 water molecules (Marx et al., 1999). Protonated water wires have also been studied using the older and less accurate PM6 potential (Stillinger and David, 1978): Pomès and Roux have considered equilibrium quantum mechanical properties (Pomès and Roux, 1996b) and (classical) free energy barriers to PT (Pomès and Roux, 1998), while Hammes-Schiffer and co-workers have performed non-equilibrium quantum/classical dynamics simulations (Decornez et al., 1999).

Aspects of these potentials have been discussed elsewhere (Schmitt and Voth, 1998, 1999a, b). Suffice it to state here that the MS-EVB model is both sufficiently accurate and numerically efficient for it to be used in the classical and quantum calculation of equilibrium and dynamical properties of a variety of protonated water systems, ranging from small clusters to larger bulk-like situations. This potential has been found to give a good description of many properties of the excess proton in liquid water. For example, it is able to reproduce experimental observables such as the hopping rate and density of vibrational states and, importantly, it provides a good description of the relative stabilities of the solvated Eigen H₉O and Zundel H₅O species (Agmon, 1999). As a result of the extensive validation studies for bulk water, the MS-EVB potential holds considerable promise for the study of PT in realistic biological systems.

Although there is considerable interest in modeling actual biological ion channels at an atomistic level of detail (see above), there is also a need for the study of even more simplified models for channels. Indeed, such models have previously been used to provide useful and more general insight into equilibrium and dynamical properties of, for example, water molecules alone in channels and cavities (Sansom et al., 1996; Hartnig et al., 1998) and also aqueous solutions of various ions in channels (Allen et al., 1999; Lynden-Bell and Rasaiah, 1996). The aim of the present work is therefore to study aqueous PT in model channels using classical MD simulation and the MS-EVB potential. We use the smooth cylindrical hydrophobic channel potential function that has been used in a comprehensive study of the solvation and mobility of ions in channels (Lynden-Bell and Rasaiah, 1996). With this approach one can consider key properties as a function of channel radius, such as the relative stabilities of solvated proton structures (particularly the Eigen and Zundel cations) and the mobility of the excess charge, so as to provide a series of reference simulations for comparison to more realistic models for ion channels.

The present study differs from earlier work on proton wires (Pomès and Roux, 1996a, 1998; Decornez et al., 1999) in at least two important ways. First, as mentioned earlier, an accurate and validated potential for an excess proton in water has been used. Second, a range of pore sizes for the model proton channel are studied instead of a single-file water chain as defined by construction. In the larger pore radius systems, the proton transport occurs over an increasingly three-dimensional hydrogen-bonding network in the channel. As a result, we are able to observe a rather distinct and quite interesting transition to a quasi-one-dimensional proton wire (with greatly enhanced proton mobility) as the channel pore size decreases. The results presented here should be of general significance for all biological proton channel systems, but they may be especially relevant to understanding proton fluxes through membranes in the absence of an explicit biomolecular channel (Nichols and Deamer, 1980; Nagle, 1987; Deamer, 1987; Marrink et al., 1996). In these systems it is believed that transient "pores" may be present in the lipid bilayer into which water wires may form, thus transporting protons through the membrane.

This paper is organized as follows. In the next section the MS-EVB potential energy surface will be introduced and a coordinate for the center of the excess protonic charge will be described. Details of the MD simulations will also be given. Static and dynamic properties of the aqueous PT channel systems are then presented in the following section and discussed, along with a consideration of the limitations of our models and the relevance of the present results. Conclusions that may be drawn from this work will be given in the last section.

	METHODOLOGY

TOP ABSTRACT INTRODUCTION METHODOLOGY RESULTS AND DISCUSSION CONCLUSIONS APPENDIX REFERENCES

Potential energy surface

The MS-EVB potential energy surface is now briefly summarized to provide the reader with a description of terms that are used in the following sections, while a detailed description of the surface can be found elsewhere (Schmitt and Voth, 1999a). In general, any EVB technique rests upon the construction and diagonalization of a [N × N] Hamiltonian matrix H, such that the potential energy is given by the lowest eigenvalue E₀ of this matrix:

(1)

where c is the ground state eigenvector with elements c_i [i = 1, N], and H is represented in a basis of N valence bond type states (Warshel, 1991). The diagonal elements of the matrix are then approximated with empirical potential functions of the nuclear coordinates, of the type more commonly used in MD simulations. The functional forms of the off-diagonal elements are parametrized so that Eq. 1 can reproduce the potential surface in regions where the standard potentials fail, e.g., the transition state region(s) of the reaction. The instantaneous forces on the system nuclei may be computed from the adiabatic ground state in Eq. 1 using the Hellmann-Feynman theorem (Schmitt and Voth, 1999a).

In the present case of an excess proton and (n) water molecules, each empirical valence bond (EVB) state is taken to consist of a hydronium cation and (n  1) solvent water molecules. For example, in the gas phase two EVB states are required to represent the Zundel cation H₅O, while four states are required to represent the Eigen cation H₉O. Illustrations of the Eigen and Zundel cations can be found in several articles by two of us (Schmitt and Voth, 1998, 1999a, b). In bulk, EVB states extending out to the second solvation shell of the so-called "pivot" hydronium are included in the construction of H. The pivot hydronium is that which belongs to the EVB state with the largest amplitude (i.e., the largest value of c [i = 1, N] in Eq. 1), and can thus be regarded as the most hydronium-like species in a given configuration. EVB states beyond the second solvation shell have only very small amplitudes and can be safely ignored, so that on average 10 states are used to construct H for an excess proton in liquid water. This multi-state approach allows for the unrestricted hopping of an excess proton through the hydrogen-bonded network of water molecules.

In the MS-EVB model the diagonal elements of H are composed of inter and intramolecular contributions from the hydronium cation and water molecules, which are modeled with the modified flexible TIP3P potential (Dang and Pettitt, 1987). The functional form of the off-diagonal matrix elements takes account of the interaction between the exchange charge distribution of distinct protonated dimers H₅O within EVB complex and all remaining solvent water molecules. This interaction was found to be especially important in bulk simulations in order to reduce the localization tendency, or self-trapping, of (diagonal) hydronium-like EVB states. Finally, parameters entering the MS-EVB model were fit to reproduce ab initio global minimum geometries and energies of selected protonated water clusters H⁺(H₂O)_n (n = 2, 4), so as to combine the accuracy obtained from high-level quantum chemistry calculations with the convenience of the MS-EVB formalism.

To confine the excess proton and water systems to cylindrical hydrophobic channels, a repulsive potential was used as introduced by Lynden-Bell and Rasaiah (LBR) (Lynden-Bell and Rasaiah, 1996). The potential, given by

(2)

was applied to all oxygen atoms, where r_O is the perpendicular distance from a given oxygen atom to the channel axis, and r_cyl is the radius of the channel. Following LBR, the same values of A (= 38.72 kcal/mol) and B (= 4.18 Å¹) were used for both solvent water and solute (e.g., hydronium) oxygen atoms. In this way the interaction of all oxygen atoms in the system with the channel wall is the same for each EVB state, so that this interaction need only be calculated once, after the construction and diagonalization of the EVB Hamiltonian matrix.

Center of excess charge

Because of the delocalized nature of the positive charge associated with the excess proton (Hynes, 1999; Agmon, 1999), it is necessary to have a coordinate that reflects this property when considering the location of the excess "proton" species. Such a coordinate can be defined as a linear combination of the amplitudes associated with each resonance EVB state that is used to represent the Hamiltonian matrix H. A "center of excess charge" (CEC) coordinate, given by

(3)

has been used in this work. Here, r_i(t) are coordinates of the center of mass of the hydronium cation in the ith EVB state at time t, and the amplitudes c [i = 1, N] are obtained from Eq. 1.

Simulation details

Some physical properties of the model channels are given in Table 1, including the average number of EVB states used per potential energy evaluation. In practice, an integer number of EVB states are used to calculate the potential energy and forces within a given configuration; the fact that the average number of EVB states is non-integer reflects the dynamic nature of the algorithm that is used to generate the EVB states (Schmitt and Voth, 1999a). As in other works (Allen et al., 1999; Lynden-Bell and Rasaiah, 1996), we have opted for the simple procedure of maintaining a constant density of water molecules in the channels (of volume rl) by treating the channel length l as an adjustable parameter.

For each channel model, initial configurations of water molecules and a hydronium cation were generated at random inside cylinders of radius r_cyl. Each random configuration was then equilibrated for 75 ps at 300 K, using a stochastic constant temperature algorithm with heat bath particles of relative mass 0.01 (Kast et al., 1994). After equilibration, microcanonical trajectories were calculated for 5-ps intervals, separated by short equilibration periods of 0.25 ps, and each trajectory was followed in this way for a total of 1 ns. Newton's equations of motion were integrated with a standard Verlet algorithm using a timestep of 0.5 fs. Periodic boundary conditions were applied along the z direction of each channel axis, and all interactions were cut off at a length of 12 Å with a shifted force potential. For a discussion of the long-range interactions that may be appropriate in the model channels considered here we refer the reader to work of LBR (Lynden-Bell and Rasaiah, 1996). The short equilibration periods between microcanonical sections were also performed with the constant temperature algorithm at 300 K, but with lighter heat bath particles of relative mass 0.005. This mixed microcanonical/canonical approach was taken to overcome slight energy drifts encountered during the microcanonical sections of the trajectory, as well as to ensure better sampling of the phase space explored by each trajectory. The use of a small relative mass for the bath particles ensures only a minimal deviation from Newtonian dynamics and does not affect the calculated diffusion constants. This was confirmed in the case of the narrowest channel by computing the CEC diffusion constant from 100 microcanonical trajectories (each with a shorter simulation time of 0.25 ns, equilibrated initial conditions, and temperature of ~300 K), which was found to be consistent with results obtained from the microcanonical/canonical procedure described above.

Quantities such as the atomic positions, the CEC coordinates, and the amplitudes of the EVB states were sampled at 25-fs intervals during the simulations so that various probability distributions P_x and their associated free energies F_x = k_BT ln[P_x] could be calculated. Data gathered from 10 trajectories as described above were used to evaluate the static properties of each channel. Mean-squared displacements of both the oxygen atoms and the CEC were also monitored during the simulations for the calculation of the associated diffusion constants, and further details of these calculations will be given at the appropriate place in the text. At various stages throughout the Results and Discussion section we refer to bulk phase results for an excess proton in liquid water. These have been obtained from two sources, the original study by Schmitt and Voth (1999a), and also some more recent MS-EVB simulations (Cuma et al., to be published).

To investigate possible effects arising from finite channel length, channels of radius 2.0 and 3.0 Å were also considered with approximately twice the length (l  60 Å) and number of water molecules. Properties of these two elongated channels are shown at the bottom of Table 1, and results obtained with them are given in the Appendix.

	RESULTS AND DISCUSSION

TOP ABSTRACT INTRODUCTION METHODOLOGY RESULTS AND DISCUSSION CONCLUSIONS APPENDIX REFERENCES

Static properties

The pivot oxygen-oxygen (O*O) and oxygen-oxygen (OO) radial distribution functions obtained from channel models with different radii are shown in Fig. 1 (where the pivot oxygen O* belongs to the pivot hydronium). In conjunction with this figure, O* and O coordination numbers are plotted as a function of channel radius in Fig. 2. It is clear from Fig. 1 that, as in bulk, the O*O equilibrium separations are shifted to shorter distances than the OO separations inside the model channels.

View larger version (17K): [in this window] [in a new window]   FIGURE 1   The pivot oxygen-oxygen O*O and oxygen-oxygen OO radial distribution functions g(r) obtained from the model channels of radius rcyl = 2, 3, 4, and 5 Å. The O*O distribution function is shown with a solid line, and the OO distribution function with a dotted line.

View larger version (17K): [in this window] [in a new window]   FIGURE 2   Coordination numbers and the number of donor hydrogen bonds per water molecule as a function of channel radius rcyl. Pivot oxygen and non-pivot coordination numbers were obtained by integrating their associated radial distribution function to 2.9 and 3.3 Å, respectively. A hydrogen bond was designated if the OO distance was <3.5 Å, and the angle between the OO and OH bond vectors was <30°.

Considering now the O*O distribution functions, it is seen that two different first solvation shell O*O equilibrium distances are evident in the channels. In the narrow channel of radius 2.0 Å the O*O distribution function displays a single peak at 2.41 Å. This peak is seen as a shoulder in the distribution function obtained from the channel of radius 3.0 Å along with another, more dominant, peak at 2.54 Å. In general, the peak at 2.54 Å becomes more dominant as channel radius increases, while the peak at 2.41 Å becomes less dominant. Gas-phase MS-EVB calculations at global minimum energy geometries predict an O*O distance of 2.40 Å in the Zundel cation H₅O, and a longer O*O distance of 2.57 Å in the Eigen cation H₉O. Hence, it appears that solvated forms of these two cations are to some extent responsible for the two O*O distances observed in the channels, with the O*O distance associated with the Eigen species becoming more apparent as channel radius increases, i.e., when there is sufficient space for the larger of the two cations to be formed. This idea is supported by the observed O* coordination number, which is shown as a function of channel radius in Fig. 2. The O* coordination number has a value of 2 in the channel of radius 2.0 Å and a value of 3 in all other channels. In the channel of radius 5.0 Å the first solvation shell peak is at 2.53 Å, which is close to the same peak in bulk calculations at 2.52 Å and reflects the ensemble average of all O* distances that contribute to the distribution function.

The OO radial distribution functions, which are obtained by omitting the pivot oxygen atom from the calculation and therefore reflect primarily water structure and not hydronium solvation, all display a single first solvation shell peak. The position of this peak shifts slightly from ~2.68 Å in the 2.0-Å-radius channel to 2.73 Å in the largest 5.0-Å-radius channel. This is because in the narrow channel a larger fraction of the total number of oxygen atoms will be in closer proximity to the charged species, and as a result will be more hydronium in character and have shorter OO separations. As channel radius is increased and the number of water molecules becomes larger, this fraction decreases and the OO equilibrium bond separation approaches the bulk phase value of 2.73 Å. This value is slightly shorter than the reported bulk phase value for TIP3P water of 2.78 Å because only the pivot oxygen was omitted in these distribution function calculations; neglecting all EVB oxygen nuclei does recover this value (Schmitt and Voth, 1999a). The coordination number of non-pivot oxygen atoms also has a value of 2 in the 2.0-Å-radius channel and this number increases with increasing channel radius to reach a value of 4.3 in the 5.0-Å-radius channel which is closer to the bulk coordination number of 4.5.

Cross-channel density profiles ((r) = P_r/(2r) where r is the perpendicular distance from the channel axis) of the oxygen and hydrogen atoms, and the CEC in various channels, are shown in Fig. 3. In all channels the oxygen atoms, hence also the water molecules, form distinct solvation layers at a distance of ~1 Å from the channel wall. At larger channel radii another solvation layer is able to form in the middle of the channel; this begins in the 4.5-Å-radius channel and can be seen in the density profile of the 5.0-Å-radius channel. Hydrogen atoms, which do not feel the repulsive channel wall directly, can protrude slightly outside the channel walls and generally follow the O atom distributions, except in the case of the narrow 2.0-Å-radius channel where hydrogen atom density is at the center of the channel. Because hydrogen atoms extend in this manner they will not be involved in hydrogen bonding. Indeed, the number of donor hydrogen bonds per water molecule in the channels, which is also shown in Fig. 2, is reduced relative to the bulk value of ~1.7 but increases with increasing channel radius. Similar oxygen and hydrogen density profiles were observed by LBR in channel simulations of water modeled with the rigid SPC/E potential (Lynden-Bell and Rasaiah, 1996), suggesting that neither OH bond flexibility nor the presence of the excess proton alters these properties significantly. (In the 4.0-Å-radius channel, however, we do not observe the anomalously low number of hydrogen bonds or large number of close OH pair interactions as seen by LBR, and this could be an artifact of the different water models.)

View larger version (21K): [in this window] [in a new window]   FIGURE 3   The cross-pore density profiles (r) obtained from model channels of radius 2, 3, 4, and 5 Å. The channel radius is indicated in the top right-hand corner of each panel. Oxygen atom density is shown with a solid line, hydrogen atom density with a dotted line, and the center of excess charge density with a dashed line.

The CEC density profile is similar to the H atom density profile, except in cases where a central solvation layer is formed. In these cases the CEC displays a distinct preference to remain within the solvation shell that is closest to the channel wall. LBR analyzed this type of behavior for a variety of charged and neutral solute particles inside the 3.0-Å-radius channel in terms of their energetic and entropic contributions to Landau solvation free energies (Lynden-Bell and Rasaiah, 1996). It was found that "small ions favor the channel center because the energetic terms dominate the Landau free energy, and large ions favor the outside, as for them entropic terms dominate," i.e., there is a favorable energetic contribution to the free energy as an ion is solvated, and this will be most favorable at the channel center where complete solvation is possible; however, at the same time there will be an entropy cost due to increased ordering of solvent water molecules. The balance between these two terms determines where a solute ion prefers to occupy the channel. A bare proton is clearly a very small and highly charged localized ion, and based on the above argument might be expected to reside near the channel center. However, due to the actual delocalization of the excess protonic charge over several water molecules, the migrating species can be considered as effectively a rather large ion with a diffuse charge distribution. This would therefore favor the outside edge of the channels, which is indeed the case in the present simulations of the 3.0-Å-radius channel. It is reasonable to expect that similar arguments apply to the observed CEC density profiles of the other channels.

The relative stabilities of aqueous solvation structures of the excess proton, in particular the Zundel and Eigen complexes, are important and timely topics in aqueous PT (Agmon, 1999). The picture from classical MS-EVB simulations of an excess proton in liquid water is one where the solvated Zundel and Eigen species are to be regarded as limiting structures, with the Eigen species being the slightly more stable of the two by ~1.1 kcal/mol, depending on the means by which it is estimated (Schmitt and Voth, 1999a). To study the solvation structures inside the channels we have considered the two-dimensional probability distributions P_q,rOO, where r_OO is the shortest O*O distance, q is the associated PT coordinate, defined as

(4)

and r_OH is the distance between O* and the transferring hydrogen. Free energy surfaces obtained from this distribution function for certain channels are shown in Fig. 4. The surface for the 2.0-Å-radius channel displays a single low-energy basin centered at q = 0.0 Å and r_OO = 2.38 Å, which corresponds to the solvated Zundel complex. At a channel radius of 3.0 Å the basin of lowest free energy extends symmetrically to regions where r_OO = 2.45 and q = ±0.15 Å, which correspond more closely to solvated forms of the Eigen complex. [Note that the lowest energy basins, i.e., the most heavily shaded regions in Fig. 4, depict configurations with a free energy in the range 0-0.25 kcal/mol, which are therefore easily mixed by thermal fluctuations at 300 K.] Once the channel radius has reached 4.0 Å the two-dimensional free energy distribution has become bimodal, with distinct low-energy basins centered at r_OO = 2.47 and q = ±0.17 Å. The free energy surfaces for the channels with radii 4.0 and 5.0 Å are very similar to the classical free energy surfaces obtained for the excess proton in bulk water (Schmitt and Voth, 1999a). This suggests that, at least in terms of these two-dimensional probability distribution functions, an ~4.0-Å-radius hydrophobic channel is sufficient to provide a bulk water-like environment for the microscopic solvation of the excess proton.

View larger version (71K): [in this window] [in a new window]   FIGURE 4   The free energy profiles obtained from the two-dimensional probability distribution PrOO,q for channels of radius 2, 3, 4, and 5 Å. The channel radius is indicated in the bottom left-hand corner of each panel. The free energy contours are shown at 0.25 kcal/mol intervals and the heaviest shading indicates regions of lowest free energy.

MS-EVB amplitude analysis

The probability distribution and value of the dominant EVB state amplitude c also provides a convenient way of looking at the proportion of the Eigen and Zundel entities in the liquid phase (Schmitt and Voth, 1999a). In bulk a c value of ~0.7 is found to correspond to a symmetrically solvated Eigen complex with three approximately equal O*O distances in the range 2.5-2.6 Å. A c value of ~0.5, where two EVB states share most of the total amplitude equally, corresponds to a solvated Zundel complex with a short O*O distance of ~2.4 Å. Clearly, if the channel wall impinges upon the solvated proton structure, then these interactions will modify the structures to which these c values correspond, as will be discussed below. The free energy curves obtained from the probability distributions given by F = k_BT ln are shown in Fig. 5 for a representative set of the model channels.

View larger version (22K): [in this window] [in a new window]   FIGURE 5   The free energy profile obtained from the probability distribution of the maximum EVB amplitudes Pcmax in a representative selection of channels. Each line is labeled with the radius of the channel in angstroms.

Consider first the narrow 2.0-Å-radius channel. Due to the single-file arrangement of molecules (recall that the pivot and non-pivot oxygen coordination numbers are both equal to 2) it is not possible for the Eigen cation to be formed in this narrow channel. The values of c can therefore be interpreted in terms of the following two limiting situations. First, at c ~0.5, amplitude is generally distributed evenly between two EVB states of hydronium cations that share a hydrogen atom. This situation corresponds to the formation of a Zundel-like complex with a short O*O distance of ~2.4 Å or, equivalently, a strong hydrogen bond. Second, at values of c  0.6, the remaining amplitude is distributed between the two EVB states that have hydronium ions on either side of the pivot hydronium, with slightly longer O*O separations of ~2.5 Å that can be considered as medium-strength hydrogen bonds. In some sense this second structure could be regarded as an Eigen cation that has one of its three first solvation shell water molecules removed. It is clear from Fig. 5 that the second of these two situations is more probable than the first, but at the same time that the free energy difference between the two situations is very small at 0.2 kcal/mol.

In 2.5-Å-radius channels there is sufficient space for threefold water solvation of the pivot hydronium to occur and the values of c can be interpreted as described earlier for the bulk liquid case. The free energy at c ~0.5, i.e., that for formation of the Zundel complex, increases with channel radius, and when it reaches a value of ~0.9 kcal/mol in the 5.0-Å-radius channel it is more comparable to the bulk phase value of 1.15 kcal/mol for this distribution function. The free energy profile for the 2.5-Å-radius channel is slightly anomalous inasmuch as the Zundel formation free energy is higher than that of the 3.0-Å-radius channel.

In general, the various c free energy curves show that the Zundel cation can be stabilized by almost 1.0 kcal/mol when compared to bulk, depending on the radius of the channel. It is interesting to note, however, that only in the case of the 2.0-Å-radius channel is the free energy of the Zundel cation significantly less than k_BT at 300 K (0.6 kcal/mol).

We now consider the free energy curves associated with the probability distribution of a PT "reaction coordinate" (Schmitt and Voth, 1999a), which is defined as

(5)

where the subscripts k and j denote the two EVB states with largest amplitude. The coordinate defined in Eq. 5 will be zero when the two largest amplitudes c and c are equal, e.g., in the case of the Zundel complex, and will take on a finite value for non-equal values of the two amplitudes. The free energy curves obtained from the probability distribution P_qreact in various model channels are shown in Fig. 6. The free energy for the formation of the Zundel cation intermediate is again lowest in the 2.0-Å-radius channel. In the 5.0-Å-radius channel the free energy barrier is ~1.0 kcal/mol, which is more comparable to the classical bulk phase value of 1.25 kcal/mol for the P_qreact distribution. Again, the free energy curve for the 2.5-Å-radius channel is qualitatively different from the other curves. In this case the Zundel cation is metastable with respect to non-symmetric distributions of the two largest EVB amplitudes.

View larger version (23K): [in this window] [in a new window]   FIGURE 6   The free energy profile obtained from the probability distribution of the proton transfer reaction coordinate Pqreact in a representative selection of channels. Each line is labeled with the radius of the channel in angstroms.

To investigate the origin of the P_qreact free energy profile in the 2.5-Å-radius channel we have visualized the MD trajectories. An example of the total configuration in the channel is illustrated in Fig. 7. The structure shown in this figure contains a Zundel cation that is oriented almost perpendicular to the channel axis and forms an integral part of two five-membered hydrogen-bonded water molecule rings. Structures similar to the one in Fig. 7 are quite persistent in the 2.5-Å-radius channel simulations and might help to explain the metastability with respect to the Zundel formation seen in Fig. 6. Although this is by no means a quantitative explanation, it is interesting to note that a similar configuration to that shown in Fig. 7 has been put forward by Agmon as a possible stable structural unit for concentrated HCl/water mixtures, where the water molecules at the top and bottom of the five-membered rings are replaced with Cl anions (Agmon, 1998). Furthermore, it appears that stable microscopic solvation structures such as the one depicted in Fig. 7 are correlated with the dynamical properties of the CEC, as we will discuss in the following section in the context of CEC diffusion along the channel axis.

View larger version (13K): [in this window] [in a new window]   FIGURE 7   An example configuration taken from a model channel with radius 2.5 Å. A solvated Zundel cation H5O is indicated at the center of the configuration, lying perpendicular to the axis of the channel, which is parallel to the z axis.

Channel diffusion

To probe dynamical properties of the model channels one can consider the diffusion of both oxygen atoms and the CEC. Diffusion constants for these two species were determined using the Einstein relation, which for a particle in one dimension with coordinate x is defined as

(6)

where D is the one-dimensional diffusion constant, x(t) is the position of the particle(s) under consideration, and the brackets ( ) denote the ensemble average (Allen and Tildesley, 1987).

Due to the "hopping" nature of a migrating excess proton in water, a water molecule is not uniquely defined throughout a simulation, i.e., it may at times become a more hydronium-like species and also exchange the oxygen atoms to which it is chemically bonded. Hence, instead of using the molecular centers of mass to estimate water self-diffusion constants, we use the coordinates of the oxygen atoms, since both approaches will lead to the same result in a simulation of pure water. The mean-squared displacements of the CEC and oxygen atoms were monitored, in the x and y directions (i.e., perpendicular to the channel axis) and the z direction (i.e., parallel to the channel axis) for up to 50 ps, with origins taken at 25-fs intervals. For a given model channel, a total of 10 trajectories, each with different initial conditions and a total simulation time of 1 ns, were used to evaluate the mean-squared displacement of the oxygen atoms. An additional 10 such trajectories were used to compute the CEC mean-squared displacement. The CEC behaves like a single particle in the system and longer simulation times were required to obtain satisfactory statistics. Finally, diffusion constants were extracted by fitting a straight line to the mean-squared displacements in the range 25-50 ps. The results of these calculations for diffusion along the channel axis are displayed in Fig. 8.

View larger version (14K): [in this window] [in a new window]   FIGURE 8   Diffusion constants of (a) the oxygen atoms DOXY and (b) the center of excess charge DCEC as a function of the channel radius rcyl. Note that the ordinate of each graph has a different scale. The DOXY error bars are the standard deviation of the diffusion constants obtained from the 10 (1-ns) trajectories. The DCEC error bars are the standard deviations calculated from four diffusion constants, each of which was obtained by dividing the 20 (1-ns) trajectories into four groups of five and calculating the diffusion constant for each group.

Consider first the axial diffusion constants of water oxygen atoms inside the channels shown in Fig. 8 a. In general, the diffusion constant of the oxygen atoms is reduced relative to the bulk phase value of 0.34 ± 0.04 Å²/ps. The diffusion constant of oxygen atoms in the 2.0-Å-radius channel is in fact too small to be accurately measured on the timescale of the present simulations; we estimate that the oxygen diffusion constant in this channel is of the order of 10⁴ Ų/ps. The water molecules in this narrow channel are restricted such that they are not able to pass around one another and the only mechanism for water diffusion is via the collective motion of all molecules along the channel (Levitt, 1984). In the 2.5-Å-radius channel it becomes possible for water molecules to pass one another so there is measurable diffusion in the axial direction. The oxygen diffusion constant in the 4.0-Å-radius channel is lowered relative to the 3.5- and 4.5-Å-radius channels; this result is in agreement with the pure water simulations of LBR, who concluded that water structure inside this channel was anomalously stable (Lynden-Bell and Rasaiah, 1996).

Turning now to the axial diffusion constants of the CEC that describe PT along the model channels, the most striking result in Fig. 8 is that in the narrow 2.0-Å-radius channel the CEC diffusion constant is greater, by an order of magnitude, than the CEC diffusion constants obtained from the other channels. The diffusion constant is lowest in the 2.5-Å-radius channel, where it is approximately equal to the diffusion constant of oxygen atoms in the same channel. In the 3.0-5.0-Å range of channel radii there is relatively little variation of the CEC diffusion constants. Except for the CEC diffusion constant in the 2.0-Å-radius channel, all CEC diffusion constants are reduced relative to the classical bulk value of 0.45 ± 0.11 Å²/ps, but they are generally all greater than the water-oxygen atom diffusion constants in the same channels.

We now propose an explanation for these CEC diffusion results, particularly the dramatically increased diffusion constant in the narrow 2.0-Å-radius channel. In this channel each water molecule participates in two hydrogen bonds, one with each of the two nearest-neighbor water molecules that lie above and below it in the channel. In addition, the pivot hydronuim also participates in two hydrogen bonds with its nearest neighbors. As outlined already in terms of limiting c values in this channel, the pivot hydronium either forms one strong hydrogen bond to form a Zundel complex with one water neighboring water molecule and a weaker hydrogen bond with the other, or two hydrogen bonds of medium strength, one with each of its nearest-neighbor water molecules. The net effect is that the channel contains a continuous one-dimensional chain (or "wire") of hydrogen-bonded molecules arranged along its axis. It is clear from Fig. 5 that the configurations that correspond to the two limiting c values inside the channel can interconvert very readily due to thermal fluctuations and, from Eq. 3, that such interconversions will induce rapid CEC motions along the channel. When channel radius is increased beyond 2.0 Å it becomes possible for the Eigen complex to be formed, and as a result the distribution of directions along which PT might occur becomes more isotropic. The direction of PT is no longer exclusively "steered" along the channel axis, and the free energy to pass through the Zundel intermediate in the hopping process is increased relative to the thermal energy. Furthermore, certain configurations in the wider channels appear to trap the CEC in the channel. For example, in the configuration shown in Fig. 7 the direction of PT is perpendicular to the channel axis, so that while the system is in this configuration the CEC cannot displace significantly in the z direction, and there is a reduction in the axial CEC diffusion constant. It is clear that solvation structures, which are determined by the width of the channel, play an important part in determining the mobility of the CEC and water molecules along the channel. Indeed, the existence of a true "proton wire" in which proton diffusion is enhanced by an order of magnitude depends critically on the geometrical constraints and resulting solvation structure. In the present case, an increase in channel radius of just 0.5 Å served to disrupt the proton wire function.

Further discussion

Before presenting conclusions that can be drawn from this study it is appropriate to discuss the limitations and significance of our results.

It has been shown that quantum mechanical effects may be non-negligible in the equilibrium and dynamic properties of the excess proton in bulk water (Schmitt and Voth, 1999a). Compared to classical results, Schmitt and Voth found that quantization of the nuclear degrees lowered the free energy barriers to Zundel formation by ~35%, and that the Eigen and Zundel structures should be regarded as limiting structures that are strongly mixed by quantum and thermal fluctuations at 300 K in the liquid phase. It is reasonable to expect that similar reductions in barrier height would be found upon quantization of the nuclei in the model channels considered here, also leading to a physical picture where limiting solvation structures become less well-defined and more easily mixed. At the level of dynamics, proton transfer rates were found to be enhanced by a factor of approximately two upon quantization of the nuclei, and a similar enhancement of PT rates can also be expected in the model channels, which would increase the CEC diffusion constants.

For simplicity we have used a smooth cylindrical hydrophobic channel to confine the excess proton and water molecules. In support of these models, we note that the relative reduction in the mobility of sodium ions inside them is in agreement with the reduction that was observed in some more realistic models for hydrophobic channels (Smith and Sansom, 1998). However, the insides of ion channels are often lined with hydrophilic residues with which water molecules can hydrogen-bond and reduce the extent of water-water hydrogen bonding inside the channel. In turn, these water-channel interactions would affect the mobility of the CEC. This process has been demonstrated by Pomès and Roux, who in neglecting water-channel electrostatic interactions during classical simulations of a protonated water wire in the gramacidin channel using the PM6 potential found an increase in the mobility of a hydronium coordinate (Pomès and Roux, 1996a). Clearly, the extent to which hydrophilic residues inside a channel will affect the CEC and water molecule diffusion will depend on the specific structure of the channel. Nevertheless, the main conclusions of this study concerning the effects of pore radius and proton solvation structure on the CEC diffusion are expected to hold qualitatively.

As stated in the Introduction, the aqueous PT results presented here, particularly the enhanced diffusion constant of the CEC when a water wire is formed, are relevant to a proposed mechanism for the permeation of protons through lipid bilayers (i.e., a hydrophobic environment), which is via the formation of transient hydrogen-bonded chains of water molecules (Nichols and Deamer, 1980; Nagle, 1987; Deamer, 1987; Marrink et al., 1996). If one assumes that the diffusion constant of a small cation or anion such as a sodium ion will be of the same order of magnitude as the water molecules in the narrow channel (10⁴ Ų/ps), then the CEC diffusion constant (3.6 Å²/ps) will exceed it by five orders of magnitude. This large difference is a direct consequence of the chemical, as opposed to hydrodynamic, mechanism by which the CEC can diffuse. The present results imply that in the narrow 2.0-Å-radius channel a period of water wire stability of ~100 ps is needed for the CEC to acquire a root-mean-squared displacement of 30 Å along the hydrogen-bonded water molecules in its interior.

This work has also highlighted the benefits of the molecular dynamics approach. Physical properties such as the cross-pore density profiles, the variation of microscopic solvation structures with channel radius and, in particular, their effect on the diffusion of the CEC would unlikely be captured with a continuum approachthese observables are a direct result of treating molecular interactions explicitly with an appropriate potential energy surface. One interesting possibility though, as suggested by Smith and Sansom (1998), is the incorporation of data generated from MD simulations (e.g., diffusion constants) into continuum approaches such as those based on, for example, the Poisson-Nernst-Planck equation (Cárdenas et al., 2000), to allow simulations over longer timescales than are currently possible with the MD method.

	CONCLUSIONS

TOP ABSTRACT INTRODUCTION METHODOLOGY RESULTS AND DISCUSSION CONCLUSIONS APPENDIX REFERENCES

In this article we have presented a classical molecular dynamics study of aqueous PT in smooth hydrophobic cylindrical channels with a range of different radii. The channel radius has been found to play a critical role in both the equilibrium and dynamical properties of the excess proton in these systems.

As one might have anticipated simply on geometric grounds, the Zundel complex, H₅O, is increasingly stabilized relative to the Eigen complex, H₉O, as channel radius is decreased. This has been confirmed by the calculation of radial distribution functions and various probability distributions for channels with different radii. Our calculations also demonstrate that these solvation structures for the aqueous excess proton, as in bulk, should also be regarded as limiting structures inside the model channels and that they are readily mixed by thermal fluctuations.

Water structure inside the model channels is found to be organized into distinct solvation layers, as also observed in pure water channel simulations (Lynden-Bell and Rasaiah, 1996). The CEC, a coordinate that reflects the electronically delocalized nature of the charge due to the excess proton, is found to preferentially occupy the water solvation shell that is closest to the channel wall. We assert that this behavior is due to a decrease in entropy that would otherwise occur if the excess charge were to adopt a more central position in the channels.

The diffusion of water molecules and the CEC along the channels display a sensitive dependence on the radius of the channels. In the narrow 2.0-Å-radius channel the water molecules and excess proton adopt a hydrogen-bonded wire-like configuration in which the diffusion constant of water is very small, but the diffusion constant of the CEC is much larger than it is in channels of greater radii. This increased diffusion constant or conductivity of the channel is a combined result of the continuous one-dimensional nature of the hydrogen-bonding network along the narrow channel, and also the low free energy barrier to proton transfer inside the channel. As soon as channels become wide enough for threefold water solvation of the hydronium ion to occur the situation becomes more complicatedthe solvated proton complex can become more stabilized and proton transfer no longer occurs exclusively along the direction of the channel axis.

The simulation results presented in this work are intended to display some general features of the equilibrium and dynamical properties of the aqueous excess proton as a function of channel radius. At the same time it will also be desirable to study these properties in more realistic biological proton channels, and such work is currently underway in our research group.