Search tips
Search criteria 


Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
Phys Chem Chem Phys. Author manuscript; available in PMC 2012 October 14.
Published in final edited form as:
PMCID: PMC3470879

Energetics and dynamics of proton transfer reactions along short water wires


Proton transfer (pT) reactions in biochemical processes are often mediated by chains of hydrogen-bonded water molecules. We use hybrid density functional calculations to study pT along quasi one-dimensional water arrays that connect an imidazolium-imidazole proton donor-acceptor pair. We characterize the structures of intermediates and transition states, the energetics, and the dynamics of the pT reactions, including vibrational contributions to kinetic isotope effects. In molecular dynamics simulations of pT transition paths, we find that for short water chains with four water molecules, the pT reactions are semi-concerted. The formation of a high-energy hydronium intermediate next to the proton-donating group is avoided by a simultaneous transfer of a proton from the donor to the first water molecule, and from the first water molecule into the water chain. Lowering the dielectric constant of the environment and increasing the water chain length both reduce the barrier for pT. We study the effect of the driving force on the energetics of the pT reaction by changing the proton affinity of the donor and acceptor groups through halogen and methyl substitutions. We find that the barrier of the pT reaction depends linearly on the proton affinity of the donor but is nearly independent of the proton affinity of the acceptor, corresponding to Brønsted slopes of one and zero, respectively.

1. Introduction

Proton transfer (pT) is an elementary chemical reaction central to enzyme catalysis and membrane transport.15 Anticipating a molecular picture of pT in water, Grotthuss6,7 suggested that it is a charge defect, rather than the proton itself, that diffuses. In the charge-defect diffusion model, an excess proton associated with a water molecule carrying the defect, changes its bonding partner, thus reorienting the polarity of the hydrogen bond and moving the location of the defect. The basic mechanisms of pT processes have been studied extensively by molecular simulations5,816 and experiments.1,17,18 Water-mediated pT in solution and in the confined environments of channels and enzyme active sites is associated with transitions between different hydration structures of the excess proton, with hydronium (H3O+), Eigen17 (H9O4+), and Zundel18 ions (H4O2+) as extreme forms.816

In contrast to electron transfer (eT), where quantum mechanical tunnelling is a major transfer mechanism,1922 the relatively large mass of the proton limits tunnelling contributions at ambient conditions except for short donor-acceptor distances.2327 In enzymes, where the distance between the initial proton donor and the ultimate acceptor can reach ~10–30 Å, a proton conduction pathway must thus be employed. Such pathways typically involve water molecules. In narrow and weakly polar protein channels, water tends to form quasi one-dimensional hydrogen bonded chains that are thought to provide excellent proton conduction pathways.2832

Water arrays similar to the ones observed inside enzymes, also form within carbon nanotubes,2833 mimicking the low-polarity environment of the protein interior. Accordingly, nanotubes have served as important model systems to explore the physical principles of biological pT. Studies of water-mediated pT in nanotubes and similar systems have demonstrated the importance of electrostatic effects15,16,2932,34,35 in forming the water chain that connects the proton donor and acceptor, and in determining the height of the proton dissociation barrier and the charge mobility.

Here we use quantum chemical density functional theory to study the pT between a proton donor and an acceptor connected by a short water chain. The calculations are performed in vacuum, and in a low-dielectric environment, to capture the essence of water-mediated pT reactions in the low-polarity protein interior. The protons are treated as classical particles, thus ignoring nuclear quantum delocalization effects.36 By halogenation and methylation of one of the imidazole substituents, we control the driving force for the pT and correlate it to the height of the energy barrier. Dynamic effects are studied by constructing transition paths on the Born-Oppenheimer surface, obtained by classical molecular dynamics simulations of the nuclear positions on the quantum mechanical energy surface.

2. Computational models and methods

Quasi one-dimensional proton wires were built by connecting a proton donor (D) and acceptor (A) with N=1–4 intervening water molecules (Figure 1 top), so that the average distance between the heavy atoms in the chain was similar to that observed in molecular dynamics (MD) simulations of water molecules in carbon nanotubes.28,29 The distances between the terminal Cβ-atoms were (2.52N + 10.28) Å for chains of N water molecules. In a symmetric arrangement that mimics pT between two histidine side chains, 4-methyl-imidazolium and its conjugate base 4-methyl imidazole were used as proton donor and acceptor, respectively. To study the effect of an increase in the driving force of the pT reaction, additional calculations were performed in which the terminal methyl-Cβ-carbon atoms of the donor group were chlorinated, fluorinated, and methylated.

Figure 1
(Top) Structure of proton wire with N=4 water molecules. (Bottom) Difference in potential energy to reactant state as a function of dtot, the average distance between donor-hydrogen pairs in the reactant structure. Energies are shown for idealized hydronium ...

To explore the energetics of the pT reaction, and to determine transition state structures, geometry optimizations were performed at the B3LYP/def2-SVP level of theory.3739 In the optimizations, the distance between the terminal Cβ-carbon atoms was fixed. In a protein, this constraint would correspond to freezing the fluctuating backbone. However, since we find that proton transfer is fast (~200 fs; see below) on the time scale of backbone motions, we expect the effect on the pT dynamics to be small.

