Article pubs.acs.org/JPCA
Quantum Mechanical Enhancement of Rate Constants and Kinetic Isotope Effects for Water-Mediated Proton Transfer in a Model Biological System James W. Mazzuca* and Chase P. Schultz Chemistry Department, Alma College, Alma, Michigan 48801, United States S Supporting Information *
ABSTRACT: Biological systems have been shown to shuttle excess protons long distances by taking advantage of tightly organized hydrogen-bonded water bridges in hydrophobic protein cavities, and similar effects have been observed in carbon nanotubes. In this theoretical study we investigate how quantum effects of proton motion impact the rate constants for charge transfer in a model system consisting of a donor and acceptor molecule separated by a single-molecule water bridge. We calculate quantum and classical rate constants for the transfer of an excess proton over two possible paths, one with an H3O+ intermediate, and one with an OH− intermediate. Quantum effects are included through ring polymer molecular dynamics (RPMD) calculations. We observe a 4-fold enhancement of reaction rate constants due to proton tunneling at temperatures between 280 and 320 K, as shown by transmission coefficient calculations. Deuteration of the donor and acceptor proton are shown to decrease the reaction rate constant by a factor of 50, and this is another indicator that tunneling plays an important role in this proton transfer mechanism.
■
INTRODUCTION Efficient, fast charge transfer in chemical systems is essential for biological life and a highly desirable feature of engineered chemical products. There are many possible ways that a chemical system can move a charge through space, and two mechanisms commonly observed in nature are proton transfer and electron transfer. These charge transfer mechanisms frequently consist of a single step when occurring over a short distance (1−5 Å), but in the case of long-distance charge transfer (10−30 Å), the mechanism can be more complicated. In the case of proton transfer in biological systems, nature must find a way to mitigate the problem of moving an excess proton over such a long distance within a reasonable time scale. Nearly all biological processes take place in an aqueous environment. In the presence of a hydrophobic protein tunnel, a “bridge” or “wire” of highly organized water molecules can form. These water bridges can aid in long-distance charge transfer through a sort of proton-hopping mechanism where an excess proton is donated to one end of the water bridge and ultimately taken up by an acceptor molecule at the other end, a process commonly referred to as the Grotthuss mechanism.1,2 There is evidence to suggest that biological systems such as the gramicidin A transmembrane channel,3−5 cytochrome C oxidase,6,7 bacteriorhodopsin protein,8,9 and carbonic anhydrase10 form hydrophobic intraprotein tunnels which take advantage of Grotthuss-like behavior and increase the rate of proton transfer. Some proteins will include hydrophilic functional groups in these tunnels, such as carboxylic acids,11 in an effort to control the rate of charge transfer. © XXXX American Chemical Society
Some biological systems, such as the aquaporin class of transmembrane proteins, will form water bridges that appear very similar to the ones mentioned above, but are impermeable to protons. It is currently understood that the protein environment arranges the water molecules in a way that, to move an excess proton along the wire, would introduce an insurmountable free energy barrier.12−14 A Grotthuss-like charge transfer is impossible in these proteins, and as a result, water molecules can move freely through the channel without disrupting the transmembrane charge gradient. The Grotthuss mechanism of proton transfer has also been studied in the context of bulk water,15 small molecules in aqueous solution,16 carbon nanotubes,17,18 and mesoporous aluminosilicate model systems.19 While there are many conceivable instances of such a charge transfer mechanism, we will focus our attention on biological systems. Understanding how biological systems transfer charge in this way may provide insight into the chemical physics that govern these processes and the steps nature has taken to optimize them. It has been shown that the quantum motion of nuclei affects the structure of water bridges20 as well as structural and dynamical properties of large hydrogen-bonded networks of water molecules.21 Theoretical studies on the hydrated proton in bulk water have shown that the free energy barriers as well as H/D kinetic isotope effects (KIE) for proton transfer between water molecules are substantially impacted by the inclusion of Received: October 12, 2016 Revised: December 7, 2016 Published: January 4, 2017 A
DOI: 10.1021/acs.jpca.6b10337 J. Phys. Chem. A XXXX, XXX, XXX−XXX
Article
The Journal of Physical Chemistry A quantum nuclear effects.22−24 We expect that the highly organized and constrained water wires that form in biological systems will amplify quantum effects such as zero-point energy (ZPE) and proton tunneling. In this study, we calculate rate constants for proton transfer in a small model system which resembles a donor−water− acceptor intraprotein environment with a water bridge consisting of a single molecule. We perform calculations in which we can include quantum effects of proton motion such as ZPE and tunneling, and we examine how rate constants for proton transfer are affected. We then investigate how deuterating the system affects the reaction rate constants and determine whether this behavior is indicative of a classical or quantum proton transfer mechanism. The results of this study may be useful in providing a basic physical understanding of this phenomenon and a direction for future work in this area.
■
Figure 2. Potential energy surface generated by varying the bond lengths indicated in Figure 1 and calculating single-point energy values at the B3LYP/6-311+G** level of theory. The four marked locations on the potential energy surface (PES) correspond to the four structures sketched below.
MODEL SYSTEM We focus our study on a model system developed by Kaila and Hummer which consists of two anchored imidazole molecules, a central water molecule, and an excess proton.25 This model was chosen to resemble charge transfer between two histidine residues at different locations within a protein. A sketch of this system can be seen in Figure 1. The terminal methyl carbon
where xj is the position in the jth degree of freedom, and x = x1, x2, ..., x6. mH is the mass of a proton, and ωj is the frequency of mode j and is defined as ⎛ j − 3/2 ⎞ ⎟ ωj = −ωc ln⎜ ⎝ 5 ⎠
(2) −1
ωc was chosen to be 3000 cm to mimic the frequency of an N−H bond and provide adequate coupling to nonreactive degrees of freedom through the coupling constant cj
Figure 1. Sketch of the model system used in this study. The two methyl carbons are anchored 12.8 Å apart, and the excess proton originates on the left imidizole molecule. The coordinates x1 and x2 correspond to the bond lengths that were varied to generate the 2-D potential energy surface in Figure 2.
cj = ωj
2mHωc 5π
(3)
The constants x1,min = 2.8a0 and x2,min = 2.7a0 were chosen to be near the top of the barrier in the center of the PES in Figure 2. The coupling between x1/x2 and the other degrees of freedom therefore increases as one moves away from the diagonal which passes through the center of this barrier, and we would expect strong coupling in the reactant and product states. As seen in Figure 2, there are two primary paths which connect the reactant state to the product state. We investigate both the 1 → 2 → 4 path which involves the H 3 O + intermediate, and the 1 → 3 → 4 path which involves the OH− intermediate. In order to determine how quantum effects directly influence the rate of proton transfer, it is desirable to employ a method which can be easily compared to a classical calculation. We employ the ring polymer molecular dynamics (RPMD) method.29−32 Using this approach, we can approximate the quantum mechanical rate constant33 by solving classical molecular dynamics equations for groups of particles in an extended phase space.34,35 This method avoids the exponential scaling of a fully quantum mechanical description, and instead results in a linear scaling with respect to the spatial dimensionality of our system.36 This approach provides the added benefit that classical dynamics is a limiting case of RPMD calculations. Variations of this method have been used to compute the dynamics of high-dimensional systems, including gas-phase reactions,37−39 intramolecular proton transfer,40 and properties of condensed-phase systems such as liquid water.41 While our current model system consists of only
atoms were fixed at 12.8 Å apart to mimic attachment to a protein backbone and provide enough space for a single water molecule. A symmetric transition state was calculated at the B3LYP/6-311+G** level of theory in Q-Chem 4.226 in which all degrees of freedom were optimized. Starting from the transition state geometry, single-point energy calculations, at the same level of theory, were performed by independently varying the values of x1 and x2 as shown in Figure 1. A 50 × 50 grid of energy values was generated in which x1 and x2 values range from 0.95 a0 (0.50 Å) to 4.35 a0 (2.30 Å) at 0.0694 a0 (0.0367 Å) increments, where a0 is the Bohr radius. Cubic spline interpolation between the grid points provided a continuous surface for dynamics calculations.27 We will refer to the two-dimensional potential energy surface as Vgrid(x1, x2) as seen in Figure 2. In order to provide energy dispersion and translational freedom for the two protons, we include four additional degrees of freedom in the form of a modified system-bath model.28 Our model therefore consists of two three-dimensional protons for a total of six degrees of freedom. The potential energy is V (x) = Vgrid(x1 , x 2) + 2⎤ ⎡ ⎛ x1 − x1,min + x 2 − x 2,min ⎞ ⎥ 1 ⎢ 2⎜ ⎟ ∑ ⎢ mHωj ⎜xj − cj ⎟⎥ mHωj2 ⎝ ⎠⎦ j=3 ⎣ 2 6
(1) B
DOI: 10.1021/acs.jpca.6b10337 J. Phys. Chem. A XXXX, XXX, XXX−XXX
Article
The Journal of Physical Chemistry A ‡
e−ΔG /kBT is the probability of a reactant state configuration reaching the dividing surface where ΔG⧧ is the free energy barrier at a temperature T. Both κ(T, t) and ΔG⧧ are evaluated numerically through dynamics calculations, and the results are combined to calculate a thermal rate constant.
two three-dimensional protons, the option to study highdimensional systems in the future was an important factor in choosing an appropriate method. The Hamiltonian of our classical system is written as pj2
6
H(x , p) =
∑ j=1
2mj
■
+ V ( x)
CALCULATIONS The free energy was calculated using the umbrella integration technique of Kästner and Thiel.45 This method provides an alternative to explicitly calculating the partition function in eq 7, and it provides a straightforward picture of system evolution from the reactant state to the dividing surface. The change in free energy ΔG⧧ can be calculated by integrating ∂G/∂x along the reaction coordinate x from the reactant state xR to the dividing surface xD
(4)
where p = p1, p2, ..., p6 and x = x1, x2, ..., x6; there are six degrees of freedom, or three for each proton. The RPMD method treats each physical particle as n replicates, or “beads”, which are connected to two neighboring beads by harmonic springs. The Hamiltonian of our two protons represented as a ring polymer is Hn(x1, x2 , ..., x n, p1, p2 , ..., pn) ⎡ n p2 mj j,i = ∑ ⎢∑ + ⎢ 2βn2ℏ2 j = 1 ⎣ i = 1 2mj 6
ΔG⧧ =
⎤ ∑ (xj ,i − xj ,i− 1)2 ⎥⎥ i=1 ⎦ n
∑ V (x i )
(5)
i=1
f (xj) ≡
where xi = x1,i, x2,i, ..., x6,i and pi = p1,i, p2,i, ..., p6,i for an n-bead ring polymer system where i is bead index. The boundary condition xj,0 = xj,n ensures that the ring polymer is continuous and each bead is connected to only two others. Quantum effects of chemical processes can be approximated by substituting ring polymers for classical particles at a temperature βn = β/n for an n-bead system at β = 1/kBT where kB is the Boltzmann constant and T is temperature. Increasing the number of beads is analogous to increasing the resolution of the calculation, and as such calculations should converge in the limit of large n. The general formula for a reaction rate constant is 1 k(T ) = lim c fs(T , t ) Q r(T ) t →∞
∫ dp0 ∫ dx 0e
g i (x ) ≡
t →∞
∑ f (xj)Δx (9)
j=1
∂G ∂x
nW
= xj
⎛ ∂Gi ⎞ ⎟ ∂x ⎠
∑ pi (xj)⎜⎝ i=1
nW
= xj
∑ pi (xj)gi(xj) i=1
b ∂Gi 1 x − xi̅ = − k W(x − xi) ∂x β (σib)2
(11)
and pi is a probability density which depends on all windows pi (x) =
Pib(x) nW b ∑i = 1 Pi (x)
(12)
where
(6)
Pib(x)
=
1 σib
⎡ b 2⎤ 1 ⎛ x − xi̅ ⎞ ⎥ ⎢ ⎟ exp − ⎜⎜ ⎢ 2 ⎝ σ b ⎟⎠ ⎥ 2π i ⎣ ⎦
(13)
The variance is written as = ⟨(x − and is the average position along the reaction coordinate during the dynamics. The superscript b indicates dynamics occur under the influence of a biasing potential
−βH(x 0, p0)
(σbi )2
δ(x 0)v(x 0, p0)h(x t )
where ℏ = h/2π and h is Planck’s constant. Subscripts 0 and t indicate time; δ(x0) is a delta function centered at x0. v(x0, p0) is velocity, and h(xt) is a Heaviside step function. This equation can be written in the same way for a ring polymer by replacing H → Hn and β → βn and integrating over all ring polymer beads.34 Solving eq 6 as it is written is not straightforward, so we employ the Bennett−Chandler approach to transition state theory which allows us to rewrite eq 6 in a way which can be solved numerically, and without any loss of generality.42−44 Within this formulation, the rate constant for a process in which reactants and products can be separated by a dividing surface can be written as ⧧
R
(10)
(7)
k(T ) = lim κ(T , t )e−ΔG
Nbins
f (x)dx ≈
Here we can approximate the derivative of the free energy in window i as
where Qr(T) is the partition function of the reactant state per unit volume, and cfs(T, t) is the flux-side correlation function33 1 c fs(T , t ) = 2π ℏ
xD
where xj is the location of a sampling bin and ∂G/∂x is approximated as
n
+
∫x
V ib(x) =
1 k W(x − xi)2 2
xbi̅ )2⟩,
xbi̅
(14)
0.5Eh/a20.
and we have chosen kW = The reaction coordinate x was defined as x1 for the 1→ 3 → 4 path, and x2 for the 1 → 2 → 4 path, and this coordinate was constrained by the potential in eq 14. Umbrella integration dynamics were performed for nW = 1024 equally spaced windows along the reaction coordinate in the presence of a Berendsen thermostant,46 and the average position xbi̅ and variance (σbi )2 were calculated for each window. The free energy was calculated by interpolating between these windows at Nbins = 30 000 equally spaced bins and then numerically integrating the resulting curve, as shown in eq 9. The difference between the free energy minimum at the reactant state and the free energy maximum is the ΔG⧧ which is used in eq 8.
/ kBT
(8)
where κ(T, t) is the transmission coefficient for a process which originates at the dividing surface xD, and the exponential term C
DOI: 10.1021/acs.jpca.6b10337 J. Phys. Chem. A XXXX, XXX, XXX−XXX
The Journal of Physical Chemistry A
■
All dynamics calculations, both in the calculation of ΔG⧧ and κ(T, t), employed a leapfrog integration algorithm with a time step size of 2.0 atomic units (0.04838 fs). In each of the 1024 umbrella integration windows, dynamics were run for 8 × 104 time steps, for a total time of 1.6 × 105 atomic units (7.740 ps), and the system was allowed to equilibrate for half of the total time before averages were tabulated for eq 9. To control for the initial conditions of each umbrella integration window, eight replicate calculations were performed in each window with different, randomly sampled initial conditions for all ring polymer beads. The location along the reaction coordinate that corresponds to ΔG⧧ is taken to be the dividing surface, xD. Once the dividing surface location is determined, the transmission coefficient κ(T, t) is calculated in two steps. First, we perform constrained dynamics at the dividing surface by pinning the centroid of the ring polymer to xD and allowing the system to evolve under the influence of a thermostat for 4000 time steps (193 fs). After this thermalization the constraints are relaxed, and the thermostat is turned off; each particle is randomly assigned a momentum from the Boltzmann distribution corresponding to T for classical particles or nT for ring polymer beads. Recrossing trajectories are evolved forward in time for 16 000 time steps (774 fs) which allows the system enough time to settle on a reactant or product state. The average speed of the ring polymer at time 0 normal to the dividing surface, v(x0, p0), is multiplied by a step function as the system evolves ⎧1 if x1 + x 2 ≥ x D h(x t ) = ⎨ ⎩ 0 if x1 + x 2 < x D
Figure 3. Here we see the change in free energy as the protonated and deuterated reactant proceeds to the H3O+ transition state at 300 K. The change in free energy, ΔG⧧, is the difference between the minimum and maximum along these curves. Not only does zero-point energy lower the effective barrier, but the spatial location of the minima and maxima are also affected by the RPMD inclusion of quantum effects. As expected, we see identical free energy profiles for the both isotopes in the classical case.
The free energy barrier is approximately 7 mEh lower in the quantum proton calculation when compared to its classical counterpart, and this has a dramatic impact on the rate constant as seen in eq 8 due to its place in the exponential term. This free energy discrepancy is largely a result of ZPE effects included in the RPMD calculations. Free energy curves for the other cases were similar, including along the 1 → 3 → 4 path. Deuteration has a substantial impact on the free energy curve when quantum effects are included, raising ΔG ⧧ by approximately 3 mEh, and the curve more closely resembles what is calculated for the classical case. Trends in the location of ΔG⧧ as a function of temperature for the quantum proton can be seen in Figure 4. As temperature is increased, the free energy maximum increases, but on a much smaller scale than what is seen upon the exclusion of quantum effects or isotopic substitution. The temperature in the exponential term of eq 8 substantially
(15)
and the x1 + x2 terms result from drawing the dividing surface at a 45° angle along the diagonal of the PES in Figure 2. Averaged over all trajectories sampled from the Boltzmann distribution, the transmission coefficient is written as κ(T , t ) = ⟨v0(x 0, p0)h(x t )⟩
RESULTS AND DISCUSSION
In every case, the calculated change in free energy led to the largest discrepancy between the quantum and classical results. The change in free energy along the reaction coordinate for the 1 → 2 → 4 path, in which the system passes through the H3O+ transition state, for proton transfer at 300 K can be seen in Figure 3.
⎪
⎪
Article
(16)
which is a measure of the speed of reactive trajectories that originate at the dividing surface xD. The long time limit of this function is used in the rate constant calculations as shown in eq 8. For the calculations reported here, 40 000 recrossing trajectories were run, and 32 beads were used in RPMD calculations to generate converged results. To sample a range of temperatures, both reaction paths, and the effect of isotopic substitution, 40 separate calculations were performed. At each temperature, T = {280, 290, 300, 310, 320} K, the transmission coefficients κ(T, t), change in free energy ΔG⧧, and rate constants k(T) were evaluated for (a) quantum proton, (b) classical proton, (c) quantum deuterium, and (d) classical deuterium transfer. All calculations were performed on the Blue Waters supercomputer at the University of Illinois at Urbana−Champaign. For each calculation, 1024 compute cores were utilized for approximately 30 min of wall-clock time. The code exhibited a linear scaling with respect to number of cores, and this behavior is expected because of the parallels to classical molecular dynamics calculations. Each umbrella integration window and recrossing trajectory were evaluated independently, and the results were consolidated to calculate the rate constants.
Figure 4. Free energy maximum for the quantum proton which proceeds through the H3O+ transition state for all temperatures. The free energy barrier slightly increases as the temperature is raised, but the effect on the rate constant is diminished by the kBT denominator in eq 8. D
DOI: 10.1021/acs.jpca.6b10337 J. Phys. Chem. A XXXX, XXX, XXX−XXX
Article
The Journal of Physical Chemistry A
the ΔG⧧ calculation in the sense that much more conformational space is explored by the recrossing trajectories, and quantum tunneling may contribute in unexpected ways. Rate constants were calculated for each of the 40 cases outlined above, and ratios of quantum to classical rate constants kQ/kC were calculated for each. Very similar results were seen for both reaction paths, so we show results for just the 1 → 2 → 4 (H3O+ transition state) path in Figure 7. Rate constants are
outweighs these small changes in free energy, and the overall rate constant increases as temperature is increased. The transmission coefficient determines which trajectories, assuming they make it to the dividing surface, are actually reactive. Normalized transmission coefficients κ(T, t)/⟨v(x0, p0)⟩ for proton transfer at 300 K in the 1 → 2 → 4 path are shown in Figure 5. We see that, in the long
Figure 5. Normalized transmission coefficient for the protonated system originating in the H3O+ transition state at 300 K. We see much less recrossing when quantum effects are included, which may indicate one way in which tunneling increases the rate constant. From this plot we see a 4-fold increase in the rate constant for quantum hydrogen, and about half that for deuterium when compared to classical calculations.
Figure 7. We observe an over 1000-fold increase in rate constants due to quantum effects for the proton transfer reaction which progresses through the H3O+ transition state. Nearly identical results are observed for the OH− transition state.
enhanced by a factor of 100 for deuterium, and over 1000 for hydrogen, and this is due to a combination of both ZPE and quantum tunneling effects. A summary of all calculated ΔG⧧, limt→∞ κ(T, t), and k(T) values can be seen in Table 1. Since the ratios shown in Figure 7 are not something that could be observed experimentally, we found it useful to calculate the kinetic isotope effects, KIE = kH/kD for each case. The results of these calculations can be seen in Figure 8. We observe a KIE ranging from approximately 20−50, and an abnormally high KIE such as this can be indicative of a proton transfer mechanism where tunneling plays a large role. An increase in temperature lowers the kinetic isotope effect, but all temperatures in the range studied show an effect far larger than the KIE ≈ 2 = mD /mH calculated for the classical case. An accurate description of the proton transfer reaction in this model system therefore requires a quantum description of nuclei. Large isotope effects such as this have been observed for proton transfer in enzymes such as the hydrogen abstraction seen in soybean lipoxygenase-1 (KIE = 81),47,48 horse liver alcohol dehydrogenase (KIE = 7.5),49 and methylamine dehydrogenase (KIE = 16.8).50 These proton transfer reactions are single-step abstractions directly from hydrocarbon donor molecules to the active site acceptor moiety, and unlike our model system, there are no intervening water molecules. In these three enzymatic systems, a proton tunneling mechanism is suspected to be responsible for the abnormally high observed KIE values. We propose a similar explanation for what is measured in our model system: a two-proton transfer that is strongly influenced by quantum tunneling and ZPE effects. One possible source of error in these calculations is an overestimation of the potential energy barrier. An artificially high potential energy barrier would lead to an overestimation of quantum tunneling effects and, as a result, KIE values that are too high. This can result from the level of theory used in electronic structure calculations, coupling to the four nonreactive vibrational modes, and the manner in which x1 and x2 were scanned through a series of single-point calculations.
time limit, the transmission coefficient is substantially larger in the quantum case. Assuming the reactant proceeds to the H3O+ transition state, the quantum proton is 4 times more reactive than the classical case, and we observe about half of this effect when the system is deuterated. We can attribute this increase in reactivity primarily to quantum tunneling effects since most ZPE effects are accounted for in the ΔG⧧ calculations. The effect of temperature on the normalized transmission coefficient for quantum proton transfer is shown in Figure 6.
Figure 6. Normalized transmission coefficient for the quantum proton which originates at the H3O+ transition state. Here we show results from all temperatures from 280 to 300 K. There is no clear relationship between the transmission coefficients and the temperature. This is largely due to the possibilities encountered by recrossing trajectories, and also small differences in the location of the dividing surface as determined by ΔG⧧ calculations.
There is no clear correlation between the temperature and the magnitude of the normalized transmission coefficient, and there are several reasons for this. First, this transmission coefficient is normalized, and κ(T, t)/⟨v(x0, p0)⟩ is the property which is being plotted, so any changes in ⟨v(x0, p0)⟩ are not accounted for in this figure due to the predictable effect temperature has on average speed of thermalized particles. Additionally, the physics of this property are much more complex than those of E
DOI: 10.1021/acs.jpca.6b10337 J. Phys. Chem. A XXXX, XXX, XXX−XXX
Article
The Journal of Physical Chemistry A Table 1. Rate Constant Results from Dynamics Calculationsa quantum
classical
TS geom
T (K)
ΔG⧧ (mEh)
κ(T)/⟨v0⟩
k(T) (a0Eh/ℏ)
ΔG⧧ (mEh)
κ(T)/⟨v0⟩
k(T) (a0Eh/ℏ)
+
280 290 300 310 320 280 290 300 310 320 280 290 300 310 320 280 290 300 310 320
19.451794 19.463857 19.537238 19.629381 19.671434 22.171916 22.200339 22.238121 22.244514 22.094649 19.248134 19.307641 19.310548 19.415242 19.441547 22.086350 22.028968 22.029280 21.995500 22.055614
0.799177 0.748446 0.781365 0.820800 0.776110 0.481504 0.597122 0.512964 0.531412 0.467219 0.820651 0.767230 0.748294 0.824548 0.654180 0.607580 0.530144 0.532456 0.568535 0.510619
6.264924(−14) 1.256574(−13) 2.505377(−13) 4.885869(−13) 8.417656(−13) 1.262030(−15) 3.745515(−15) 7.224883(−15) 1.544796(−14) 3.384926(−14) 8.175811(−14) 1.539651(−13) 3.050870(−13) 6.123826(−13) 8.827796(−13) 1.785677(−15) 3.983323(−15) 9.455949(−15) 2.151990(−14) 3.838502(−14)
26.542126 26.535161 26.530839 26.514990 26.505857 26.542665 26.524988 26.526337 26.509018 26.500151 26.870381 26.864571 26.856133 26.844141 26.829343 26.831776 26.825912 26.819191 26.813807 26.800064
0.230357 0.238611 0.246253 0.263330 0.263707 0.226805 0.245884 0.244451 0.255352 0.274668 0.140163 0.155317 0.164957 0.173052 0.197335 0.129676 0.142541 0.159019 0.173347 0.181416
6.273344(−18) 1.866137(−17) 5.143339(−17) 1.395623(−16) 3.322309(−16) 4.398938(−18) 1.387286(−17) 3.665174(−17) 9.740193(−17) 2.493685(−16) 2.633831(−18) 8.477742(−18) 2.445591(−17) 6.557346(−17) 1.805668(−16) 1.813906(−18) 5.793291(−18) 1.750798(−17) 4.844435(−17) 1.224912(−16)
H3O (H)
H3O+ (D)
OH− (H)
OH− (D)
a Change in free energy ΔG⧧, normalized transmission coefficients κ(T)/⟨v0⟩ where κ(T) ≡ limt→∞ κ(T, t) and ⟨v0⟩ ≡ ⟨v(x0, p0)⟩, and rate constants k(T). H/D values in parentheses indicate which isotope is present. H3O+ corresponds to the 1 → 2 → 4 reaction path, and OH− corresponds to the 1 → 3 → 4 path.
calculating reaction rate constants. We have shown that, in this model system, quantum effects of proton motion at temperatures near 300 K can increase the rate constants by over 1000-fold, and this is due to a combination of zero-point energy and tunneling effects. The transmission coefficient, which is a measure of the degree of recrossing from the dividing surface, was shown to be 4 times larger in quantum calculations. This is an indicator that tunneling plays an important role in this reaction. The kH/kD kinetic isotope effect ranged from 20 to 50 in the temperature range 280−320 K, indicating that quantum effects dramatically influence the proton transfer mechanism in this model system. We would like to investigate how the number of bridge water molecules influence both rate constants and the importance of quantum effects. Since the RPMD method will scale linearly with the number of protons, we can examine much larger systems within this framework. We would like to expand our model system, especially in cases where many water molecules are present, to include a hydrophobic tunnel in which the water bridge can be placed while varying the donor and acceptor molecules. We expect that the structure of this tunnel will have a profound impact on the rigidity of the water bridge, and therefore on the importance of quantum effects in proton motion. In an effort to accurately simulate the motion of protons within this environment, we would like to explicitly treat all degrees of freedom on a high-dimensional potential energy surface, and on-the-fly electronic structure calculations may be most suitable for these cases. The self-consistent charge density functional tight binding method is a semiempirical electronic structure method which has been used to calculate the electronic structure of large biological systems51 as well as properties of bulk protonated water systems,52,53 and may therefore be suitable for exploring larger chemical systems in future studies of this phenomenon.
Figure 8. Kinetic isotope effect kH/kD for both reaction pathways when quantum effects are included. We see nearly identical behavior, with a very large KIE of 30−35 near room temperature. This indicates that deuterating the system would dramatically decrease the rate, and this is due primarily to the quantum nature of the transferring nuclei. For purely classical rate calculations, the KIE was calculated to be √2 for every case, as we expect due to the Boltzmann distribution of velocities for particles of different masses.
■
SUMMARY AND CONCLUSION Rigid, highly organized water wires in biological systems have been observed to move excess protons over very long distances in a relatively short time. An understanding of this proton transfer mechanism can lead to a more accurate description of the chemical physics at this level, and it could potentially lead to the development of molecules which interact with these biological systems, or the design of materials which take advantage of this type of proton transfer. We set out to investigate how quantum effects of proton motion influence the rate of charge transfer in biological systems such as this by studying a small model system similar to the one proposed by Kaila and Hummer.25 We account for quantum effects in the dynamics calculations by employing the ring polymer molecular dynamics method and directly F
DOI: 10.1021/acs.jpca.6b10337 J. Phys. Chem. A XXXX, XXX, XXX−XXX
Article
The Journal of Physical Chemistry A
■
upon QB Photoreduction in Rhodobacter sphaeroides Reaction Center Mutants at Asp-L213 and Glu-L212 Sites. Biochemistry 2004, 43, 7236−7243. (12) Ilan, B.; Tajkhorshid, E.; Schulten, K.; Voth, G. A. The Mechanism of Proton Exclusion in Aquaporin Channels. Proteins: Struct., Funct., Genet. 2004, 55, 223−228. (13) Burykin, A.; Warshel, A. What Really Prevents Proton Transport through Aquaporin? Charge Self-Energy Versus Proton Wire Proposals. Biophys. J. 2003, 85, 3696−3706. (14) Chakrabarti, N.; Tajkhorshid, E.; Roux, B.; Pomès, R. Molecular Basis of Proton Blockage in Aquaporins. Structure 2004, 12, 65−74. (15) Rivard, U.; Thomas, V.; Bruhacs, A.; Siwick, B.; Iftimie, R. Donor-Bridge-Acceptor Proton Transfer in Aqueous Solution. J. Phys. Chem. Lett. 2014, 5, 3200−3205. (16) Markovitch, O.; Chen, H.; Izvekov, S.; Paesani, F.; Voth, G. A.; Agmon, N. Special Pair Dance and Partner Selection: Elementary Steps in Proton Transport in Liquid Water. J. Phys. Chem. B 2008, 112, 9456−9466. (17) Hassan, S. A.; Hummer, G.; Lee, Y. Effects of Electric Fields on Proton Transport through Water Chains. J. Chem. Phys. 2006, 124, 204510. (18) Vaitheeswaran, S.; Rasaiah, J. C.; Hummer, G. Electric Field and Temperature Effects on Water in the Narrow Nonpolar Pores of Carbon Nanotubes. J. Chem. Phys. 2004, 121, 7955−7965. (19) Li, H.; Mahanti, S. D.; Pinnavaia, T. J. Water Mediated Proton Transfer in a Mesostructured Aluminosolicate Framework: An ab initio Molecular Dynamics Study. J. Phys. Chem. B 2005, 109, 21908−21914. (20) Pomes, R.; Roux, B. Structure and Dynamics of a Proton Wire: a Theoretical Study of H+ Translocation Along a Single-File Water Chain in the Gramicidin A Channel. Biophys. J. 1996, 71, 19−39. (21) Ceriotti, M.; Fang, W.; Kusalik, P. G.; McKenzie, R. H.; Michaelides, A.; Morales, M. A.; Markland, T. E. Nuclear Quantum Effects in Water and Aqueous Systems: Experiment, Theory, and Current Challenges. Chem. Rev. 2016, 116, 7529−7550. (22) Lobaugh, J.; Voth, G. A. The Quantum Dynamics of an Excess Proton in Water. J. Chem. Phys. 1996, 104, 2056−2069. (23) Pavese, M.; Chawla, S.; Lu, D.; Lobaugh, J.; Voth, G. A. Quantum Effects and the Excess Proton in Water. J. Chem. Phys. 1997, 107, 7428−7432. (24) Schmitt, U. W.; Voth, G. A. The Isotope Effect on the Hydrated Proton. Chem. Phys. Lett. 2000, 329, 36−41. (25) Kaila, V. R. I.; Hummer, G. Energetics and Dynamics of Proton Transfer Reactions Along Short Water Wires. Phys. Chem. Chem. Phys. 2011, 13, 13207−13215. (26) Shao, Y.; Molnar, L. F.; Jung, Y.; Kussmann, J.; Ochsenfeld, C.; Brown, S. T.; Gilbert, A. T.; Slipchenko, L. V.; Levchenko, S.; O’Neill, D. P.; et al. Advances in methods and algorithms in a modern quantum chemistry program package. Phys. Chem. Chem. Phys. 2006, 8, 3172− 3191. (27) Press, W. H.; Flannery, B. P.; Teukolsky, S. A.; Vetterling, W. T. Numerical Recipes: The Art of Scientific Computing, 2nd ed.; Cambridge University Press: Cambridge, 1992. (28) Topaler, M.; Makri, N. Quantum Rates for a Double Well Coupled to a Dissipative Bath: Accurate Path Integral Results and Comparison with Approximate Theories. J. Chem. Phys. 1994, 101, 7500−7519. (29) Habershon, S.; Manolopoulos, D. E.; Markland, T. E.; Miller, T. F. Ring-Polymer Molecular Dynamics: Quantum Effects in Chemical Dynamics from Classical Trajectories in an Extended Phase Space. Annu. Rev. Phys. Chem. 2013, 64, 387−413. (30) Craig, I. R.; Manolopoulos, D. E. Quantum Statistics and Classical Mechanics: Real Time Correlation Functions and Ring Polymer Molecular Dynamics. J. Chem. Phys. 2004, 121, 3368−3373. (31) Richardson, J. O.; Althorpe, S. C. Ring-Polymer Moleular Dynamics Rate-Theory in the Deep-Tunneling Regime: Connection with Semiclassical Instanton Theory. J. Chem. Phys. 2009, 131, 214106. (32) Perez, A.; Tuckerman, M. E.; Muser, M. H. A Comparative Study of the Centroid and Ring-Polymer Molecular Dynamics
ASSOCIATED CONTENT
S Supporting Information *
The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpca.6b10337. Description for use of codes and data (PDF) Necessary codes for reproducing the calculations found in this paper (ZIP) Q-Chem electronic structure data used to generate the 2500 grid points that were used in the potential energy surface (ZIP)
■
AUTHOR INFORMATION
Corresponding Author
*E-mail:
[email protected]. ORCID
James W. Mazzuca: 0000-0002-4563-6606 Notes
The authors declare no competing financial interest.
■
ACKNOWLEDGMENTS This work was supported by the Alma College Provost’s Office. This research is part of the Blue Waters sustained-petascale computing project, which is supported by the National Science Foundation (Awards OCI-0725070 and ACI-1238993) and the state of Illinois. Blue Waters is a joint effort of the University of Illinois at Urbana−Champaign and its National Center for Supercomputing Applications. We would like to thank Shodor and the National Computational Science Institute (NCSI) for organizing the Blue Waters Student Internship Program (BWSIP).
■
REFERENCES
(1) Nagle, J. F.; Morowitz, H. J. Molecular Mechanisms for Proton Transport in Membranes. Proc. Natl. Acad. Sci. U. S. A. 1978, 75, 298− 302. (2) Agmon, N. The Grotthuss Mechanism. Chem. Phys. Lett. 1995, 244, 456−462. (3) Wraight, C. A. Chance and Design - Proton Transfer in Water, Channels and Bioenergetic Proteins. Biochim. Biophys. Acta, Bioenerg. 2006, 1757, 886−912. (4) Rosenberg, P. A.; Finkelstein, A. Interaction of Ions and Water in Gramicidin A Channels. J. Gen. Physiol. 1978, 72, 327−340. (5) Pomès, R.; Roux, B. Free Energy Profiles for H+ Conduction Along Hydrogen-Bonded Chains of Water Molecules. Biophys. J. 1998, 75, 33−40. (6) Xu, J.; Voth, G. A. Computer Simulation of Explicit Proton Translocation in Cytochrome C Oxidase: The D Pathway. Proc. Natl. Acad. Sci. U. S. A. 2005, 102, 6795−6800. (7) Lee, H. J.; Svahn, E.; Swanson, J. M. J.; Lepp, H.; Voth, G. A.; Brzezinski, P.; Gennis, R. B. Intricate Role of Water in Proton Transport Through Cytochrome C Oxidase. J. Am. Chem. Soc. 2010, 132, 16225−16239. (8) Baudry, J.; Tajkhorshid, E.; Molnar, F.; Phillips, J.; Schulten, K. Molecular Dynamics Study of Bacteriorhodopsin and the Purple Membrane. J. Phys. Chem. B 2001, 105, 905−918. (9) Kandt, C.; Schlitter, J.; Gerwert, K. Dynamics of Water Molecules in the Bacteriorhodopsin Trimer in Explicit Lipid/Water Environment. Biophys. J. 2004, 86, 705−717. (10) Maupin, C. M.; Castillo, N.; Taraphder, S.; Tu, C.; McKenna, R.; Silverman, D. N.; Voth, G. A. Chemical Rescue of Enzymes: Proton Transfer in Mutants of Human Carbonic Anhydrase II. J. Am. Chem. Soc. 2011, 133, 6223−6234. (11) Nabedryk, E.; Breton, J.; Okamura, M. Y.; Paddock, M. L. Identification of a Novel Protonation Pattern for Carboxylic Acids G
DOI: 10.1021/acs.jpca.6b10337 J. Phys. Chem. A XXXX, XXX, XXX−XXX
Article
The Journal of Physical Chemistry A Methods for Approximating Quantum Time Correlation Functions from Path Integrals. J. Chem. Phys. 2009, 130, 184105. (33) Miller, W. H.; Schwartz, S. D.; Tromp, J. W. Quantum Mechanical Rate Constants for Bimolecular Reactions. J. Chem. Phys. 1983, 79, 4889−4898. (34) Craig, I. R.; Manolopoulos, D. E. Chemical Reaction Rates from Ring Polymer Molecular Dynamics. J. Chem. Phys. 2005, 122, 084106. (35) Craig, I. R.; Manolopoulos, D. E. A Refined Ring Polymer Molecular Dynamics Theory of Chemical Reaction Rates. J. Chem. Phys. 2005, 123, 034102. (36) Collepardo-Guevara, R.; Suleimanov, Y. V.; Manolopoulos, D. E. Bimolecular Reaction Rates from Ring Polymer Molecular Dynamics. J. Chem. Phys. 2009, 130, 174713. (37) Suleimanov, Y. V.; Collepardo-Guevara, R.; Manolopoulos, D. E. Bimolecular Reaction Rates from Ring Polymer Molecular Dynamics: Applicaiton to H + CH4 → H2 + CH3. J. Chem. Phys. 2011, 134, 044131. (38) Perez de Tudela, R. P.; Aoiz, F. J.; Suleimanov, Y. V.; Manolopoulos, D. E. Chemical Reaction Rates from Ring Polymer Molecular Dynamics: Zero Point Energy Conservation in Mu + H2 → MuH + H. J. Phys. Chem. Lett. 2012, 3, 493−497. (39) Li, Y.; Suleimanov, Y. V.; Green, W. H.; Guo, H. Quantum Rate Coefficients and Kinetic Isotope Effect for the Reaction Cl + CH4 → HCl + CH3 from Ring Polymer Molecular Dynamics. J. Phys. Chem. A 2014, 118, 1989−1996. (40) Wong, K. F.; Sonnenberg, J. L.; Paesani, F.; Yamamoto, T.; Vanicek, J.; Zhang, W.; Schlegel, H. B.; Case, D. A.; Cheatham, T. E.; Miller, W. H.; et al. Proton Transfer Studied Using a Combined Ab Initio Reactive Potential Surface with Quantum Path Integral Methodology. J. Chem. Theory Comput. 2010, 6, 2566−2580. (41) Marsalek, O.; Markland, T. E. Ab initio Molecular Dynamics with Nuclear Quantum Effects at Classical Cost: Ring Polymer Contraction for Density Functional Theory. J. Chem. Phys. 2016, 144, 054112. (42) Frenkel, D., Smit, B. Understanding Molecular Simulation: From Algorithms to Applications, 2nd ed.; Academic Press: San Diego, CA, 2002. (43) Chandler, D. Statistical Mechanics of Isomerization Dynamics in Liquids and the Transition State Approximation. J. Chem. Phys. 1978, 68, 2959−2970. (44) Voth, G. A.; Chandler, D.; Miller, W. H. Rigorous Formulation of Quantum Transition State Theory and its Dynamical Corrections. J. Chem. Phys. 1989, 91, 7749−7760. (45) Kästner, J.; Thiel, W. Bridging the Gap Between Thermodynamic Integration and Umbrella Sampling Provides a Novel Analysis Method: “Umbrella Integration. J. Chem. Phys. 2005, 123, 144104. (46) Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; DiNola, A.; Haak, J. R. Molecular Dynamics with Coupling to an External Bath. J. Chem. Phys. 1984, 81, 3684−3690. (47) Rickert, K. W.; Klinman, J. P. Nature of Hydrogen Transfer in Soybean Lipoxygenase 1: Separation of Primary and Secondary Isotope Effects. Biochemistry 1999, 38, 12218−12228. (48) Knapp, M. J.; Rickert, K.; Klinman, J. P. TemperatureDependent Isotope Effects in Soybean Lipoxygenase-1: Correlating Hydrogen Tunneling with Protein Dynamics. J. Am. Chem. Soc. 2002, 124, 3865−3874. (49) Klinman, J. P. An Integrated Model for Enzyme Catalysis Emerges from Studies of Hydrogen Tunneling. Chem. Phys. Lett. 2009, 471, 179−193. (50) Basran, J.; Sutcliffe, M. J.; Scrutton, N. S. Enzymatic H-Transfer Requires Vibration-Driven Extreme Tunneling. Biochemistry 1999, 38, 3218−3222. (51) Elstner, M. The SCC-DFTB Method and its Application to Biological Systems. Theor. Chem. Acc. 2006, 116, 316−325. (52) Maupin, C. M.; Aradi, B.; Voth, G. A. The Self-Consistent Charge Density Functional Tight Binding Method Applied to Liquid Water and the Hydrated Excess Proton: Benchmark Simulations. J. Phys. Chem. B 2010, 114, 6922−6931.
(53) Goyal, P.; Elstner, M.; Cui, Q. Application of the SCC-DFTB Method to Neutral and Protonated Water Clusters and Bulk Water. J. Phys. Chem. B 2011, 115, 6790−6805.
H
DOI: 10.1021/acs.jpca.6b10337 J. Phys. Chem. A XXXX, XXX, XXX−XXX
