Critical Role of Deep Hydrogen Tunneling to Accelerate the

Article pubs.acs.org/JPCB

Critical Role of Deep Hydrogen Tunneling to Accelerate the Antioxidant Reaction of Ubiquinol and Vitamin E Taichi Inagaki and Takeshi Yamamoto* Department of Chemistry, Graduate School of Science, Kyoto University, Kyoto 606-8502, Japan S Supporting Information *

ABSTRACT: In biomembranes a variety of antioxidants work to suppress oxidative damage. Vitamin E and ubiquinol are among the most important lipid-soluble antioxidants, which trap lipid peroxyl radicals directly or work cooperatively in the regeneration of vitamin E radicals by ubiquinol. Here, we investigate the latter regeneration reaction by using variational transition-state theory with multidimensional tunneling corrections. The result shows that the system forms a compact H-bonded complex by significantly rearranging the donor and acceptor moieties, which leads to a rather narrow potential barrier for H transfer and a very large tunneling effect with a transmission coefficient >4000. In accord with experiment, the Arrhenius activation energy is found to be very small (∼1 kcal/mol), which is interpreted here in terms of mean tunneling energy through the barrier. Regarding the electronic structure, we demonstrate that the present reaction proceeds via a proton-coupled electron transfer (PCET) mechanism and suggest that the PCET character also contributes to the large tunneling effect by sharpening the potential barrier. Finally, a systematic comparison is made among relevant reactions and it is indicated that the antioxidant defense of biomembranes may benefit rather significantly from quantum tunneling to enhance the reaction efficiency.

1. INTRODUCTION Living organisms are exposed to severe oxidative stress due to free radicals, which may cause damages in biomembranes and DNA and eventually lead to various diseases.1−4 Living systems are thus equipped with a variety of antioxidant mechanisms,5 particularly against lipid peroxidation (namely, the oxidative degradation of biomembranes).6,7 Once a lipid radical (L·) has been generated by heat, light, metals, or other radicals, lipid peroxidation proceeds as follows: L· + O2 → LOO·

(1)

LOO· + LH → LOOH + L·

(2)

Figure 1. Chemical structure of (a) ubiquinol and (b) vitamin E (tocopheroxyl) radical. This paper addresses the case R1 = R2 = isopropyl group (i.e., 5,7-diisopropyl-tocopheroxyl radical) rather than R1 = R2 = methyl group (i.e., α-tocopheroxyl radical) to make quantitative comparison with experiment (refs 37 and 38). The isoprenoid chain of ubiquinol and the phytil chain of tocopherol are replaced by methyl groups for computational efficiency.

where the lipid radical L· reacts extremely rapidly with oxygen to form a lipid peroxyl radical LOO·, which in turn attacks another lipid molecule (LH) to produce lipid peroxide (LOOH) and a new radical L·. The rates of reactions 1 and 2 are typically on the order of 109 and 101 M−1 s−1, respectively, i.e., reaction 2 is much slower than reaction 1.7 The above chain reactions transform a number of lipid molecules to lipid peroxide, causing degradation of lipid membranes. The above reaction terminates when the chain-propagating LOO· radicals are trapped (or scavenged) by antioxidants such as vitamins, ubiquinol, flavonoids, and so on. Among others, vitamin E is known to be the most important lipid-soluble antioxidant, with α-tocopherol being the most abundant and active form in humans (see Figure 1).8−12 Vitamin E (tocopherol, TocH) stops the chain reaction by reducing the lipid peroxyl radical as follows: © XXXX American Chemical Society

TocH + LOO· → Toc· + LOOH

(3)

where Toc· is a tocopheroxyl radical. It should be noted that because TocH is a phenolic compound, the Toc· radical is relatively stable because of the resonance effect. However, it Received: October 16, 2013 Revised: December 24, 2013

A

dx.doi.org/10.1021/jp410263f | J. Phys. Chem. B XXXX, XXX, XXX−XXX

The Journal of Physical Chemistry B

Article

KIE value of 22. Antioxidant reactions may proceed via concerted hydrogen transfer or sequential electron/proton (or proton/electron) transfer mechanism.43 However, because several previous studies suggest that the sequential mechanism is favored in aqueous solution at proper pH,31−35 here we concentrate only on the concerted H-transfer mechanism. With this assumption, we study the kinetics of reaction 4 by combined use of variational transition-state theory with multidimensional tunneling corrections and integral equation theory for molecular solvents. The obtained results are then used to obtain insights into the condition under which quantum tunneling plays a crucial role in the antioxidant efficiency. Because ubiquinol and vitamin E are phenolic compounds, it is also of interest whether the reaction proceeds via a hydrogen atom transfer (HAT) or proton-coupled electron transfer (PCET) mechanism.44−51 This problem is addressed by performing ground- and excited-state calculations and examining the electronic character of the system. In the following, we denote the complex of ubiquinol and vitamin E derivative in Figure 1 simply as the UQ/Toc complex, while more specific names such as those in eq 4 are used where necessary.

may eventually react with lipid molecules and continue or even enhance the chain reaction (called the prooxidant action).11,13 To suppress this, the Toc· radical is recycled back to the active reduced form (TocH) through reduction by other antioxidants, such as ubiquinol and vitamin C.10,11 The regeneration of Toc· by ubiquinol proceeds as follows: UQH 2 + Toc· → UQH· + TocH

(4)

where UQH2 is the doubly reduced form of ubiquinone (see Figure 1), which is an important lipid-soluble antioxidant synthesized in the body.14 The UQH· radical thus produced is recycled back to UQH2 by different redox mechanisms.10,15 Ubiquinol thus contributes to the antioxidant defense of biomembranes in two ways, i.e., by directly trapping lipid peroxyl radicals (as in eq 3) or by regenerating vitamin E radicals to the active reduced form.10,13,16−23 A variety of theoretical studies have been performed to obtain mechanistic insights into antioxidant reactions.24−36 Structure−activity relationships have been studied extensively based on quantum chemical calculations. In some cases, reaction kinetics have been examined to identify the most favorable reaction mechanism. Because many of the antioxidant reactions involve a hydrogen transfer step, a natural question is whether nuclear quantum effects play some role in the antioxidant efficiency. As expected, it is found that the degree of quantum effects varies considerably depending on systems. For example, a recent theoretical study has demonstrated that the trapping of methyl peroxyl radical by epicatechin (i.e., a flavonoid rich in green tea) is dramatically accelerated by quantum tunneling28 with the transmission coefficient being 103−105. Large tunneling effects have also been obtained for the antioxidant reaction of quercetin (i.e., a flavonoid rich in fruits and vegetables) which exhibits a 500-fold increase in the reaction rate.29 On the other hand, it is also known that quantum effects are not particularly large for other systems; for example, the trapping of an OH or OOH radical by vitamin E or ubiquinol exhibits minor or essentially no tunneling effect.24−27 That is, the reaction proceeds via the traditional “over the barrier” mechanism. Only a limited number of experimental studies have been reported to date concerning nuclear quantum effects in the antioxidant reactions. Of particular interest here is the work of Nagaoka, Mukai, and co-workers,37−41 who investigated a series of model antioxidant reactions in solution and micellar dispersion. They demonstrated via stopped-flow spectroscopy that the regeneration reaction of vitamin E by ubiquinol exhibits a large deuterium kinetic isotope effect (KIE) of ∼20, which exceeds the maximum semiclassical ratio,42 and that the Arrhenius activation energy is about 1−3 kcal/mol, which is fairly small compared to typical barrier heights of hydrogen transfer reactions. Large KIE values have also been obtained for the antioxidant reaction of flavonoids in solution.40 On the other hand, the KIE is found to be relatively small for the regeneration of vitamin E by vitamin C (ascorbate monoanion).41 This indicates again that the role of quantum effects are rather different among relevant antioxidant reactions. On the basis of the above observation, here we investigate the reaction mechanism of the regeneration of vitamin E radical by ubiquinol in eq 4 with particular focus on nuclear tunneling effects. To make a quantitative comparison with experiment,37,38 we consider the regeneration reaction of a derivative of vitamin E radical by ubiquinol in ethanol (see Figure 1 for chemical structures), which has been shown to exhibit a large