As an intermediate step following the optimization of reactant and product state structures, we optimized the structures of idealized hydronium (H3O+) and Zundel cations (H5O2+) at different positions along the pT pathway. The three d(O-H) distances of the idealized hydronium cation were constrained to 0.9854 Å, and the two central d(O-H) distances of the Zundel cation were constrained to 1.20 Å, with d(O-O)=2.40 Å. Starting from the optimized structures with idealized Zundel ions, the constraints on the d(O-H) distances were released and transition state optimizations were performed.

To identify dynamically important species of the pT reactions, transition paths were constructed between reactant and product states using Born-Oppenheimer molecular dynamics simulations. Trajectory pairs were initiated from the transition state structures, and the classical equations of motion on the quantum mechanical energy surface at the B3LYP/def2-SVP level of theory were integrated until either the reactant state or the product state was reached, respectively. The leapfrog algorithm was used in these simulations with an integration time step of 0.5 fs. Initial particle velocities were assigned from a Maxwell-Boltzmann distribution at a temperature of 310 K, with velocity vectors of the same magnitude but inverted signs used in the trajectory pairs. The two trajectory segments can thus be thought of as evolving forward and backward in time, respectively.40 Trajectory pairs that end in opposite states can be stitched together to form a transition path after inverting the time order and velocities of the segment ending in the reactant state.41

Before the final analysis, single-point calculations were performed at the B3LYP/def2-TZVP level of theory42 on all structures to extract energies and electrostatic potential charges. The accuracy of the B3LYP calculations were assessed by comparison to spin-component scaled Møller-Plesset perturbation theory (SCS-MP2) calculations.43 For reference, Table 1 also lists results for Hartree-Fock (HF) calculations, MP2 second-order approximate coupled cluster theory (CC2),44 and density functional theory with the following functionals: Becke-Parr 86 (BP86),45,46 Tao-Perdew-Staroverov-Scuseria (TPSS),47 and dispersion corrected B3LYP (D-B3LYP).48,49 The results in Table 1 suggest that the hybrid density functional B3LYP accurately describes the barriers for the pT reaction, with a deviation of less than 0.4 kcal mol−1 relative to ab initio SCS-MP2 theory, whose estimated accuracy is expected to be comparable to that of coupled cluster theory with single and double excitations, CCSD.43 HF theory severely overestimates the reaction barrier, whereas the gradient generalized functional (GGA) BP86 and the meta-GGA TPSS somewhat underestimate the barrier. Dispersion effects obtained by the dispersion corrected B3LYP (D-B3LYP) are found to be small. Accordingly, we consider only pure B3LYP energies in the remaining study.

Table 1
Assessment of computational methodology on the pT barrier ΔE* with one water molecule.

During method testing, we observed that the use of sufficiently large basis sets of triple-ζ quality was crucial for an accurate energetic description of the protonated intermediates. SVP-quality basis sets (comparable to 6-31G*) led to underestimation of protonated intermediates by up to ~5 kcal mol−1.

Dielectric effects were studied using the conductor-like screening model (COSMO).50 In these calculations, the system is embedded in a molecule-shaped cavity determined for a probe radius of 1.3 Å within an otherwise homogeneous medium of dielectric constants ε=4 and 80, respectively.

Vibrational contributions to kinetic isotope effects (KIE) were estimated by calculating zero-point energies (ZPE) of reactant and transition state species. We evaluated the molecular Hessian, defined as the second derivative matrix of the energy with respect to displacements in the coordinates, at the B3LYP/def2-TZVP level of theory to compute ZPE and entropic contributions to the activation free energy for the systems with N=1–4 water molecules. KIEs were obtained from the difference of reactant and transition state ZPEs with hydrogen and deuterium masses,5153


where ΔZPED* and ΔZPEH* are the differences in the ZPEs at transition and reactant states for the deuterated and hydrogenated systems, respectively, and 1/β=kBT ~ 0.6 kcal mol−1 at 310 K. Additional tunnelling corrections are estimated from Wigner’s semi-classical approximation involving the imaginary frequencies (iω*) at the transition state by multiplying kH/kD in Eqn. 1 with54


However, one should keep in mind the approximate nature of this KIE estimate, without a fully quantized proton and without quantum dynamics.36

All structure optimizations and single-point energy calculations were performed using TURBOMOLE 6.0 and 6.2.55 VMD57 was used for analysis and visualization of molecular structures. ORCA56 was used for computing ZPEs. To calculate reaction free energies ΔG*, we included contributions from translational, rotational, and vibrational entropies and enthalpies, evaluated within harmonic approximations.56

To quantify the progress of the reaction, we defined dtot as the average distance di of the hydrogen atoms involved in hydrogen bonds in the reactant proton wire and their acceptor atoms N (i=0) and O (i=1 to N),


During the pT process the individual di increase as the initially covalent bond transitions into a hydrogen bond. We found this coordinate to resolve the reaction dynamics well on the time scales of individual pT events, despite the evident shortcomings at longer times because of changes in the identity of the protons involved in hydrogen bonds. In previous studies of pT, also other reaction coordinates had been used, e.g., the dipole moment of the water chain13,14 and the centre of excess charge.31,33