2. COMPUTATIONAL METHODS The gas-phase potential energy surface for this reaction was first calculated by using density functional theory (DFT). The unrestricted MPW1K functional52 was used for this purpose (unless otherwise noted), which has been developed with the aim of accurately describing the barrier height of H-transfer reactions. It has been indicated that the MPW1K method gives energy profiles for typical PCET reactions that are more accurate than the traditional B3LYP method.45 The basis set employed was of double-ζ plus polarization quality ((10s5p1d/ 4s1p)/[3s2p1d/2s1p]) with the polarization functions on nonpolar hydrogens removed, which led to a total of 579 basis functions. Geometry optimization was performed to locate the equilibrium and transition-state geometries. We have confirmed via normal-mode analysis that the number of imaginary frequency is 0 and 1 for equilibrium and transitionstate geometries, respectively. The minimum energy path was calculated in mass-scaled Cartesian coordinates by starting from the saddle point and going downhill to the reactant or product side. For this purpose, the Gonzalez−Schlegel method53 was used with a step size of 0.02 bohr. To examine the PCET character of the system, we also calculated the excited-state potential energy profile using multiconfiguration wave function theories. Specifically, the state-averaged complete active space self-consistent field (SACASSCF) method was used to calculate the ground and excited states with equal weights. Dynamic electron correlation was taken into account by using multiconfiguration quasidegenerate perturbation theory (MC-QDPT).54,55 The active space was constructed by distributing 11 electrons over 10 orbitals. More details of the excited-state calculations including the choice of active space are given in the Supporting Information. Because the obtained results suggest that the energy gap between the ground and excited states is relatively large (∼14 kcal/mol, see PCET Character and Electronic Nonadiabaticity), we performed the subsequent reaction rate calculation using an adiabatic formulation of transition-state theory (see Reaction Rates and Kinetic Isotope Effect). Next, the free-energy surface in ethanol solution was calculated using the three-dimensional reference interaction B



Article

site model self-consistent field (3D-RISM-SCF) method.56,57 This method calculates the Helmholtz free energy of the solvated molecule with a geometry R as follows: G(R) = ⟨Ψ|H0|Ψ⟩ + Δμ(R)

3. RESULTS AND DISCUSSION 3.1. Reaction Profile in the Gas Phase and in Solution. We first consider reaction 4 in the gas phase. The reaction proceeds by first making a reactant complex (RC) in the entrance channel, followed by a transition state (TS), a product complex (PC) in the product channel, and finally leads to the separated products. The optimized geometries obtained from the MPW1K calculation are shown in Figure 2, and selected

(5)

where the first term is the solute electronic energy and the second term is the solvation free energy. In this method, the Δμ was calculated by solving the 3D-RISM integral equation57 together with the Kovalenko−Hirata closure relation58 at a temperature of 298 K and a solvent density of 0.789 g/cm3 (which correspond to the experimental conditions).37−39 The ethanol molecule was described with a four-site united-atom model,59 while the solvent distribution functions were expressed in terms of 128 × 128 × 128 grid points spanned over a cubic simulation box of side length 38.4 Å. We confirmed that the above simulation box is sufficient for converging the solvation free energy as well as solvent distribution functions. As in the case of the gas-phase reaction, the equilibrium and TS geometries were located on the free-energy surface G(R), and the minimum free-energy path was calculated by starting from the TS geometry. Generalized normal-mode analyses were then performed for every two points along the minimum free-energy path to obtain transversal normal modes. The rate constant for the solution reaction was calculated using variational transition-state theory with multidimensional tunneling (VTST/MT) method.60−68 In this method, the absolute reaction rate at temperature T is calculated as

k(T ) = κ(T )kCVT(T )

(6)

where kCVT(T) is the canonical variational theory (CVT) rate constant evaluated under the equilibrium solvation model and κ(T) is the semiclassical transmission coefficient that accounts for multidimensional tunneling. The κ(T) is obtained as κ (T ) =

∫ dE P(E) exp(−E/kBT ) ∫ dE PC(E) exp(−E/kBT )

(7)

where kB is the Boltzmann constant, P(E) the quantummechanical transmission probability at energy E, and PC(E) the classical counterpart of P(E). The transmission probability was calculated with several different approximations, including the small-curvature tunneling (SCT) approximation,64 largecurvature tunneling (LCT) approximation (here evaluated with the LCG4 scheme),66,68 and the microcanonical optimized multidimensional tunneling (OMT) approximation.67 Note that in the OMT method the P(E) is evaluated by accepting whichever LCT or SCT scheme gives the larger tunneling probability at each energy, i.e., P(E) = max{PLCT(E), PSCT(E)}. The transmission coefficient was calculated by using the freeenergy surface in eq 5 as an effective potential for H tunneling (called the zeroth-order canonical mean-shape approximation).69 The standard-state partition function in the CVT rate constant was evaluated using the separable rotation/vibration approximation. All the VTST/MT calculations were performed with an in-house program, which was validated by comparison to published results.61,64,66 The electronic structure calculation was performed with a locally modified version of the GAMESS program package70,71 in which the 3D-RISM-SCF method has been implemented.

Figure 2. Optimized geometries in the gas phase obtained with the MPW1K method: (a) reactant complex, (b) transition state, and (c) product complex.

geometrical parameters are listed in Table 1. A notable feature of Figure 2 is that the molecular planes of the UQ and Toc moieties are nearly orthogonal to each other so as to avoid the steric hindrance of side chains. This is qualitatively different from the phenol/phenoxyl complex (i.e., a prototypical PCET system), which has no side chains and takes a nearly coplanar TS geometry.44 To characterize the orientation of the UQ and Toc moieties, it is useful to introduce two dihedral angles as follows (see Figure 2a for atomic labels):

C

ϕ = ∠(C2−C−C′−C′2)

(8)

ω = ∠(O2 −O−O′−O′2 )

(9)



Article

Table 1. Selected Geometrical Parameters of the Reactant Complex (RC), Transition State (TS), and Product Complex (PC) in the Gas Phase and Ethanol Solutiona gas H−O H−O′ O−O′ O−H−O′ ϕ ω

ethanol

RC

TS

PC

RC

TS

PC

0.966 1.899 2.773 149.1 83.0 10.0

1.203 1.146 2.343 171.6 101.7 38.7

1.808 0.969 2.775 176.1 75.5 47.7

0.968 1.933 2.767 142.9 123.0 −42.3

1.210 1.142 2.345 171.2 138.5 13.1

1.744 0.971 2.698 166.8 143.6 13.9

a See Figure 2 for atomic labels. Torsional angles are defined as ϕ = C2−C−C′−C2′ and ω = O2−O−O′−O2′ . Bond lengths are in angstroms and angles are in degrees.

Figure 3. Energy diagram (in kcal/mol) for the reaction in the gas phase and in solution obtained with the MPW1K method. The free energy in solution was obtained from the MPW1K/3D-RISM calculation. Figures in parentheses include zero-point energies.

The physical interpretation of these variables is illustrated in Figure S1 of the Supporting Information. Qualitatively, ϕ specifies the extent to which the molecular planes of UQ and Toc are twisted around the donor−acceptor O···H···O axis, with ϕ = 90° corresponding to the orthogonal limit. On the other hand, ω describes to what extent the UQ is tilted toward one of the isopropyl side chains of Toc, with ω = 0° corresponding to the bisector of two isopropyl groups. With this definition the reactant complex assumes ϕ = 83° and ω = 10°, indicating that UQ and Toc are nearly orthogonal to each other (Figure 2a). Note that the donor hydrogen atom forms two H bonds with the methoxyl oxygen of UQ and the acceptor oxygen of Toc. That is, the transferring H atom is attracted in two directions, which makes the H transfer energetically unfavorable at this configuration. The system then makes a substantial rearrangement of UQ and Toc to reach the transition state with ϕ = 101° and ω = 38° (Figure 2b). This rearrangement of UQ and Toc plays the role of breaking the intramolecular H bond (with the methoxyl oxygen) while increasing the linearity of the O···H···O bond. Once the H transfer has completed, the system makes a further rearrangement of molecular planes and reaches the product complex with ϕ = 75° and ω = 47° (Figure 2c). Another important feature of Figure 2 is that the O···H···O bond is rather compact. Indeed, the donor−acceptor distance at the transition state (2.34 Å) is similar to that of the phenol/ phenoxyl complex (2.40 Å).48 This is somewhat unexpected because of the apparently large steric hindrance of side chains. In fact, the UQ/Toc complex avoids the steric hindrance by significantly rotating the molecular planes while making a very compact H bond between the donor and acceptor. The energy profile in the gas phase that corresponds to the optimized geometries is shown in Figure 3. The binding energy of the reactant complex (as measured from the separated reactants) is 4.9 kcal/mol, which is essentially attributed to the donor−acceptor H bond. The barrier height from the separated reactants is 12.1 kcal/mol, while the intrinsic reaction barrier (i.e., the energy difference between the transition state and the reactant complex) is 17.0 kcal/mol. The MPW1K calculation predicts the reaction to be slightly endothermic in the gas phase; the reaction energy is estimated to be 4.3 kcal/mol. This is consistent with the experimental values of bond dissociation enthalpy for the reactive OH bonds, which are 80.5 ± 1.5 and 78.2 ± 0.25 kcal/mol for UQH2 and TocH in benzene solution, respectively.72,73 We next consider the reaction in ethanol solution. Figure S3 of the Supporting Information displays the optimized geo-

metries in solution obtained from the 3D-RISM-SCF calculation, with selected geometrical parameters listed in Table 1. We find that the optimized geometries in solution are qualitatively similar to those in the gas phase. The UQ and Toc moieties form an H-bonded complex of significantly distorted form, with a rather compact O···H···O bond (Table 1). The complex is somewhat more distorted in solution than in the gas phase (e.g., ϕ = 123° and ω = −42° at the reactant complex) so as to increase the solvation of the methoxyl groups. Figure 3 displays the free-energy profile in solution obtained from the 3D-RISM-SCF calculation. The free-energy barrier in solution (as measured from the separated reactants) is 3−4 kcal/mol higher than the potential energy barrier in the gas phase. The reason for the increased barrier can be understood from the solute−solvent radial distribution functions (Figure S4 of the Supporting Information). This figure shows that the donor hydroxyl group and the acceptor oxygen atom are solvated more strongly for the separated reactants than for the reactant complex or the transition state, thus lowering the (relative) free energy of the separated reactants. The same conclusion can be obtained from the energy decomposition analysis of G(R) in eq 5 into the solute and solvent terms (see Figure S5 of the Supporting Information and its caption for details). Figure S5 shows that the solvation free energy Δμ of the separated reactants is about 3 kcal/mol smaller than that of the other states, resulting in an increased free-energy barrier in solution. 3.2. PCET Character and Electronic Nonadiabaticity. H-transfer reactions are sometimes categorized into hydrogen atom transfer (HAT) and proton-coupled electron transfer (PCET) reactions depending on the electronic character of the system.44−51 Prototypical examples of HAT and PCET reactions are the self-exchange reaction of the toluene/benzyl (PhCH3/PhCH2·) and phenol/phenoxyl (PhOH/PhO·) complexes, respectively. While the precise definition of the term PCET varies in the literature,44,48−50 we employ the broad definition that it refers to simultaneous transfer of a proton and an electron between different sets of orbitals.44 In the following we examine the PCET character of the system using several different diagnostics. Mayer et al.44 proposed the distinction between the HAT and PCET mechanisms based on the character of the singly occupied molecular orbital (SOMO) at the transition state. According to their criterion, the H transfer of the phenol/ phenoxyl system is classified as PCET because the SOMO is D



Article

PCET mechanisms. As for the TS geometry, the donor− acceptor distance of the UQ/Toc system (2.35 Å) is closer to that of the phenol/phenoxyl system (2.40 Å) than that of the toluene/benzyl system (2.72 Å) by reflecting the donor− acceptor hydrogen bond. PCET reactions have also been characterized in terms of electronic nonadiabaticity.48−50 To do so, we calculated the potential energy profiles for the ground and excited states in the transition-state region. Figure 5 displays the gas-phase potential

delocalized over the donor and acceptor moieties, whereas the H transfer of the toluene/benzyl system is classified as HAT because the SOMO is localized in the reactive C···H···C bond. Figure 4 displays the SOMO of the UQ/Toc system in the gas

Figure 5. Ground- and excited-state potential energy profiles in the TS region obtained from the SA-CASSCF and MC-QDPT calculations. The active space was chosen as (11e,10o). The abscissa is the Htransfer coordinate xH in eq 10. The dashed lines represent the diabatic potential curves obtained by fitting the 2 × 2 model Hamiltonian to the energy profiles at the MC-QDPT level. Zero of energy is set at the TS geometry (xH = 0) in the ground state.

energy profiles obtained from the SA-CASSCF and MC-QDPT calculations with the active space of (11e,10o). The energy profile is plotted as a function of one-dimensional H-transfer coordinate, xH, which specifies the position of the transferring H atom as follows: r(H) = r⧧(H) + x H Figure 4. Singly occupied molecular orbitals (SOMOs) of the UQ/ Toc system in the gas phase obtained from the MPW1K calculation: (a) reactant complex, (b) transition state, and (c) product complex. Note that the SOMO is localized on the Toc side in the reactant state, whereas it is localized on the UQ side in the product state. At the transition state, the SOMO is delocalized over both the UQ and Toc moieties, which is a signature of PCET reaction. The highest occupied molecular orbitals (HOMOs) are shown in the Supporting Information.

r⧧(O′) − r⧧(O) |r⧧(O′) − r⧧(O)|

(10)

⧧

where r (X) denotes the Cartesian coordinates of atom X at the transition state, while O and O′ are the donor and acceptor oxygen atoms, respectively. Note that the potential energy profiles are obtained by moving only the H atom while keeping the other heavy atoms fixed. It is seen from Figure 5 that the adiabatic energy gap between the ground and excited states is 12.2 and 14.4 kcal/mol at the SA-CASSCF and MC-QDPT levels, respectively. The energy gap for the UQ/Toc system is thus intermediate between those of the phenol/phenoxyl and toluene/benzyl systems (4 and 80 kcal/mol, respectively).48 The relatively large energy gap of the UQ/Toc system (as compared to the phenol/phenoxyl system) is probably due to the distorted conformation of the transition state. This in turn suggests that the 2p orbitals of the donor−acceptor oxygens possess an overlap for the UQ/Toc system that is somewhat greater than the overlap for the phenol/phenoxyl system. The energy profiles in Figure 5 allow us to calculate the proton adiabaticity parameter pad introduced by Skone et al.48 This parameter characterizes the degree of electronic adiabaticity of a given H-transfer reaction as follows: τp pad = τe (11)

phase calculated at the MPW1K level. As seen, the SOMO at the transition state is significantly delocalized over the donor and acceptor, indicating that the H transfer proceeds via a PCET mechanism (more specifically, the proton transfer occurs between the donor and acceptor oxygens, while at the same time the electron transfer occurs between the 2p−π orbitals of UQ and Toc). The qualitative pattern of the SOMO is similar to that of the phenol/phenoxyl complex, except that the overlap between the 2p orbitals of the donor−acceptor oxygens has some σ-type character. This is because the molecular planes of UQ and Toc are distorted substantially from the coplanar limit, and as a result the O···H···O bond makes some angle with those planes. The above observation suggests that the present system may possess a mixed character between HAT and E