In addition to visual inspection, Zundel, hydronium and water species were identified from the trajectories by computing the average over the distances rij of oxygen atom i to its three nearest hydrogen atoms,


with rij=d(Oi−Hj) and ri1< ri2<…. The average distance ravei increases from ~1.05 Å for hydronium ions to ~1.08 Å for Zundel ions and 1.16 Å for water.

3. Results and discussion

3.1 Energetics

Figure 1 compares the potential energies of idealized hydronium and Zundel-ions obtained by constrained minimization to the energies of fully optimized transition state structures. For a chain with N=4 water molecules, we considered four possible hydronium ions and three Zundel ions. Three transition states were found by saddle-point searches. Energies relative to the optimized reactant state are plotted as a function of the progress coordinate dtot. Ideal hydronium ions are 3–4 kcal mol−1 above the respective Zundel ions, and thus of limited relevance at ambient conditions. In contrast, ideal Zundel ions are energetically close to the optimized transition states.

Figure 2 shows the structures of the symmetric transition states for chains with N=1 to 4 water molecules. For chains with even numbers of water molecules, N=2 and 4 (Figure 2B and E), this central transition state has Zundel-like character, with a proton at the centre of the chain shared by the two flanking water molecules at a location halfway between the deprotonated donor and acceptor imidazoles. In contrast, for odd-numbered chains, N=1 and 3 (Figure 2A and C), the central transition state has a hydronium-like character. Also shown is the second transition state in the chain with N=4 water molecules, in which a hydronium-like species is between water molecules 3 and 4 (or, by symmetry, 1 and 2). The small deviations from perfectly symmetric bond distances in Fig. 2A–C and E indicate numerical uncertainties in the globally optimized transition state structures.

Figure 2
Transition state structures and energies relative to reactant state with a protonated imidazole for (A–D) N=1 to 4 water molecules, and (E) the second transition state for N=4. Electrostatic potential charges are indicated for the protonated water ...

In contrast to pT in long, open-ended water wires, where protonation of the central water molecule is energetically favorable,31 pT to the centre of the short water chain terminated by relatively strong proton acceptors is energetically uphill (Figures 1, ,2;2; Table 1). However, protonation of the central water molecule becomes 3 kcal mol−1 more favourable when the chain length is increased from two to four water molecules, reflecting the favourable charge-dipole interaction of the proton with two oriented water chains on either side (see below).

Moreover, the effective charge of the protonated water cluster increases from +0.38e to +0.70e upon increasing the chain length from N=1 to 4 water molecules (Figure 2), resulting in a polarizing effect that strengthens the charge-dipole interactions. The remaining positive charge is transferred to the D/A groups, which thus exert an electrostatic repulsion on the protonated water cluster. Increased charge localization on the protonated cluster, a decreased charge transfer to the D/A groups, and a more favourable charge-dipole interaction will thus lower the pT barrier with increasing chain length, consistent with the energetic trends in Figure 2.

To study quasi one-dimensional pT reactions, constraints were imposed on the distances between the terminal donor and acceptor groups, mimicking the effect of a relatively stiff protein backbone. Removal of these distance constraints lowered the pT barriers, but changed the mechanism. For a system with one water molecule, upon re-optimization without constraints the imidazole planes tilted by 50 degrees, resulting in a short hydrogen bond of 1.47 Å distance and a low pT barrier of only ~4 kcal mol−1. However, since such a model significantly deviates from the quasi-one dimensional geometry of interest, we will refer only to the constrained systems in the following.

Embedding the water chains in a dielectric environment increases the pT barrier heights. For a dielectric constant of ε=4, a value thought to be typical of the protein interior, the barriers increased by ~2 kcal mol−1 in calculations using the COSMO50 model (Table 2). A dielectric constant of 80 increased the barriers by an additional ~1 kcal mol−1 above the ε=4 values. As a consequence, low-dielectric environments favour fast pT reactions along short proton wires. We found that the dielectric lowers the energy of the reactant state more than that of the transition state. A possible explanation for this difference in solvation energies is that in the reactant state all water dipoles are aligned, in contrast to the symmetric transition state, where they oppose each other. Note that in our calculations, we only considered a static dielectric, with dielectric relaxation dynamics in proteins covering a broad range of time scales.5860

Table 2
Effect of the dielectric constant on the pT barriers ΔE* in units of kcal mol−1, calculated with the COSMO50 model.

The results for the dielectric-constant and chain-length dependence of pT have important implications for biological systems. Assuming that the homogenous dielectric model used here is at least qualitatively relevant, the increased barrier in a polar environment suggests that the low-dielectric environments observed in many enzymes not only help organize the water chains, but also lower the barrier for pT by effectively increasing the relative affinity of the water chain for the proton. Examples for such systems include the non-polar cavity above Glu-242 in cytochrome c oxidase35,6165 and the pT cavity between Asp-251 and the active site heme in cytochrome P450.6668 As mentioned above, short water chains seem to have different properties compared to the long open-ended water wires, which may capture elements of pT reactions in longer biological channels and pores such as aquaporin,6971 gramicidin,13,14 and the D-channel of cytochrome c oxidase.7277 According to our results, protonation of long water wires is energetically more favourable than protonation of short ones. Favourable protonation of a water cluster in the D-channel of cytochrome c oxidase was recently suggested.74 However, proton uptake through the long water wire in this channel is most probably linked to deprotonation of residues at the top of the channel,61,62,78 and might thus affect the stability of a protonated water cluster.