Article

where τp is the characteristic time spent by the tunneling proton in the barrier region (i.e., the proton tunneling time), and τe is the characteristic time required to change the electronic state (i.e., the electronic transition time). When τp ≫ τe, the proton tunneling is considered as electronically adiabatic because the electrons are able to respond virtually instantaneously to the proton motion. When τp ≪ τe, the proton tunneling is regarded as electronically nonadiabatic because the electrons are unable to rearrange fast enough to follow the proton motion. Following ref 48, we have calculated the value of τp and τe by fitting a 2 × 2 model Hamiltonian to the ground and excited potential energy curves in Figure 5 and utilizing the relevant semiclassical formulas. The obtained results are τp = 0.39 fs, τe = 2.11 fs, and hence pad = 0.19 (see Table 2). Because Figure 6. Population of the transferring electron as a function of the H coordinate xH in eq 10. It is seen that the population varies more rapidly at xH ≃ 0 for the UQ/Toc and phenol/phenoxyl complexes than for the toluene/benzyl complex by reflecting the PCET character of the former.

Table 2. Adiabaticity Parameter pad, Proton Tunneling Time τp, Electronic Transition Time τe, and Electronic Energy Gap ΔE for the UQ/Toc System Calculated at the MCQDPT(11e,10o) Levela UQ/Toc phenol/phenoxyl toluene/benzyl

pad

τp (fs)

τe (fs)

ΔE (kcal/mol)

0.19 0.01 3.45

0.39 0.10 1.28

2.11 7.60 0.37

14.4 4.0 81.7

is well-correlated with the energy gap at the TS for the individual systems. That is, a system with smaller energy gap exhibits more rapid population transfer in the TS region. All the above diagnostics suggest that the present UQ/Toc reaction may be regarded as a PCET reaction with some electronic nonadiabaticity. On the other hand, the energy gap at the TS geometry is relatively large (14 kcal/mol at the MCQDPT level), indicating that the effect of electronic nonadiabaticity on the reaction rate may be rather moderate. Therefore, in the following sections we utilize an adiabatic formulation of transition-state theory to approximately calculate the reaction rate. Including the effect of electronic nonadiabaticity remains the subject of future study.74 3.3. Reaction Rates and Kinetic Isotope Effect. In this section we calculate the rate constant for reaction 4 using variational transition-state theory with multidimensional tunneling corrections. To this end, we first calculated the minimum energy path on the free-energy surface G(R) by using the 3D-RISM-SCF method at the MPW1K level. Figure 7 displays the vibrationally adiabatic free-energy profile along the minimum energy path,75 i.e.

a

For comparison, the reference values for the phenol/phenoxyl and toluene/benzyl systems (ref 48) are also shown. The energy gap for the phenol/phenoxyl and toluene/benzyl systems are assumed to be twice the magnitude of electronic coupling reported in ref 48.

pad is smaller than 1, electronic nonadiabaticiy may be nonnegligible for the present UQ/Toc system. For comparison, the value of pad is 0.01 and 3.45 for the H transfer of the phenol/phenoxyl and toluene/benzyl complexes, respectively. The above result suggests that the UQ/Toc system is electronically intermediate between the adiabatic and nonadiabatic limits, in accord with the picture obtained from the energy gap. As an additional measure of the electronic character, we have calculated the population of the transferring electron as a function of xH. Specifically, we calculated the population of πelectron residing on the donor or acceptor moiety as follows: PM(x H) =

∑ μ∈M

(PS)μμ

(M = D, A) (12)

where P and S are the density matrix constructed from the highest occupied α molecular orbital and the overlap matrix at the MPW1K level, respectively. The index μ runs over the atomic orbitals belonging to the donor (M = D) or acceptor (M = A) moiety. Figure 6 displays the PM(xH) thus obtained for the UQ/Toc system. When the proton is bonded to the donor oxygen atom (xH < −0.3 Å), the π electron resides essentially on the UQ side with PD(xH) ≃ 1. When the proton is bonded to the acceptor oxygen atom (xH > 0.3 Å), the π electron is transferred almost completely to the Toc side with PA(xH) ≃ 1. Note that the electron population of the UQ/Toc system varies very rapidly when the proton passes the TS region, i.e., −0.1 < xH < 0.1 Å. For comparison, we also calculated the electron population for the phenol/phenoxyl and toluene/benzyl systems using the same computational protocol. It is seen that the PM(xH) for the phenol/phenoxyl system varies rapidly than the UQ/Toc system, while the PM(xH) for the toluene/ benzyl system varies much more slowly by reflecting the greater electronic adiabaticity. Notably, the slope of PM(xH) at xH ≃ 0

Figure 7. Vibrationally adiabatic free-energy profile for the reaction in solution obtained from the MPW1K/3D-RISM calculation. The abscissa s is the signed distance along the minimum free-energy path. This figure depicts the relative value of ga(s) in eq 13 as measured from the reactant complex (RC), namely ga(s) − ga(RC). The horizontal lines indicate the mean tunneling energy. F



Article

Table 3. Rate Constants (M−1 s−1) for the Reaction in Solution Obtained from the VTST/MT Calculation at the MPW1K Levela

a

T (K)

kTST

kCVT

kCVT/ZCT

288 298 308

3.27(−4) 7.12(−4) 1.48(−3)

3.20(−4) 6.97(−4) 1.45(−3)

3.49(−1) 5.33(−1) 8.00(−1)

288 298 308

4.92(−5) 1.13(−4) 2.48(−4)

4.80(−5) 1.11(−4) 2.42(−4)

3.39(−2) 5.41(−2) 8.37(−2)

kCVT/SCT H-transfer 1.43 2.06 2.91 D-transfer 1.17(−1) 1.73(−1) 2.52(−1)

kCVT/LCT

kCVT/OMT

1.72 2.40 3.30

2.06 2.89 4.00

3.70(4) 4.11(4) 4.93(4)

4.94(−2) 7.44(−2) 1.09(−1)

1.19(−1) 1.76(−1) 2.57(−1)

1.73(3) 1.90(3) 2.02(3)

expt

Figures in parentheses indicate powers of 10. The experimental values are obtained from ref 37.

ga (s) = GMEP(s) +

∑ εi(s) i

Table 4. Transmission Coefficient for the Reaction in Solution Obtained from the VTST/MT Calculation at the MPW1K Levela

(13)

where s is the signed distance along the minimum energy path (with s = 0 corresponding to the transition state), and {εi(s)} are zero-point energies for orthogonal vibrational modes. It should be noted that zero of energy for the free-energy profile in Figure 7 is set at the reactant complex. The sharp barrier in the region −0.5 < s < 0.5 corresponds to the local migration of the H atom between the donor and acceptor oxygens. On the other hand, the gradual decrease in ga(s) in the region |s| > 0.5 corresponds to the slow rearrangement of the donor and acceptor moieties. Figure S6a of the Supporting Information depicts the variation of the donor−acceptor distance r(O−O′) against the reaction coordinate r(O−H) − r(H−O′), which shows that the local migration of the H atom is well-separated from the donor−acceptor motions. Figure S6b of the Supporting Information depicts the variation of other internal coordinates along the minimum energy path. This figure shows that the system first makes a large change in the dihedral angles ϕ and ω, which corresponds to the rearrangement of UQ and Toc. The system then evolves by compressing the donor− acceptor distance to reach a precursor state at s ≃ −0.5 (called the “tunneling-ready” configuration).76−78 Subsequently, the local migration of the H atom occurs with the other heavy atoms almost fixed. With the free-energy profile thus obtained, we calculated the TST and CVT rate constants, kTST and kCVT, which are listed in Table 3. The difference between kTST and kCVT is quite small, i.e., the variational effect is almost negligible for the present reaction. The deuterium kinetic isotope effect (KIE) obtained from the CVT rate constant is 6.3 at 298 K, which is much smaller than the corresponding experimental value, 21.6. We emphasize that the above CVT rate constant accounts for zeropoint energies but not for tunneling effects. There are several approximate ways to estimate the tunneling correction. The simplest approach is to use the zero-curvature tunneling (ZCT) approximation.62 This approximation assumes that the tunneling occurs along the minimum energy path while neglecting the curvature of the reaction path. The obtained value of κ at the ZCT level is 765 (see Table 4). In reality, it is known that the hydrogen atom exhibits the cornercutting effect, i.e., the tunneling path deviates from the minimum energy path toward the concave side on the potential energy surface. The small-curvature tunneling (SCT) approximation accounts for corner cutting by using a path-dependent effective mass and is suitable for reactions with weak to moderate tunneling. The obtained SCT value of κ is 2950, which is about four times greater than the ZCT value and

T (K)

a

κZCT

288 298 308

1.09(3) 7.65(2) 5.52(2)

288 298 308

7.06(2) 4.87(2) 3.46(2)

κSCT H-transfer 4.48(3) 2.95(3) 2.01(3) D-transfer 2.43(3) 1.56(3) 1.04(3)

κLCT

κOMT

5.39(3) 3.45(3) 2.28(3)

6.44(3) 4.15(3) 2.76(3)

1.03(3) 6.70(2) 4.52(2)

2.47(3) 1.59(3) 1.06(3)

Figures in parentheses indicate powers of 10.

points to the importance of corner cutting. Another way to evaluate κ is the large-curvature tunneling (LCT) approximation, which relies on a straight-line path connecting the reactant and product valleys and is more suitable for systems with strong corner cutting. The LCT scheme predicts κ to be 3450. We have also utilized the microcanonical optimized multidimensional tunneling (OMT) approximation, which improves upon the other schemes by choosing an optimal tunneling path between the SCT and LCT ones at each energy. The OMT scheme predicts more significant tunneling with κ = 4150. Overall, the above results suggest that the tunneling effect plays a crucial role for accelerating the present UQ/Toc reaction. Figure 8 displays the thermal distribution of tunneling probability, P(E) exp(−E/kBT), with zero of energy set at the reactant complex. We recall that the transmission coefficient κ is obtained by integrating the above distribution function over E (see eq 7). It is seen from Figure 8 that the dominant tunneling mechanism changes according to the mass of the transferring hydrogen. In the case of H transfer (Figure 8a), the LCT mechanism dominates at low energies because the cornercutting effect is significant, while the SCT mechanism becomes more important at higher energies. In the case of D transfer (Figure 8b), the SCT mechanism dominates over the entire energy region because the tunneling path remains close to the minimum energy path. The representative tunneling energy (i.e., the energy at which the distribution takes on the maximal value) is 4.1 and 5.5 kcal/mol for H and D transfer, respectively, as measured from the reactant complex. In practice, we find that the mean tunneling energy (MTE) is more relevant for the Arrhenius activation energy (see Discussion). The MTE is found to be 5.2 and 6.1 kcal/mol for H and D transfer, respectively, which are depicted as horizontal lines in Figure 7. It is noteworthy that the MTE is G



Article

We first consider the first step of the reaction, i.e., the formation of the reactant complex. Figure 9 displays the

Figure 8. Thermal distribution of the transmission probability, P(E) exp(−E/kBT), obtained from the small-curvature tunneling (SCT) and large-curvature tunneling (LCT) calculations at the MPW1K level. Zero of energy is set at the reactant complex.

Figure 9. Comparison of binding energy in the gas phase obtained with various electronic structure methods. Geometry optimization was performed at the individual levels of theory, except for the MRMP calculation where the optimized geometry at the MPW1K level was used.79

lower than the top of the barrier by 6−7 kcal/mol, indicating that the sharp barrier in the region −0.5 < s < 0.5 is essentially eliminated by deep tunneling. Table 3 summarizes the reaction rates obtained by including the transmission coefficient at various levels of approximation. By reflecting the extremely large values of κ, the rate constant with tunneling correction for H transfer is 3 orders of magnitude greater than the corresponding CVT rate constant. The KIE value increases up to 12 and 32 at the SCT and LCT levels, respectively (Table 5). The OMT scheme gives an

binding energy of the donor−acceptor complex calculated with various electronic structure theories. Here, the binding energy, EB, is defined as the energy difference between the separated reactants and the reactant complex. The molecular geometries were fully optimized at the individual levels of theory, except for the MRMP calculation.79 The MPW1K and B3LYP methods predict a similar binding energy of about 5 kcal/ mol, which is consistent with the previous finding for typical PCET reactions.45 More notably, all the other methods tested here predict a much greater binding energy for the present system. For example, the ROHF/MP2 method gives EB = 9.8 kcal/mol, which is about 5 kcal/mol greater than the B3LYP and MPW1K results. The difference is probably attributed to the inadequate description of the weak noncovalent interactions (such as dispersion) considering the fact that the B3LYP method is not capable of describing such interactions.80 To obtain more insight, we calculated the binding energy using dispersion and/or long-range corrected DFT methods81−83 (here, the BHHLYP-D3, LC-BLYP, and LC-BOP + LRD methods were used). These methods are found to give a rather similar value of EB in the range of 11−12 kcal/mol, which is close to the ROHF/MP2 result. We also tested the M06-2X method that is known to provide a well-balanced description of the reaction barrier and weak noncovalent interactions.84 The obtained result is EB = 9.9 kcal/mol, which is again close to the ROHF/MP2 result. Therefore, it turns out that the binding energy increases up to 10−12 kcal/mol by taking into account the weak noncovalent interactions. We have also compared the optimized geometries obtained with the respective methods (Table S1 of the Supporting Information). Unlike the comparison of binding energy, we find that the optimized geometry is not particularly sensitive to the choice of the electronic structure methods. For example, the variation of the softest degrees of freedom (i.e., the dihedral angles ϕ and ω) is at most ∼7°. This may be reasonable because the intra- and intermolecular H bonds play a more crucial role than weak noncovalent interactions in determining

Table 5. Deuterium Kinetic Isotope Effect (KIE) for the Reaction in Solution Obtained from the VTST/MT Calculation at the MPW1K Levela

a

T(K)

CVT

ZCT

SCT

LCT

OMT

expt

288 298 308

6.67 6.28 5.99

10.3 9.85 9.56

12.2 11.9 11.5

34.8 32.3 30.3

17.3 16.4 15.6

21.4 21.6 24.4

The experimental values are obtained from ref 37.

intermediate KIE value, 16.4, which compares most favorably with the experimental result, 21.6. Nevertheless, it should be noted that the absolute reaction rate differs substantially from the experimental result. For example, the rate constant at 298 K obtained at the OMT level (2.89 M−1 s−1) differs by 4 orders of magnitude from the experimental value (4.11 × 104 M−1s−1). The error may be attributed to an insufficient accuracy of the potential energy surface calculated at the MPW1K level, although the latter is known to predict accurate barrier heights for H-transfer reactions. In the next section we thus re-examine the accuracy of the MPW1K method and make an appropriate energy correction on the reaction rate. 3.4. Effect of Donor−Acceptor Binding on the Overall Reaction Barrier. In this section we re-examine the accuracy of the energy profile by using a variety of electronic structure methods. To do so, we separate the overall reaction into two steps, namely, (i) the formation of the reactant complex from the separated reactants and (ii) the activation of the reactant complex to the transition state. H



Article