3.2 Thermodynamic properties and kinetic isotope effects

Vibrational contributions to the free energy barriers of pT can be estimated from nuclear zero-point energies (ZPE). This simplified analysis neither considers the proton as a quantum particle, nor does it treat the dynamics quantum mechanically.36 More refined treatments were used to study, e.g., pT in soybean lipo-oxygenase, a system with a large KIE ~ 80.79

Table 3 lists the ZPEs and thermodynamic contributions to the free energy barrier computed from a harmonic approximation to the energy surface via the molecular Hessian, i.e., the second derivative of the potential energy with respect to nuclear coordinate displacements. The resulting KIEs range from 2 to 10, which are values typical of systems without significant proton tunnelling, with a classical limit considered to be 7.9.51,80 The KIEs are also similar to the values of 1–7 obtained for different steps of allylic hydroxylation of propane by a model system of cytochrome P450 compound I at the B3LYP-level of theory.52 Both KIEs and ZPEs decrease with increasing chain length. The reaction free energies (ΔG*) are 2.9–4.3 kcal mol−1 lower relative to the electronic energy barriers, corresponding to rate accelerations by 2–3 orders of magnitude on the basis of transition state theory. The computed reaction free energy barriers of 9–13 kcal mol−1 correspond to pT rates of ~1/μs to 1/ms for a kinetic prefactor of 1/ps, a range consistent with experimentally observed pT rates in carbonic anhydrase and some of its mutants,81 and in cytochrome c oxidase.35,61

Table 3
Activation energies (ΔE*), zero-point energies (ΔZPE*), translational, rotational, and vibrational contributions to the internal energy (thermal) and entropy (TΔS*) at T=298 K, and resulting reaction free energies (ΔG*) ...

The difference between the reaction free energies and the electronic energies are mainly due to the entropic term, TΔS*, shifting the barriers by 1–3 kcal mol−1, and the zero-point effects, decreasing the barrier by 2.4 – 4.5 kcal mol−1. The entropic contributions TΔS* lower the barrier for N>2, consistent with a softening of the transition state structures with increasing chain length.

3.3. Dynamics of proton transfer

The energetic analysis identified key structures along the putative pT pathway. To explore their dynamic relevance, we constructed transition pathways for the pT reaction between the donor and acceptor groups that pass exactly through the central transition state structure. Pairs of Newtonian trajectories on the Born-Oppenheimer surface were initiated from the symmetric transition state structure for N=4 water molecules with Maxwell-Boltzmann starting velocities of opposite sign. The two trajectory segments can thus be thought of as evolving forward and backward in time, respectively.40 After inverting the velocities of the backward segment, the two trajectories are stitched together to form a transition path.41 Nine of the ten trajectory pairs went to opposite sides, statistically consistent with a committor value of ½ for the symmetric transition state and a transmission coefficient of nearly 1. On the basis of these tests we conclude that the initial structure is a dynamically relevant transition state. However, to construct a properly weighted transition path ensemble, we would have to include additional pathways that do not pass exactly through the saddle-point structure, and to re-weight the pathways according to the time spent at the surface from which they were started.40

Figure 3 shows three transition trajectories passing through the transition state. The plot of dtot versus time indicates some barrier recrossing events for this coordinate. These recrossings, and the small loops in the E-versus-dtot plot in Figure 1, highlight the stochastic nature of the pT process at finite temperature. The transition is complete within 200 to 300 fs. The gradual progression of dtot between 1.2 and 1.6 Å is consistent with the relatively flat barrier top implied by the heights of the three transition states in Figure 1.

Figure 3
pT transitions paths along a chain of N=4 water molecules. The figure shows three transition path trajectories (red, blue, and purple lines) passing through the Zundel-like transition state, as a function of the average progression coordinate, dtot. The ...

As a consequence, almost the entire transition path duration is spent with the proton progressing along the N=4 water chain, whereas the deprotonation and protonation events of the donor and acceptor groups, respectively, are comparably fast.

Figure 1 compares the potential energies along one of the transition paths to the energies of the optimized transition states. The energies of the dissociation trajectory are relative to the state where the proton reaches the terminal imidazole residue to correct for the thermally induced increase in potential energy (resulting in some ambiguity). With this correction, the dissociation trajectory remains energetically close to the idealized Zundel structures and transition state intermediates. As the transition states in the system with N=4 water molecules consist of a Zundel- and two hydronium-like structures (Figure 2), the potential energy curves suggest that pT proceeds by alternation between Zundel-like and hydronium-like species. Snapshots from the MD trajectories and analysis of the dynamically significant species further support this interpretation (Figure 4). The Zundel intermediate at the centre of the water chain dissociates into a hydronium-like ion localized at the second water molecule, which further protonates the imidazole via a transient Zundel ion.

Figure 4
Dynamics of pT obtained from transition paths. A) Average distance rave of water oxygen atoms to the three nearest hydrogen atoms for one transition path as a function of time. Zundel-like species have two neighbouring water molecules near rave~1.08 Å; ...

The charge dissociation/recombination reaction follows a semi-concerted reaction mechanism, where formation of a high-energy hydronium ion at the first water molecule is avoided by concertedly extending its distal hydrogen bond towards the second water molecule (Figure 4). To quantify the degree of concertedness of the pT reaction, we examined the correlation between the lengths d0 and d1 of the N-H and O-H bonds, respectively, of the donor and the first water molecule in the hydrogen-bonded chain. As shown in Figure 5, we performed linear fits of d1(t) as a function of d0(t) at equal times t for the 20 dissociation trajectories in the window of imidazole-water d(O-H) distances of 1.3 Å < d0 < 1.7 Å, the transition region between a covalent bond and a hydrogen bond. The line resulting from this fit has a slope of 0.81 (r2 = 0.74). We also found that two lines fitted separately for the “covalent” region of 1.3 Å < d0 < 1.5 Å, and the “hydrogen bonding” region of 1.5 Å < d0 < 1.7 Å coincide, consistent with a semi-concerted nature of the pT process. For a “pure” stepwise process, the two lines would cross at an angle to each other. We observed that along the transition paths even the bond lengths d0 and d2 are weakly correlated, with a slope of 0.19 (r2 = 0.53; not shown), further supporting the notion of a semi-concerted pT process.

Figure 5
Concertedness of the pT reaction. A) Free energy plot –kBT ln p(d0,d1) of the joint probability density of the N-H distance d0 of the proton donor, and the O-H distance d1 of the first water molecule in the chain, averaged over 20 MD trajectories ...

The semi-concerted reaction mechanism leads to a larger charge delocalization in the water wire in comparison to a fully stepwise-mechanism. A semi-concerted pT process might explain the dynamical behaviour observed, e.g., in QM/MM simulations of pT between Asp-95 and the photoisomerized Schiff base in bacteriorhodopsin,82 where pre-organization of a water chain connecting the donor and acceptor pair occurs prior to pT. A highly ordered water chain with tight hydrogen bonds is expected to be a prerequisite for a semi-concerted pT process, as suggested also for pT between an acid and a base in liquid water.83

It was observed in recent Car-Parrinello molecular dynamics simulations on water-mediated proton dissociation from 2,4,6-tricyanophenol that a sequential proton hopping process changed to a concerted dissociation process through reorganization of the water solvation shell.84 Concerted/semi-concerted pT processes have also been suggested by time-resolved mid-infrared measurements where the distance dependence of the pT rate was studied in a photoexcited acid-base pair.85 Semi-empirical tight-binding DFT calculations revealed that water-mediated pT between two protonation sites in fluorescein, shows a high degree of concertedness, most likely due a highly strained transition state, which also contributed to a high activation energy of ~ 15 kcal/mol.86 These findings are consistent with the reaction barriers and transition state structures found for the quasi-one dimensional proton wires studied here. Since we observed correlation between the motions of neighbouring hydrogen bonds, the pT process shares aspects with the “delocalized soliton transition” of Stuchebrukhov.87

A stepwise mechanism can become competitive with a semi-concerted mechanism by stabilization of the protonated water molecule next to the proton donor. For example, we found in calculations for a one-water system that the proton dissociation barrier is lowered by ~7 kcal mol−1 when the water molecule/hydronium ion is hydrogen bonded to a carbonyl group of formaldehyde, added to mimic the effect of a backbone carbonyl. The backbones of Val-247 and Val-243 in cytochrome P450 and cytochrome c oxidase,88 respectively, might serve such functions.

By using a multistate-empirical valence bond (MS-EVB) approach, Cao et al.33 recently observed that pT reactions in long open-ended water wires inside nanotubes are dominated by transitions between Zundel-like species that are possibly connected by a H7O3+ intermediate, in contrast to pT in solution where transitions between Eigen cations via intermediate Zundel cations dominate.5,8,9 The present simulations suggest that pT in short water wires proceeds through alternating hydronium-like and Zundel-like intermediates.

The presence of the donor and acceptor (D/A) groups affects the energetic ordering of the intermediate species. To quantify this effect, we performed single-point energy calculations of previously optimized transition states from which the D/A groups were removed. In the resulting open water chains, the hydronium transition state was destabilized (Figure 2D) relative to the Zundel intermediate (Figure 2E) by ~2 kcal mol−1, in accordance with simulations of open-ended pT wires, where protonation of the central water molecule is thermodynamically favourable.3032,34 The charge of the protonated species in the system without the D/A groups increases from 0.7e to 0.79e, indicating a more localized charge distribution relative to the systems with end groups.

3.4 Driving force effects