the optimized geometry. The same trend can be observed for the optimized geometry of the transition state, i.e., the geometry is not significantly affected by weak noncovalent interactions except for the B3LYP method (Table S2 of the Supporting Information). We next consider the second step of the reaction, i.e., the activation of the reactant complex to the transition state. Figure 10 compares the energy profiles obtained with the individual

Figure 11. Energy diagram (in kcal/mol) for the reaction in the gas phase and in solution obtained with energy correction at the M06-2X level. Figures in parentheses include zero-point energies.

intrinsic barrier. Consequently, the overall barrier as measured from the separated reactants becomes ∼5 kcal/mol lower than at the MPW1K level. Because we are interested in the secondorder reaction rate that corresponds to the experimental condition,37 we assume that the reaction rate is governed by the free-energy barrier measured from the separated reactants rather than from the reactant complex. With this assumption we have calculated the modified VTST/MT rate with singlepoint energy correction at the M06-2X level (Table 6). The

Figure 10. Gas-phase potential energy profile in the TS region obtained with various electronic structure methods. Zero of energy is set at the reactant complex (RC). The minimum energy path (MEP) was recalculated at the individual levels of theory, except for the ROHF/MP2 and LC-BLYP calculations where the MEP at the MPW1K level was used.85

Table 6. Modified Rate Constant (M−1 s−1) for the Reaction in Solution Obtained with the Single-Point Energy Correction at the M06-2X Levela T(K)

methods.85 These profiles differ from each other mainly in the region −0.5 < s < 0.5 that corresponds to the local H migration between the donor and acceptor oxygens. We find that the ROHF/MP2 method gives the highest barrier of 21 kcal/mol (as measured from the reactant complex). This may be reasonable because the ROHF/MP2 method does not account for the diabatic coupling between the reactant and product states and thus exhibits a cusp at s ≃ 0. The B3LYP method gives a very low barrier height of 10 kcal/mol, as expected. This is consistent with the previous observation that the B3LYP method greatly underestimates the reaction barrier for typical PCET reactions.45 The MPW1K and M06-2X methods yield very similar energy profiles with an intermediate barrier height of 16−17 kcal/mol. To obtain further insight, we also performed mulireference wave function theory calculation of the energy profile (see the Supporting Information). While the latter calculation was met with some difficulty because of the limited size of active space, we find that it provides a rough estimate of the intrinsic barrier to be ∼16 kcal/mol. Therefore, we expect that the true intrinsic barrier will be around 16−18 kcal/mol based on the relatively good agreement between several different methods. On the basis of the above results, we have performed singlepoint energy correction on the reaction rate. For this purpose we have utilized the M06-2X method considering the wellbalanced description of the binding energy and the intrinsic barrier. Figure 11 displays the revised energy diagram in the gas phase and in solution obtained by the energy correction (see the Supporting Information for details). As seen, the energy profile at the M06-2X level has a binding well that is deeper than that at the MPW1K level while retaining the height of the

a

kCVT mod

288 298 308

5.41 8.51 13.0

288 298 308

0.81 1.36 2.17

kCVT/OMT mod H-transfer 3.48 × 3.53 × 3.60 × D-transfer 2.00 × 2.16 × 2.29 ×

expt

104 104 104

3.70 × 104 4.11 × 104 4.93 × 104

103 103 103

1.73 × 103 1.90 × 103 2.02 × 103

The experimental values are obtained from ref 37.

modified CVT rate thus obtained (denoted as kCVT mod ) is 4 orders of magnitude greater than that at the MPW1K level by reflecting the lower reaction barrier by ∼5 kcal/mol. To estimate the tunneling effect, we further assume that the transmission coefficient κ at the M06-2X level is approximately equal to that at the MPW1K level. This approximation is motivated by the fact that κ is determined essentially by the local shape of the potential barrier, which is very similar between the MPW1K and M06-2X methods (Figure 10). To test the validity of this approximation, we have recalculated the value of κ at the ZCT and SCT levels using the M06-2X method (see Table S3 of the Supporting Information). The obtained results are found to agree reasonably well with each other, supporting the use of the MPW1K values of κ in Table 4 (see the Supporting Information for more details). On the basis of this observation we have calculated the modified reaction rate as follows: CVT/OMT CVT k mod ≃ κOMTk mod

I

(14)



Article

where kCVT mod is the modified CVT rate at the M06-2X level and κOMT is the OMT transmission coefficient at the MPW1K level (Table 4). Table 6 displays the results thus obtained, which are now in fairly good agreement with experiment.37 This level of agreement may be partly fortuitous considering a variety of approximations involved. Nevertheless, the above result suggests that accounting for the weak noncovalent interactions (such as dispersion) can be crucial for an accurate evaluation of reaction rate because they may lower the overall reaction barrier by shifting the energy profile downward in the TS region. 3.5. Discussion. Relation between the Arrhenius Activation Energy and Mean Tunneling Energy. In the preceding sections we have seen that the regeneration reaction of vitamin E by ubiquinol is facilitated by two factors, namely, quantum tunneling and donor−acceptor binding. The tunneling effect increases the reaction rate by penetration of the potential barrier, while the donor−acceptor binding increases the rate by lowering the overall energy profile in the TS region. The situation is depicted schematically in Figure 12, where the

activation energy by using the above equation along with the calculated values of kCVT/OMT in Table 6 . The obtained value is mod Ea = 0.3 kcal/mol, which compares favorably with the experimental results (1−3 kcal/mol). It should be noted that those values of Ea are significantly smaller than typical barrier heights of hydrogen abstraction reactions (e.g., >10 kcal/mol). Experimentally, Ea has been interpreted as the energy of a quantum-mechanical “tunnel” through the barrier,39 as illustrated in Figure 12. This interpretation is appealing and can be rationalized as follows. First, we insert the VTST/MT rate expression (k(T) = κkCVT) into eq 15, which gives

Ea = EaCVT + ΔE(tun)

(16)

ECVT a

where is the activation energy corresponding to the CVT rate constant, namely EaCVT = −

d ln kCVT(T ) d(1/kBT )

(17)

and ΔE is the tunneling contribution (see below). The ECVT a represents the activation energy of a fictitious system with no tunneling, because kCVT(T) accounts only for the effect of zeropoint energies. The ECVT obtained from the data in Table 6 is a 7.8 kcal/mol, which is close to the adiabatic free-energy barrier, 6.3 kcal/mol (see Figure 11). This indicates that the following relation holds: (tun)

EaCVT ≃ ga⧧

(18)

where g⧧a ≡ ga(s⧧) is the height of the adiabatic free-energy barrier as measured from the separated reactants. As for the tunneling contribution ΔE(tun) in eq 16, one can readily show that the following relation holds (see the Supporting Information): ΔE(tun) ≃ ⟨E⟩tun − ga⧧

(19)