Chemical modifications of the proton-donating group alter the barrier height. As shown in Figure 6A, we found that the activation energy ΔE* increases linearly with the energetic driving force ΔE upon chemical modification of the donor group, where ΔE is defined as the energy difference of the product and reactant states. The corresponding Brønsted slope89 is 1. This linear dependence extends over a range of ~10 kcal/mol, from halogenations that increase the driving force (negative ΔE) to methylations that decrease the driving force (positive ΔE). In contrast, the barrier height was found to be insensitive to chemical modifications of the acceptor group, with a Brønsted slope of 0 (Figure 6B).

Figure 6
Activation energy ΔE* as a function of the energetic driving force ΔE for a system with N=4 water molecules. A) ΔE* as a function of ΔE for proton donors chemically modified by halogenation or methylation. B) ΔE ...

The energy diagrams shown in Figure 7 provide a qualitative explanation for the observed behaviour. In the underlying model, we assume that chemical modifications of the donor (acceptor) do not significantly alter the proton affinity of the acceptor (donor) and of the intervening water chain. With this simplifying assumption, the barrier height is determined completely by the difference in energy between the protonated water chain and the donor, consistent with the Brønsted slopes of 1 and 0 in Figure 6A and B, respectively.

Figure 7
Schematic energy level diagrams for driving force effects. A) Destabilization of the proton donor (D) decreases the reaction barrier, consistent with the behaviour observed in Figure 6A. B) Destabilization of the acceptor (A) leaves the reaction barrier ...

The Brønsted slopes of 1 and 0 are consistent with the relatively small changes in the structure of the protonated water chain in Figure 6C and the small differences in the energies of the three transition states in Figure 1. As result of the halogenations, the character of the transition state changes from a Zundel ion to a hydronium-like species, which can be explained by the increased proton affinity of the water molecule forming the hydronium ion. This pKa-perturbation is largely due to electrostatic effects. Halogenation results in a charge transfer of 0.05e from the acceptor to the donor that electrostatically moves the transition state closer to the reactant state.

4. Conclusions

We showed that pT reactions in short water wires connecting proton donating and accepting groups differ qualitatively from pT reactions in long and open-ended proton wires. For single-file chains of N=1 to 4 water molecules between methyl-imidazole and methyl-imidazolium, high pT reaction barriers of ~13 kcal mol−1 were obtained by using hybrid density functional theory and a triple-ζ basis set. Entropic effects and zero-point energies (ZPE) decrease the reaction barriers by 3–4 kcal mol−1. An increased water chain length and a low-dielectric environment also reduced the barrier. KIEs were found to be in the range from 2–10, and to decrease with increased chain length. The quantum chemical cluster models used in our calculations are expected to capture the essence of biological pT processes, where the proton donor and acceptor have moderate pKa values, and the donor-acceptor groups can be separated by ~10 Å, and connected by 1–4 water molecules.

Electrostatic interactions are important factors in determining the barrier height. In our calculations, the effective charge of the protonated water clusters in the short chains varied between +0.4e and +0.7e for short and long water chains. Enhanced charge-dipole interactions in the longer chain help solvate the increased excess charge, together with weakened charge-charge repulsion from the end-groups, which accept the remaining positive charge from the protonated water cluster.

In our molecular dynamics simulations, the pT reaction in the N=4 water chain followed a semi-concerted reaction mechanism. In dynamic transition paths on the Born-Oppenheimer surface, the proton motions in the water chain were found to be correlated, with tightly coupled transitions between hydronium and Zundel-like intermediates. This concerted motion reflects in part the high energetic cost of forming ideal hydronium ions, which are avoided in a concerted transfer.

Changes in the energetic driving force as a result of chemical modifications affected the barrier height only from the donor side, but not from the acceptor side. For donor modifications and a chain of N=4 water molecules, we observed a linear decrease in the reaction barrier with increasing driving force corresponding to a Brønsted slope of one. In contrast, acceptor modifications did not change the reaction barrier, corresponding to a Brønsted slope of zero. This unusual behaviour implies that the proton affinity of the N=4 water chain is not perturbed significantly by chemical modifications of the donor or acceptor groups. In biological systems, the rate of enzymatic pT would thus be most sensitive to donor perturbations.


V.R.I.K. acknowledges Dr. Jürgen Köfinger and Prof. Mårten Wikström for discussions on pT. This research was supported by the Intramural Research Program of the National Institute of Diabetes and Digestive and Kidney Diseases, National Institutes of Health. V.R.I.K. also acknowledges the European Molecular Biology Organization (EMBO) for a Long-Term Fellowship. CSC - the Finnish IT Center for Science and Biowulf cluster at NIH are acknowledged for computer time.


1. Agmon N. Chem Phys Lett. 1995;244:456.
2. Nagle JF, Morowitz HJ. Proc Natl Acad Sci USA. 1978;75:298. [PubMed]
3. Onsager L. In: The Neurosciences. Quarton G, Melnechuk T, Schmitt FO, editors. Rockefeller University Press; New York: 1967. p. 75.
4. Wraight CA. Biochim Biophys Acta. 2006;1757:886. [PubMed]
5. Marx D, Tuckerman ME, Hutter J, Parrinello M. Nature. 1999;397:601.
6. de Grotthuss CJT. Ann Chim (Paris) 1806;58:54.
7. de Grotthuss CJT. Biochim Biophys Acta. 2006;1757:871. [PubMed]
8. Lapid H, Agmon N, Petersen MK, Voth GA. J Chem Phys. 2005;122:014506. [PubMed]
9. Markovitch O, Chen H, Izvekov S, Paesani F, Voth G, Agmon N. J Phys Chem B. 2008;112:9456. [PubMed]
10. Rosta E, Woodcock HL, Brooks BR, Hummer G. J Comput Chem. 2009;30:1634. [PMC free article] [PubMed]
11. Cui Q, Karplus M. J Am Chem Soc. 2002;124:3093. [PubMed]
12. Dieterich JM, Werner HJ, Mata RA, Metz S, Thiel W. J Chem Phys. 2010;132:035101. [PubMed]
13. Pomès R, Roux B. Biophys J. 1996;71:19. [PubMed]
14. Pomès R, Roux B. Biophys J. 1998;75:33. [PubMed]
15. Chakrabarti N, Roux B, Pomès R. J Mol Biol. 2004;343:493. [PubMed]
16. Kato M, Pisliakov AV, Warshel A. Proteins Struct Funct Bioinformatics. 2006;64:829. [PubMed]
17. Eigen M. Angewandte Chemie International Edition. 1964;3:1.
18. Zundel G, Metzger H, Physik Z. Chem (NF) 1968;58:225.
19. Marcus RA, Sutin N. Biochim Biophys Acta. 1985;811:265.
20. Beratan DN, Onuchic JN, Winkler JR, Gray HB. Science. 1992;258:1740. [PubMed]
21. Moser CC, Keske JM, Warncke K, Farid RS, Dutton PL. Nature. 1992;355:796. [PubMed]
22. Page CC, Moser CC, Chen X, Dutton PL. Nature. 1999;402:47. [PubMed]
23. Borgis D, Hynes JT. Chem Phys. 1993;170:315.
24. Borgis D, Hynes JT. J Chem Phys. 1991;94:3619.
25. Villa J, Warshel A. J Phys Chem B. 2001;105:7887.
26. Hwang J-K, Chu ZT, Yadav A, Warshel A. J Phys Chem. 1991;95:8445.
27. Hammes-Schiffer S, Stuchebrukhov AA. Chem Rev. 2010;110:6939. [PMC free article] [PubMed]
28. Hummer G, Rasaiah JC, Noworyta JP. Nature. 2001;414:188. [PubMed]
29. Vaitheeswaran S, Rasaiah JC, Hummer G. J Chem Phys. 2004;121:7955. [PubMed]
30. Hassan SA, Hummer G, Lee Y-S. J Chem Phys. 2006;124:204510. [PMC free article] [PubMed]
31. Dellago C, Naor MM, Hummer G. Phys Rev Lett. 2003;90:105902. [PubMed]
32. Dellago C, Hummer G. Phys Rev Lett. 2006;97:245901. [PubMed]
33. Cao Z, Peng Y, Yan T, Li S, Li A, Voth GA. J Am Chem Soc. 2010;132:11395. [PubMed]
34. Köfinger J, Hummer G, Dellago C. Proc Natl Acad Sci USA. 2008;105:13218. [PubMed]
35. Kaila VRI, Verkhovsky MI, Wikström M. Chem Rev. 2010;110:7062. [PubMed]
36. Kiefer PM, Hynes JT. J Phys Org Chem. 2010;23:632.
37. Becke AD. J Chem Phys. 1993;98:5648.
38. Lee C, Yang W, Parr RG. Phys Rev B. 1988;37:785. [PubMed]
39. Schäfer A, Horn H, Ahlrichs R. J Chem Phys. 1992;97:2571.
40. Dellago C, Bolhuis PG, Csajka FS, Chandler D. J Chem Phys. 1998;108:1964.
41. Hummer G. J Chem Phys. 2004;120:516. [PubMed]
42. Weigend F, Ahlrichs R. Phys Chem Chem Phys. 2005;7:3297. [PubMed]
43. Grimme S. J Chem Phys. 2003;118:9095.
44. Christiansen O, Koch H, Jørgensen P. Chem Phys Lett. 1995;243:409.
45. Perdew JP. Phys Rev B. 1986;33:8822. [PubMed]
46. Becke AD. Phys Rev A. 1988;38:3098. [PubMed]
47. Tao J, Perdew JP, Staroverov VN, Scuseria GE. Phys Rev Lett. 2003;91:146401. [PubMed]
48. Grimme S. J Comput Chem. 2004;25:1463. [PubMed]
49. Grimme S. J Comput Chem. 2006;27:1787. [PubMed]
50. Klamt A, Schüürmann G. J Chem Soc Perkin Trans. 1993;2:799.
51. Melander L, Saunders WH., Jr . In: Reaction Rates of Isotopic Molecules. Robert E, editor. Krieger Publishing Company; Malabar, Florida: 1987.
52. de Visser SP, Ogliaro F, Sharma PK, Shaik S. J Am Chem Soc. 2002;124:11809. [PubMed]
53. Ogliaro F, Filatov M, Shaik S. Eur J Inorg Chem. 2000:2455.
54. Bell RP. Chem Soc Rev. 1974;3:513.
55. Ahlrichs R, Bär M, Häser M, Horn H, Kölmel C. Chem Phys Lett. 1989;162:165. Current version: see
56. Neese F, Becker U, Ganyushin D, Kossmann S, Hansen A, Liakos D, Petrenko T, Riplinger C, Wennmohs F. ORCA, version 270. University of Bonn; Bonn, Germany: 2009.
57. Humphrey W, Dalke A, Schulten K. J Mol Graphics. 1996;14:33. [PubMed]
58. Jordanides X, Lang MS, Fleming G. J Phys Chem B. 1999;103:7995.
59. Abbyad P, Childs W, Shi X, Boxer SG. Proc Natl Acad Sci USA. 2007;104:20189. [PubMed]
60. Leontyev IV, Stuchebrukhov AA. J Chem Phys. 2009;130:085103. [PubMed]
61. Belevich I, Bloch DA, Belevich N, Wikström M, Verkhovsky MI. Proc Natl Acad Sci USA. 2007;104:2685. [PubMed]
62. Kaila VRI, Verkhovsky MI, Hummer G, Wikström M. Proc Natl Acad Sci USA. 2008;105:6255. [PubMed]
63. Pisliakov AV, Sharma PK, Chu ZT, Haranczyk M, Warshel A. Proc Natl Acad Sci USA. 2008;105:7726. [PubMed]
64. Wikström M, Verkhovsky MI, Hummer G. Biochim Biophys Acta. 2003;1604:61. [PubMed]
65. Xu J, Voth GA. Biochim Biophys Acta. 2008;1777:196. [PMC free article] [PubMed]
66. Schlichting I, Berendzen J, Chu K, Stock AM, Maves SA, Benson DE, Sweet BM, Ringe D, Petsko GA, Sligar SG. Science. 2000;287:1615. [PubMed]
67. Denisov IG, Makris TM, Sligar SG, Schlichting I. Chem Rev. 2005;105:2253. [PubMed]
68. Taraphder S, Hummer G. J Am Chem Soc. 2003;125:3931. [PubMed]
69. Murata K, Mitsuoka K, Hirai T, Walz T, Agre P, Heymann JB, Engel A, Fujiyoshi Y. Science. 2000;407:599. [PubMed]
70. de Groot B, Grubmuller H. Science. 2001;294:2353. [PubMed]
71. Tajkhorshid E, Nollert P, Jensen M, Miercke L, Stroud RM, Schulten K. Science. 2002;296:525. [PubMed]
72. Tsukihara T, Aoyama H, Yamashita E, Tomizaki T, Yamaguchi H, Shinnzawha-Itoh K, Nakashima R, Yaono R, Yoshikawa S. Science. 1995;269:1069. [PubMed]
73. Iwata S, Ostermeier C, Ludwig B, Michel H. Nature. 1995;376:660. [PubMed]
74. Xu J, Sharpe MA, Qin L, Ferguson-Miller S, Voth GA. J Am Chem Soc. 2007;129:2910. [PMC free article] [PubMed]
75. Riistama S, Hummer G, Puustinen A, Dyer RB, Woodruff WH, Wikström M. FEBS Lett. 1997;414:275. [PubMed]
76. Olkhova E, Hutter MC, Lill MA, Helms V, Michel H. Biophys J. 2004;86:1873. [PubMed]
77. Henry RM, Yu CH, Rodinger T, Pomès R. J Mol Biol. 2009;387:1165. [PubMed]
78. Gorbikova EA, Belevich NP, Wikström M, Verkhovsky MI. Biochemistry. 2007;46:13141. [PubMed]
79. Olsson MHM, Siegbahn PEM, Warshel A. J Am Chem Soc USA. 2004;126:2820. [PubMed]
80. Huynh MHV, Meyer TJ. Chem Rev. 2007;107:5004. [PMC free article] [PubMed]
81. Silverman DN. Biochim Biophys Acta. 2000;1458:88. [PubMed]
82. Lee YS, Krauss M. J Am Chem Soc. 2004;126:2225. [PubMed]
83. Siwick BJ, Bakker HJ. J Am Chem Soc. 2007;129:13412. [PubMed]
84. Maurer P, Thomas V, Iftimie R. J Chem Phys. 2011;134:094505. [PubMed]
85. Cox MJ, Timmer RLA, Bakker HJ, Park S, Agmon N. J Phys Chem A. 2009;113:6599. [PubMed]
86. Friedman R, Fischer S, Nachliel E, Scheiner S, Gutman M. J Phys Chem B. 2007;111:6059. [PubMed]
87. Stuchebrukhov AA. Phys Rev E. 2009;79:031927. [PubMed]
88. Amino acid numbering refers to the sequences of cytochrome P450 from Pseudomonas putida and cytochrome c oxidase from Bos taurus.
89. Brønsted JN, Pedersen K. Z Phys Chem. 1924;108:185.