where ⟨E⟩tun is the mean tunneling energy (MTE) given by Figure 12. Schematic illustration of the adiabatic free-energy profile for systems (a) with no donor−acceptor binding (e.g., the toluene/ benzyl complex) and (b) with binding on both the reactant and product sides (e.g., the phenol/phenoxyl complex). Note that the donor−acceptor binding lowers the overall reaction barrier as measured from the separated reactants. Ea denotes the Arrhenius activation energy, which satisfies Ea ≃ g⧧a + ΔE(tun), where g⧧a is the height of the adiabatic free-energy barrier and ΔE(tun)(4000. We have attributed the large tunneling effect to the narrowness of the potential barrier associated with a compact H-bonded complex. The rather distorted form of the donor−acceptor complex avoids the steric hindrance of side chains while allowing the formation of a compact O···H···O interface. This is in contrast to the phenol/phenoxyl complex (i.e., a prototypical PCET system) which assumes a nearly coplanar optimized geometry. In terms of the electronic character, the present reaction may be viewed as a PCET reaction with a relatively large energy gap at the TS (14 kcal/mol) because of the distorted form of the donor−acceptor complex. The absolute reaction rates and the activation energy obtained from the VTST/MT calculation are found to be in good agreement with the experiments of Nagaoka and co-workers.37−39 A systematic comparison was then made among relevant antioxidant reactions, and it is

■

REFERENCES

(1) Sies, H. Angew. Chem., Int. Ed. Engl. 1986, 25, 1058−1071. (2) Bokov, A.; Chaudhuri, A.; Richardson, A. Mech. Ageing Dev. 2004, 125, 811−826. (3) Finkel, T.; Holbrook, N. J. Nature 2000, 408, 239−247. (4) Muller, F. L.; Lustgarten, M. S.; Jang, Y.; Richardson, A.; Remmen, H. V. Free Radical Biol. Med. 2007, 43, 477−503. (5) Chaudiere, J.; Ferrari-Iliou, R. Food Chem. Toxicol. 1999, 37, 949−962. (6) Girotti, A. W. J. Lipid Res. 1998, 39, 1529−1542. (7) Yin, H.; Xu, L.; Porter, N. A. Chem. Rev. 2011, 111, 5944−5972. (8) Burton, G. W.; Ingold, K. U. Acc. Chem. Res. 1986, 19, 194−201. (9) Burton, G. W.; Ingold, K. U. Ann. N.Y. Acad. Sci. 1989, 570, 7− 22. (10) Wang, X.; Quinn, P. J. Prog. Lipid Res. 1999, 38, 309−336. (11) Niki, E.; Noguchi, N. Acc. Chem. Res. 2004, 37, 45−51. (12) Atkinson, J.; Epand, R. F.; Epand, R. M. Free Radical Biol. Med. 2008, 44, 739−764. (13) Bowry, V. W.; Stocker, R. J. Am. Chem. Soc. 1993, 115, 6029− 6044. L



Article

(14) Wang, Y.; Hekimi, S. Crit. Rev. Biochem. Mol. Biol. 2013, 48, 69− 88. (15) Arroyo, A.; Navarro, F.; Navas, P.; Villalba, J. M. Protoplasma 1998, 205, 107−113. (16) Frei, B.; Kim, M. C.; Ames, B. N. Proc. Natl. Acad. Sci. U.S.A. 1990, 87, 4879−4883. (17) Kagan, V. E.; Quinn, P. J. Coenzyme Q: Molecular Mechanism in Health and Disease; CRC Press: Boca Raton, FL, 2001. (18) Ernster, L.; Dallner, G. Biochim. Biophys. Acta 1995, 1271, 195− 204. (19) Kagan, V. E.; Serbinova, E. A.; Packer, L. Biochem. Biophys. Res. Commun. 1990, 169, 851−857. (20) Kagan, V. E.; Serbinova, E. A.; Stoyanovsky, D. A.; Khwaja, S.; Packer, L. Methods Enzymol. 1994, 234, 343−354. (21) Mukai, K.; Kikuchi, S.; Urano, S. Biochim. Biophys. Acta 1990, 1035, 77−82. (22) Mukai, K.; Itoh, S.; Morimoto, H. J. Biol. Chem. 1992, 267, 22277−22281. (23) Ingold, K. U.; Bowry, V. W.; Stocker, R.; Walling, C. Proc. Natl. Acad. Sci. U.S.A. 1993, 90, 45−49. (24) Espinosa-Garcia, J.; Gutierrez-Merino, C. J. Phys. Chem. A 2003, 107, 9712−9723. (25) Espinosa-Garcia, J. J. Am. Chem. Soc. 2004, 126, 920−927. (26) Navarrete, M.; Rangel, C.; Corchado, J. C.; Espinosa-Garcia, J. J. Phys. Chem. A 2005, 109, 4777−4784. (27) Navarrete, M.; Rangel, C.; Espinosa-Garcia, J.; Corchado, J. C. J. Chem. Theory Comput. 2005, 1, 337−344. (28) Tejero, I.; Gonzalez-Garcia, N.; Gonzalez-Lafont, A.; Lluch, J. M. J. Am. Chem. Soc. 2007, 129, 5846−5854. (29) Chiodo, S. G.; Leopoldini, M.; Russo, N.; Toscano, M. Phys. Chem. Chem. Phys. 2010, 12, 7662−7670. (30) Markovic, Z.; Amic, D.; Milenkovic, D.; Dimitric-Markovic, J. M.; Markovic, S. Phys. Chem. Chem. Phys. 2013, 15, 7370−7378. (31) Di Meo, F.; Lemaur, V.; Cornil, J.; Lazzaroni, R.; Duroux, J.; Olivier, Y.; Trouillas, P. J. Phys. Chem. A 2013, 117, 2082−2092. (32) Iuga, C.; Alvarez-Idaboy, J. R.; Vivier-Bunge, A. J. Phys. Chem. B 2011, 115, 12234−12246. (33) Iuga, C.; Alvarez-Idaboy, J. R.; Russo, N. J. Org. Chem. 2012, 77, 3868−3877. (34) Leopoldini, M.; Chiodo, S. G.; Russo, N.; Toscano, M. J. Chem. Theory Comput. 2011, 7, 4218−4233. (35) Wright, J. S.; Johnson, E. R.; DiLabio, G. A. J. Am. Chem. Soc. 2001, 123, 1173−1183. (36) Medina, M. E.; Iuga, C.; Alvarez-Idaboy, J. R. Phys. Chem. Chem. Phys. 2013, 15, 13137−13146. (37) Nagaoka, S.; Nishioku, Y.; Mukai, K. Chem. Phys. Lett. 1998, 287, 70−74. (38) Nagaoka, S.; Inoue, M.; Nishioka, C.; Nishioku, Y.; Tsunoda, S.; Ohguchi, C.; Ohara, K.; Mukai, K.; Nagashima, U. J. Phys. Chem. B 2000, 104, 856−862. (39) Ouchi, A.; Nagaoka, S.; Mukai, K. J. Phys. Chem. B 2010, 114, 6601−6607. (40) Kakiuchi, T.; Mukai, K.; Ohara, K.; Nagaoka, S. Bull. Chem. Soc. Jpn. 2009, 82, 216−218. (41) Nagaoka, S.; Kakiuchi, T.; Ohara, K.; Mukai, K. Chem. Phys. Lipids 2007, 146, 26−32. (42) Kwart, H. Acc. Chem. Res. 1982, 15, 401−408. (43) Litwinienko, G.; Ingold, K. U. Acc. Chem. Res. 2007, 40, 222− 230. (44) Mayer, J. M.; Hrovat, D. A.; Thomas, J. L.; Borden, W. T. J. Am. Chem. Soc. 2002, 124, 11142−11147. (45) Lingwood, M.; Hammond, J. R.; Hrovat, D. A.; Mayer, J. M.; Borden, W. T. J. Chem. Theory Comput. 2006, 2, 740−745. (46) Hammes-Schiffer, S.; Soudackov, A. V. J. Phys. Chem. B 2008, 112, 14108−14123. (47) Hammes-Schiffer, S. J. Phys. Chem. Lett. 2011, 2, 1410−1416. (48) Skone, J. H.; Soudackov, A. V.; Hammes-Schiffer, S. J. Am. Chem. Soc. 2006, 128, 16655−16663.

(49) Tishchenko, O.; Truhlar, D. G.; Ceulemans, A.; Nguyen, M. T. J. Am. Chem. Soc. 2008, 130, 7000−7010. (50) Cembran, A.; Provorse, M. R.; Wang, C.; Wu, W.; Gao, J. J. Chem. Theory Comput. 2012, 8, 4347−4358. (51) Weinberg, D. R.; Gagliardi, C. J.; Hull, J. F.; Murphy, C. F.; Kent, C. A.; Westlake, B. C.; Paul, A.; Ess, D. H.; McCafferty, D. G.; Meyer, T. J. Chem. Rev. 2012, 112, 4016−4093. (52) Lynch, B. J.; Fast, P. L.; Harris, M.; Truhlar, D. G. J. Phys. Chem. A 2000, 104, 4811−4815. (53) Gonzalez, C.; Schlegel, H. B. J. Chem. Phys. 1989, 90, 2154− 2161. (54) Nakano, H.; Uchiyama, R.; Hirao, K. J. Comput. Chem. 2002, 23, 1166−1175. (55) Ebisuzaki, R.; Watanabe, Y.; Nakano, H. Chem. Phys. Lett. 2007, 442, 164−169. (56) Sato, H.; Kovalenko, A.; Hirata, F. J. Chem. Phys. 2000, 112, 9463−9468. (57) Hirata, F. Understanding Chemical Reactivity: Molecular Theory of Solvation; Kluwer: Dordrecht, 2003. (58) Kovalenko, A.; Hirata, F. J. Chem. Phys. 1999, 110, 10095− 10112. (59) Jorgensen, W. L. J. Phys. Chem. 1986, 90, 1276−1284. (60) Truhlar, D. G.; Isaacson, A. D.; Garrett, B. C. Theory of Chemical Reaction Dynamics; CRC Press: Boca Raton, FL, 1985; Vol. 4; pp 65− 137. (61) Lu, D.; Truong, T. N.; Melissas, V. S.; Lynch, G. C.; Liu, Y.-P.; Garrett, B. C.; Steckler, R.; Isaacson, A. D.; Rai, S. N.; et al. Comput. Phys. Commun. 1992, 71, 235−262. (62) Kuppermann, A.; Truhlar, D. G. J. Am. Chem. Soc. 1971, 93, 1840−1851. (63) Garrett, B. C.; Truhlar, D. G.; Grev, R. S.; Magnuson, A. W. J. Phys. Chem. 1980, 84, 1730−1748. (64) Liu, Y. P.; Lynch, G. C.; Truong, T. N.; Lu, D. H.; Truhlar, D. G.; Garrett, B. C. J. Am. Chem. Soc. 1993, 115, 2408−2415. (65) Garrett, B. C.; Joseph, T.; Truong, T. N.; Truhlar, D. G. Chem. Phys. 1989, 136, 271−293. (66) Fernandez-Ramos, A.; Truhlar, D. G. J. Chem. Phys. 2001, 114, 1491−1496. (67) Liu, Y. P.; Lu, D. H.; Gonzalez-Lafont, A.; Truhlar, D. G.; Garrett, B. C. J. Am. Chem. Soc. 1993, 115, 7806−7817. (68) Fernandez-Ramos, A.; Truhlar, D. G. J. Chem. Theory Comput. 2005, 1, 1063−1078. (69) Truhlar, D. G.; Liu, Y.-P.; Schenter, G. K.; Garrett, B. C. J. Phys. Chem. 1994, 98, 8396−8405. (70) Schmidt, M. W.; Baldridge, K. K.; Boatz, J. A.; Elbert, S. T.; Gordon, M. S.; Jensen, J. H.; Koseki, S.; Matsunaga, N.; Nguyen, K. A.; Su, S.; et al. J. Comput. Chem. 1993, 14, 1347−1363. (71) The code for the MPW1K method was obtained from the GAMESSPLUS package (ref S8 of the Supporting Information). (72) de Heer, M. I.; Korth, H.; Mulder, P. J. Org. Chem. 1999, 64, 6969−6975. (73) Lucarini, M.; Pedrielli, P.; Pedulli, G. F.; Cabiddu, S.; Fattuoni, C. J. Org. Chem. 1996, 61, 9259−9263. (74) Nonadiabatic extensions of transition-state theory, e.g., based on the Landau−Zener or Zhu−Nakamura formula for nonadiabatic transition, may be useful for this purpose. (75) More precisely, the minimum energy path (MEP) on the freeenergy surface G(R) in solution is referred to as the minimum freeenergy path. (76) Pang, J.; Hay, S.; Scrutton, N. S.; Sutcliffe, M. J. J. Am. Chem. Soc. 2008, 130, 7092−7097. (77) We note that the “tunneling-ready” configuration is not necessarily a stationary point on the free-energy surface. In the present case, the tunneling-ready configuration corresponds roughly to s ≃ −0.5, at which the dominant component of the reaction coordinate changes qualitatively from the slow donor−acceptor motion to the fast hydrogenic motion. A similar type of configuration is discussed in the context of electron transfer reactions (see ref 78, for example). M



Article

(78) Kreevoy, M. M.; Ostovic, D.; Truhlar, D. G.; Garrett, B. C. J. Phys. Chem. 1986, 90, 3766−3744. (79) In the MRMP calculation of binding energy, we used the optimized geometry at the MPW1K level because analytical gradients for the MRMP method were not available. We expect that the present estimate of binding energy at the MRMP level will increase somewhat if the geometry is fully optimized at the MRMP level. (80) Grimme, S. Wiley Interdiscip. Rev.: Comput. Mol. Sci. 2011, 1, 211−228. (81) Grimme, S.; Antony, J.; Ehrlich, S.; Krieg, H. J. Chem. Phys. 2010, 132, 154104/1−19. (82) Iikura, H.; Tsuneda, T.; Yanai, T.; Hirao, K. J. Chem. Phys. 2001, 115, 3540−3544. (83) Sato, T.; Nakai, H. J. Chem. Phys. 2009, 131, 224104−224115. (84) Zhao, Y.; Truhlar, D. G. Acc. Chem. Res. 2008, 41, 157−167. (85) The reason for this is that the ROHF/MP2 and LC-BLYP methods exhibit a derivative discontinuity (or a cusp) at s ≃ 0. (86) More specifically, the experimental value of Ea for the regeneration reaction of 5,7- diisopropyl tocopherol by ubiquinol in ethanol and in hexane is 2.5 and 2.7 kcal/mol, respectively (refs 37 and 38). The Ea for the regeneration reaction of α-tocopherol by ubiquinol in ethanol is 1.1 kcal/mol (ref 39.). (87) Tolman, R. C. J. Am. Chem. Soc. 1920, 42, 2506−2528. (88) Truhlar, D. G. J. Chem. Educ. 1978, 55, 309−311. (89) Masgrau, L.; Gonzalez-Lafont, A.; Lluch, J. M. J. Comput. Chem. 1999, 20, 1685−1692. (90) The local barrier associated with the H migration is not identical to the intrinsic barrier (i.e., the energy difference between the TS and the reactant complex), which is not directly relevant to the degree of quantum tunneling. This is because the intrinsic barrier includes the energy for rearranging the donor−acceptor moieties to the precursor (or tunneling-ready) state. In reality, tunneling is governed by the local potential barrier from the precursor state. (91) Reichardt, C. Solvents and Solvent Effects in Organic Chemistry, 3rd ed.; Wiley-VCH: Weinheim, Germany, 2002. (92) In the case of weak tunneling, this can be understood analytically from the Wigner approximation to the transmission coefficient, i.e., κ ≃ 1 + c|ω⧧|2, where c is a constant and ω⧧ is the imaginary frequency. In the case of deep tunneling, the same conclusion can be obtained from path-integral quantum transitionstate theory. That is, the sharper the potential barrier, the lower than the classical potential barrier the centroid potential barrier becomes. (93) Zhou, F.; Schulten, K. J. Phys. Chem. 1995, 99, 2194−2207. (94) Marquardt, D.; Williams, J. A.; Kucerka, N.; Atkinson, J.; Wassal, S. R.; Katsaras, J.; Harroun, T. A. J. Am. Chem. Soc. 2013, 135, 7523− 7533. (95) Fiorini, R.; Ragni, L.; Ambrosi, S.; Littarru, G. P.; Gratton, E.; Hazlett, T. Photochem. Photobiol. 2008, 84, 209−214.

N


Critical Role of Deep Hydrogen Tunneling to Accelerate the

Recommend Documents