Tunneling Kinetics and Nonadiabatic Proton-Coupled Electron

Article pubs.acs.org/JPCB

Tunneling Kinetics and Nonadiabatic Proton-Coupled Electron Transfer in Proteins: The Effect of Electric Fields and Anharmonic Donor−Acceptor Interactions Bridget Salna, Abdelkrim Benabbas, Douglas Russo, and Paul M. Champion* Department of Physics and Center for Interdisciplinary Research on Complex Systems, Northeastern University, Boston, Massachusetts 02115, United States S Supporting Information *

ABSTRACT: A proper description of proton donor−acceptor (D−A) distance fluctuations is crucial for understanding tunneling in proton-coupled electron transport (PCET). The typical harmonic approximation for the D−A potential results in a Gaussian probability distribution, which does not appropriately reflect the electronic repulsion forces that increase the energetic cost of sampling shorter D−A distances. Because these shorter distances are the primary channel for thermally activated tunneling, the analysis of tunneling kinetics depends sensitively on the inherently anharmonic nature of the D−A interaction. Thus, we have used quantum chemical calculations to account for the D−A interaction and developed an improved model for the analysis of experimental tunneling kinetics. Strong internal electric fields are also considered and found to contribute significantly to the compressive forces when the D−A distance distribution is positioned below the van der Waals contact distance. This model is applied to recent experiments on the wild type (WT) and a double mutant (DM) of soybean lipoxygenase-1 (SLO). The compressive force necessary to prepare the tunneling-active distribution in WT SLO is found to fall in the ∼ nN range, which greatly exceeds the measured values of molecular motor and protein unfolding forces. This indicates that ∼60−100 MV/cm electric fields, aligned along the D−A bond axis, must be generated by an enzyme conformational interconversion that facilitates the PCET tunneling reaction. Based on the absolute value of the measured tunneling rate, and using previously calculated values of the electronic matrix element, the population of this tunneling-active conformation is found to lie in the range 10−5−10−7, indicating this is a rare structural fluctuation that falls well below the detection threshold of recent ENDOR experiments. Additional analysis of the DM tunneling kinetics leads to a proposal that a disordered (high entropy) conformation could be tunneling-active due to its broad range of sampled D−A distances.

■

INTRODUCTION Coupled proton and electron transfer occur in numerous chemical reactions,1−6 many of which have critical biological functions, such as cellular respiration7,8 and photosynthesis.9 For example, cytochromes P450, which utilize a thiolate-ligated heme and a proton transport system10−12 to cleave molecular oxygen, also mediate a key proton-coupled electron transfer step during the R−H to R−OH rebound replacement reaction.13−15 A variety of studies have been concerned with the role of proton tunneling14,16−21 in these reactions. Such studies involve experimental measurement1,3,22−24 of the kinetics and the kinetic isotope effect (KIE), as well as computational18,25,26 and analytical27−29 modeling of the tunneling rates. For example, the proton coupled electron transfer (PCET) occurring in wild-type (WT) soybean lipoxygenase-1 (SLO)1 exhibits a large room temperature KIE (∼81), which indicates that proton tunneling is involved. Recently, a double mutant (DM) of SLO was also studied22 and found to have an extremely large room temperature KIE © XXXX American Chemical Society

(∼500−700). Attempts to fully explain this observation, as well as other experimental measurements on the SLO system, are still ongoing.22,30,31 In this study, we present an analysis that focuses on the anharmonic donor−acceptor interaction as well as on the magnitude and nature of the protein-induced forces that facilitate the PCET reaction in SLO. One important conclusion is that strong protein electric fields32,33 are needed in order to generate the large forces necessary to initiate the enzymatic tunneling reaction. Theoretical calculations of the proton tunneling rate27,28,34,35 typically focus on the significant effect of the proton donor (D) and acceptor (A) atom distance fluctuations. In nearly all prior studies,27,28,34 the D−A distance is assumed to be modulated within a harmonic approximation, resulting in a Gaussian D−A distance distribution. In contrast, the tunneling coordinate (often a C-H or C-D mode in enzyme reactions) has been Received: June 6, 2017 Published: June 19, 2017 A

DOI: 10.1021/acs.jpcb.7b05570 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B treated using a Morse potential.22,36,37 In one instance25 the D−A coordinate was modeled using anharmonic terms derived from gas-phase valence bond matrix elements and a computational-based multistate continuum theory. As noted previously,36 such an approach requires many parameters to describe the empirical valence bond (EVB) Hamiltonian and this makes the resulting anharmonic effects difficult to interpret systematically in comparison to a simple analytic approach that is more useful to experimentalists. The primary goal of the present work is to explore in more detail the forces that govern the D−A coordinate and to examine how this might affect the analysis of the observed tunneling kinetics. The relative D−A motion dynamically changes the proton tunneling distance and this alters the wave function overlap between the reactant and product states. A Gaussian D−A distance distribution does not properly account for the repulsive electronic interactions that take place as the distance between D and A is reduced to smaller values during the compression stage of a D−A vibrational cycle or during a compressive protein fluctuation. This repulsive interaction is important because the compressive region explored by the D− A coordinate is where most of the tunneling amplitude is found. In order to better understand the effect of the anharmonic D−A electronic repulsion at short-range, we have used a quantum chemical calculation to construct a D−A interaction potential that appropriately reduces the probability amplitude at smaller D−A distances. As might be expected, this also reduces the calculated tunneling rate and minimizes the need for proton-level adiabatic corrections.38,39 Under some conditions where strong external forces are applied, the inclusion of the anharmonic electronic repulsion can tightly localize the D−A distribution. Such localization helps to improve the linear Taylor’s expansion approximation (“linear approximation”) that is sometimes used to treat the “coupling” of the D−A motion to the tunneling reaction.22,38,40−42 In this study we develop and apply a simple anharmonic D− A interaction potential that allows an improved physical picture to be extracted from the measured tunneling kinetics. This potential models the enzyme’s ability (through active site electronic effects, local electric fields, and conformational forces) to carry out a PCET reaction using the tunneling mechanism. As an example of an application to enzyme systems, we explore the interesting tunneling kinetics measured for WT soybean lipoxygenase-1 (SLO) and its double mutant (DM).1,22 This analysis offers new insights derived from both the tunneling kinetic data and from recent ENDOR studies of the Mn-substituted SLO system.31 More specifically, we find for the WT enzyme that ∼ nN forces are needed to position the D−A distance distribution within a range that is compatible with the observed tunneling kinetics. Such large forces are more than can be realized by protein conformational forces (i.e., nonelectrostatic forces) alone. However, the inclusion of a strong protein electric field, which acts on the D−A subsystem, provides a mechanism33 that allows the protein to reduce the D−A distance below the van der Waals contact threshold so that efficient tunneling can occur.

studied electronically nonadiabatic limit.27,28,43 The most commonly used analytical models27,28,34 factor out and treat the proton (or deuteron) and the donor−acceptor (vibrational) coordinates independently from the generalized “environmental” coordinate, which accounts (classically) for the remaining vibrational degrees of freedom of the enzymesolvent system. Thermal fluctuations along the environmental coordinate bring the reactant and product electronic states into degeneracy so that tunneling can take place.34,41 The high frequency (∼100 THz or 10 fs) tunnel particle vibration must be treated quantum mechanically,44 while the slower heavy atom (D-A) motion is usually treated within the Born− Oppenheimer approximation. These approximations have been tested previously in the harmonic limit using a correlation function method.38,45−47 When the D−A motion is treated in the classical limit, it can be considered as a “gating” coordinate so that the reaction rate depends upon an integral of the proton (or deuteron) wave function overlap between reactant and product states (i.e., the Franck−Condon factor), which is weighted by a classical temperature-dependent distribution of D−A tunneling distances.1,27,46,48,49 The Franck−Condon factor depends on the distance, S , that the proton or deuteron tunnels and this can be well-approximated using the distance between the donor and acceptor atoms, R, minus the covalent bonds (S = R − 2 lc̅ ), where lc̅ is the average of the donor and acceptor bond lengths. Parameters describing the D−A interaction potential energy control the temperature dependent distribution of D−A distances, PDA(R). When the kinetic data are taken near room temperature, as is often done for enzyme systems,1,22,23,36,50−52 the distribution can be constructed using a classical statistical limit approach.27,38,46,53 However, when a broad temperature range is studied, it is necessary to treat the D−A vibrational motion at the quantum level.38 When donor and acceptor oxygen atoms form a OH---O hydrogen bond in a protein environment, the position, width, and shape of the D−A distance distribution depends on both the “bare” hydrogen bonding potential as well as on the forces from the surrounding protein. Enzyme reactions, on the other hand, often involve an aliphatic carbon donor (substrate) and an active site oxygen atom as the acceptor, where an actual hydrogen bond may not always be formed. However, negative charge (q < 0) build-up on the oxygen atom at the active site of a metallo-enzyme can polarize the CH bond of a substrate so that CH---Oq bound states are formed. The D−A interactions can have differing strength, depending upon the enzymeinduced structural, chemical, and electric field modifications at the active site. When compared to more rigid protein systems,39 enzymes also have the ability to undergo much larger conformational changes that are associated with substrate binding and product release.16,51 Thus, it is possible that the average distance between D and A is too large, or |q| is too small, so that a stable D−A bound state is not always formed. In such cases, uncorrelated fluctuations (rather than vibrations within a D− A binding potential) determine the D−A distance distribution. This situation should be contrasted with a strong hydrogen bond, such as OH---O, within a rigid protein interior where conformational disorder can sometimes be minimized. In such a case, the local D−A vibrational motion is the dominant factor that modulates the tunneling distance.39

■

GENERAL CONSIDERATIONS In many enzyme systems, proton and electron transport occur simultaneously and involve different acceptor sites. This is denoted28,30 as concerted proton coupled electron transport (PCET) or electron−proton transfer (EPT) and usually allows the tunneling reaction rate to be evaluated within the wellB


Article

The Journal of Physical Chemistry B

■

ANHARMONIC POTENTIAL FOR THE LOCAL D−A INTERACTION In most studies that use analytical models to describe tunneling, the D−A motion is implicitly taken to be harmonic and the distribution of tunnel lengths in the statistical limit is assumed to be a Gaussian function.1,27,34,37,40,41 However, when the D− A equilibrium position is less than ∼3 Å and the temperature is increased to ∼298 K, the harmonic approximation will begin to allow overly deep penetration by the tail of the D−A wave function (or its classical probability distribution) into regions where the tunnel distance becomes artificially small. It is precisely this region where most of the tunneling amplitude for the transport rate is generated. Thus, in order to move beyond the simple D−A harmonic approximation and improve analytical models for proton tunneling kinetics, the D−A distribution should include the electronic repulsion between the donor and acceptor atoms that occurs at smaller distances. This can be accomplished by using a more realistic potential for the local, or “bare”, atomic DH---A interaction. In addition, the surrounding protein can “dress” the “bare” potential and exert forces associated with dynamical conformational changes or with internal electric fields. Each of these factors has the potential to play an important role in determining the overall D−A distribution. As discussed in the Supporting Information (SI) section S1, we have calculated a variety of “bare” DH---A binding potentials, Vbare(R), at the full quantum level by utilizing density functional theory (DFT). A hybrid density functional (APF-D54) was used, which includes empirical dispersion terms and has been shown to perform well for systems where hydrogen bonds play an important role. The high accuracy basis set, aug-cc-pVTZ,55,56 was also used. In most cases we used relatively small “test” donors and acceptors because the primary goal was to approximate the repulsive part of the atomic interaction potential. In order to describe OH---O hydrogen bonding, we used water as the acceptor and either water or the hydroxyl group of serine as a donor. Methane was used as a test donor to mimic a simple aliphatic CH bond interacting with an enzyme active site oxygen atom acceptor in a H3CH---[OH]δ− “bonding” arrangement (denoted as CH--[OH]δ− throughout the text). To model the acceptor, the proton attached to the hydroxyl oxygen atom carries an adjustable nuclear charge that controls the net value of δ. As shown in the SI (section S1), alternative methods for acceptor charge variation confirm that this approach captures the basic electronic effects of negative charge build-up on the oxygen atom. The overall negative charge associated with the test acceptor polarizes the donor CH bond so that attraction and “binding” can occur. As discussed in the SI (section S2), the full quantum chemical calculation of a more realistic enzyme active site model complex determines that δ ∼ 0.6 for SLO, which is based on the average charge, q∼ − 0.9e, that is localized on the oxygen atom. It is noteworthy that as δ → 1, the test acceptor evolves into a hydroxyl anion and the CH---[OH]1.0− binding energy scale becomes similar to that of a OH---O hydrogen bond. However, the potential energy minimum occurs at ∼3.0 Å for the CH---[OH]1.0− interaction, rather than ∼2.7−2.8 Å as found for the OH---O hydrogen bond. Once Vbare(R) is selected to mimic the oxygen atom acceptor associated with the enzyme active site, it can be “dressed’ with additional potentials that mimic the effect of the protein

surroundings. This approach saves enormous computation time and yields a simple set of parameters that can be used directly to interpret experiments. As shown in the SI section S1, the analysis of the test potentials demonstrates that a Morse function can be used as an effective analytical approximation to the “bare” D−A potential (see section S1 and Table S1.1 of the SI). An example potential for δ = 0.6 along with its Morse fit is shown in Figure 1 where the D−A equilibrium distance, RDA 0 , is

Figure 1. DFT-generated D−A interaction potential, Vbare(R) is shown for δ = 0.6 fit with a Morse potential. The inset shows the effect of an applied electric field, ER, on the equilibrium D−A distance; the red points are calculated using quantum chemical methods and the analytic expression is described in section 5.2 of the SI. The field component (ER) along R, the D−A bond axis, points from the carbon donor (D) to the oxygen acceptor (A).

to found at ∼3.3 Å. Generally, increasing δ will shift RDA 0 smaller values along with the repulsive part of the potential (see Figure S1.2 of SI), which is the most important region to characterize when calculating the tunneling rate. The Morse parameters derived from a variety of hydrogen bonding test calculations are given in Table S1.1 of the SI, so that different donor and acceptor hydrogen bonding combinations can be compared. Once the bare D−A interaction potential, Vbare(R), has been determined, a potential, V⟨i⟩ conf(R), due to a specific global protein conformation, ⟨i⟩, can be added in order to include the effect of nonelectrostatic conformational forces from the surrounding protein. Such forces can either extend or compress the D−A subsystem with respect to the “bare” equilibrium, RDA 0 , favored by the isolated atomic forces. In order to fully account for other important active site forces, it is also advisible to explicitly consider the effect of protein electric fields, since they have been suggested to play an important role in enzyme catalysis.32,33 For an enzyme in global conformation ⟨i⟩, a simple analytical electrostatic potential, V⟨i⟩ field(R), can be added to Vbare(R) to test the effect of an applied electric field (see SI section S5.2). The analytical results compare very well to a full quantum chemical calculation as can be seen from the inset of Figure 1 where the D−A equilibrium distance as a function of applied field is plotted using both methods (see eq 5.3 in SI). When the component of the electric field along the Rcoordinate axis, ER, is oriented to induce stronger binding in the CH---[OH]δ− complex, it naturally leads to a reduction of the D−A equilibrium distance as shown in Figure 1. It is apparent that, for a given global protein conformation, ⟨i⟩, with charge, C


Article

The Journal of Physical Chemistry B q⟨i⟩, on the oxygen acceptor, the protein electric field acts with a ⟨i⟩ ⟨i⟩ force: F⟨i⟩ field ≌q ER (see section S53 of SI).

that the tail of the D−A distribution can penetrate to unrealistically short distances. In addition to the conformational force, it is possible that a strong electric field, E⟨i⟩ R , exists along the R-coordinate (from ⟨i⟩ ⟨i⟩ CH to O) and introduces a force, F⟨i⟩ field = q ER . This force can be combined with the conformational force to describe the overall protein compression force:

■

PROTEIN CONFORMATIONAL FORCES AND THE EFFECT OF ELECTRIC FIELDS Turning back to the effect of different protein conformations, we consider a protein global conformation ⟨i⟩ that is approximated by a simple harmonic conformational potential, ⟨i⟩ Vconf (R). In the absence of Vbare(R), the conformational potential would separate the D and A atoms by an average ⟨i⟩ DA equilibrium distance denoted as R⟨i⟩ 0 . When R0 > R0 , the protein conformational potential keeps D and A apart by introducing an energy penalty for moving the donor and acceptor atoms into a range where the local atomic D−A DA binding is optimized. On the other hand, if R⟨i⟩ 0 < R0 , the protein conformational force compresses the D−A equilibrium distance to a position that is shorter than RDA 0 . Thus, the addition of the protein conformational potential has the effect of embedding the bare atomic forces, associated with the isolated D−A subsystem, into the protein conformational DA landscape. In the limit that R⟨i⟩ 0 is significantly larger than R0 , binding by the local atomic potential is prevented and uncorrelated fluctuations of the donor and acceptor atoms within the conformational potential play a major role in determining the D−A distribution. The conformational restoring force constant is denoted by f⟨i⟩ conf as discussed in section S3.2 of the SI. As noted above (see Figure 1), the anharmonic nature of the D−A interaction potential (and its important electronic repulsion at short distance) is adequately approximated by a Morse function. Thus, using a harmonic potential to mimic the conformational forces for a given global protein conformation, ⟨i⟩, along with its electric field potential, V⟨i⟩ field(R), we can express the overall D−A potential as ⟨i⟩ ⟨i⟩ ⟨i⟩ ⟨i⟩ V DA (R ) = V bare (R ) + V conf (R ) + V field (R )

⟨i⟩ ⟨i⟩ ⟨i⟩ Fcompress = Ffield + Fconf

where the potential V⟨i⟩ conf(R) is now assumed to be nearly linear across the sub-Å width of a typical D−A vibrational distribution (e.g., see Figure 3b,c). Under this circumstance, its derivative, F⟨i⟩ conf, is nearly a constant that is limited by the maximum nonelectrostatic (i.e., conformational) force that can be exerted by a protein. We gauge this limit based on the forces measured for conformational unfolding (100−300 pN)57−59 or molecular motor activity (≲ 50 pN).60,61 In the event that ∼ nN forces are needed to produce a tunneling-active D−A distribution centered below the van der Waals contact distance, the dominant force must be supplied by F⟨i⟩ field. For a ∼ nN force, this requires protein electric fields on the order of ∼60−100 MV/cm, depending on the magnitude of q⟨i⟩. It is of interest that such strong fields have recently been documented for enzymatic systems.32,62,63 Importantly, when eq 2 is integrated, it can be used to “dress” Vbare(R) using only a single parameter, ⟨i⟩ F compress , to determine the D−A distribution (the Rindependent constant of integration can be dropped). When R⟨i⟩ 0 locates at an equilibrium position that differs from DA R0 in eq 1, or if protein electric field forces are present, the donor and acceptor atoms will respond to all interaction forces in eq 1a. It is apparent that this interplay can affect the overall shape of the D−A distribution as well as its important temperature dependence. Some examples of how the protein conformational forces might affect the OH---O hydrogen bonding interaction are given in section S4 of the SI. In the classical limit, the D−A distribution for global protein conformation, ⟨i⟩, is found using

(1a)

Setting the electric field to zero for the moment, we have ⟨i⟩ V DA (R ) = DDA [1 − exp(−bDA (R − R 0DA ))]2 1 ⟨i⟩ + f conf (R − R 0⟨i⟩)2 2

(2)

⎧ −V ⟨i⟩ (R ) ⎫ ⟨i⟩ DA ⎬ PDA (R ) = N exp⎨ ⎩ kBT ⎭ (1b)

⎪

⎪

⎪

⎪

(3a)

where N normalizes the distribution so that it integrates to unity over the domain of R at temperature T. In eq 3a, we are focusing on a level of the protein hierarchy that is associated with discrete global protein conformations, ⟨i⟩ (e.g., “open” and “closed” or “active” and “inactive”), where the associated V⟨i⟩ DA(R) might be very different (e.g., large differences in average D−A equilibrium distances, force constants, or electric fields). In the ergodic limit, the relative populations of these differing global conformations are determined by their equilibrium free energies so that the overall distribution of D−A distances found within the ensemble can be written as

where R is the distance between D and A. The parameters (bDA, DDA, and RDA 0 ) apply to a particular bare DH---A atomic interaction, as found in Table S1.1 (see also Figure S1.2 of the SI for more details). For simplicity in eq 1b, the index ⟨i⟩ is dropped for the bare atomic interaction potential. However, when there are major changes in moving from one global protein conformation to another (e.g., ⟨i⟩ → ⟨j⟩), changes in V⟨i⟩ bare(R) can be considered. For example, we tested this possibility for SLO by increasing δ for the tunneling active WT state (vide infra). Finally, we note that terms in the overall potential, which are independent of R, appear as energy offsets and are dropped because they play no role in determining the distribution of D−A distances. After a quantum chemical calculation has established δ for a specific enzyme active site, there are two unknown parameters ⟨i⟩ (R⟨i⟩ 0 and fconf) that determine the D−A probability distribution ⟨i⟩ if Vconf(R) is used to dress the bare potential. The more commonly used harmonic D−A potential (with a Gaussian D− A distance distribution) also invokes two parameters, but it does not account for electronic repulsion at small separation so

PDA(R ) =

⟨i⟩ (R ) ∑ p⟨i⟩ PDA i

(3b)

where p⟨i⟩ denotes the population of a global protein conformational state, ⟨i⟩, which is composed of a subset of associated conformational substates, i. As discussed in more detail in section S3.2 of the SI, this simple method for treating the hierarchical nature of protein conformational states allows us to consider an average tunneling rate for the global state ⟨i⟩ so long as the time scale for its average rate is slower than the D


Article


0⟩R|2, is discussed elsewhere42 (and in section S3 of SI) and uses symmetric back-to-back Morse potentials that are parametrized by an average of the donor and acceptor covalent bonds. In the following section, we explore the anharmonic P⟨i⟩ DA(R) distributions of SLO that are found to be consistent with the tunneling kinetics of the WT enzyme as well as with the interesting double mutant (DM), which has an unusually large KIE.22,30 A recent ENDOR analysis31 is also considered, so that the model developed is consistent with a wide range of experimental results.

time scale for conformational interconversions among the substates composing ⟨i⟩. This condition very likely holds in the case of SLO1,64 where the overall reaction rate, which is a product of p⟨i*⟩ and its associated PCET rate constant, is in the (ms)−1 regime for WT and much slower for the mutant protein.1,22 Clearly, if the average rate associated with the global conformation ⟨i⟩ increases and begins to approach or exceed the local substate fluctuation rates, the ergodic assumption breaks down and individual substates react independently, resulting in nonexponential kinetics.65−67

■

■

RATE EXPRESSIONS FOR NONADIABATIC PCET TUNNELING REACTIONS IN ENZYMES When the D−A motion is treated classically, the nonadiabatic tunneling rate can be found using the Golden Rule and the expressions for the Franck−Condon overlap associated with the tunnel particle zero-point vibrational states.28,34,38,42 If we assume for simplicity that only a single global conformation, ⟨i*⟩, is active in the tunneling process, the rate can be written as 2π |VRP|2 p⟨i*⟩ ⟨|⟨0 |0 ⟩|2 ⟩ DL AL R −ΔGs‡ / kBT ⟨k ⟩R = e ℏ 4πkBTλs

SOYBEAN LIPOXYGENASE TUNNELING KINETICS Experimental Tunneling Kinetics of SLO. The tunneling kinetics of the enzyme soybean lipoxygenase-1 have been studied extensively.1,22,30,36,76 In addition, structural measurements utilizing electron−nuclear double resonance (ENDOR)31 have provided independent constraints on the distances between the carbon donor of the linoleic acid substrate and the active site metal atom, as well as on the relative widths of the WT and DM distance distributions. Both the WT and the DM (with amino acid replacements L546A and L754A) were examined31 and are the focus of this analysis. The four key tunneling kinetic observables1,22,76 for these SLO species are listed in Table 1. A single mutant, I553G, was also examined with ENDOR31 but is not considered here (see section S5.4 of SI).

L

(4a)

with ⟨|⟨0DL|0AL⟩|2 ⟩R =

⟨i*⟩ (R )dR ∫ |⟨0|0⟩R |2 PDA

(4b)

where, for convenience on the right-hand side of eq 4b, we have dropped the subscripts for donor-state (D), acceptor-state (A), and tunnel particle (L = H or D). In eq 4a, the single tunnelingactive global conformation68,69 is denoted by the probability ΔH⟨i*⟩ ⎤ ⎡ ΔS⟨i*⟩ 1 factor p⟨i*⟩ ≈ 5 exp⎣⎢ k − k T ⎦⎥, as discussed in section S3.1 B B of the SI. We note that, if multiple global conformations are considered to be tunneling-active, the overall rate will involve their weighted average. When the free energy, ΔG0, of the PCET reaction is primarily enthalpic (as indicated for the SLO system1,25), eq 4a can be written using the approximation70−72 ΔG‡s ≈ ΔH‡s = (ΔH0 + λs)2/4λs, so that the temperature-independent contribution to the effective Arrhenius slope in eq 4a is given by ΔH‡s + ΔH⟨i*⟩. In addition, there is a sizable temperaturedependent contribution to the slope that arises from P⟨iDA*⟩(R) in eqs 3a and 4b. Here, we also note that when entropy is explicitly included in the reaction free energy, the temperature dependence of ΔG‡s becomes more complex.70 Thus, following others,71−73 we neglect such effects in this analysis. On the other hand, the entropic term associated with p⟨i*⟩ enters eq 4 as a “prefactor” and helps to account for the statistical probability of finding the tunneling-active conformation. Frictional effects74 can also be incorporated into eq 4, but they are neglected here as elsewhere.27,28,43 section S3 of the SI contains more details concerning eq 4. Finally, under many circumstances it is sufficient28,41 to consider only the vibrational overlap of the zero-point levels of the tunneling proton or deuteron; however, the effect of excited levels has also been studied.75 Inclusion of only the zero-point overlap generally applies41 when |ΔG0|≤ λs + ℏωL, where ℏωL ≫kBT is the quantum of vibrational energy for the tunneling particle of mass mL. However, as discussed previously,53 the higher vibrational states must be considered under certain conditions. The evaluation of the Franck−Condon factor, |⟨0|

Table 1. Experimental Observables for Wild Type SLO and Its Double Mutant SLO species

KIE (30 °C)

EHa (kcal/mol)

ΔEaa (kcal/mol)

kHcat (s−1) (30 °C)

ref.

WT DM

81 ± 5 729 ± 26,b

2.1 ± 0.2 9.9 ± 0.2

0.9 ± 0.2 n/ac

297 ± 12 0.021 ± 0.001

1 22

ΔEa = EDa − EHa . bAt 35 °C the KIE is 537 ± 55. This was determined by presteady-state-kinetics in ref 22. cThe value of ΔEa is not presented in ref 22 but, based on the similar D−A frequencies extracted from the analysis, it is estimated to be similar to the corresponding value in WT. a

In the WT form of the enzyme, the PCET reaction from the aliphatic carbon substrate to the ferric iron-hydroxide catalytic center reveals1 a KIE of 81. This demonstrates that proton tunneling is associated with the rate limiting step.64 More recently, a KIE ∼ 500−700 has been observed22 for the DM, which also displays a much larger effective Arrhenius activation energy, EHa , compared to the WT. As in this study, earlier attempts to understand the SLO kinetics have also used factorized tunneling models to analyze the KIE and its temperature dependence.1,22,30 In some of these approaches,22 the “linear approximation” was invoked within the context of relatively large “effective” harmonic D−A force constants (∼200−500 N/m), which tightly localize the D−A distribution. More recently, the quadratic terms associated with the D−A coupling have been considered,30 still using a harmonic model for D−A motion, and a range of solutions was found that reproduces the WT and DM kinetics. The corresponding interpretation30 was that the DM has a slightly elongated average D−A distance relative to the WT and a similar D−A frequency (i.e., distribution width). Oxygen Activation. Increased negative charge on the active site oxygen atom will move the onset of electronic E


Article


then used to determine Vbare(R) for the SLO system. A test was also done where the FeOH proton was shifted to the Ile839 COO− and the OH bond length of the resulting COOH was fixed in order to prevent the proton from transferring back. In this case, the oxygen atom on the iron developed a charge of −1.08e, which corresponds to δ ∼ 0.85. This results in a inward shift of the Vbare(R) minimum, RDA 0 , from ∼3.35 Å to ∼3.15 Å. Donor−Acceptor Distance Distributions for the WT and DM of SLO. All of the tunneling kinetic observables (within their listed experimental uncertainty) have been calculated using eq 4. The range of protein potentials for δ =

repulsion inward so that smaller D−A tunneling distances can be achieved. In order to quantify this effect, we performed DFT calculations that used a more detailed model of the SLO active site as described in section S2 of the SI. As can be seen in Figure 2, the X-ray structure (PDB: 3PZW)77 places a

0.6, which generate “tunneling-active” P⟨iDA*⟩(R) distributions that fit (within error) the four experimental observables in Table 1, are listed in Table 2. Tables for δ = 0.8 can be found in the SI (Tables S5.1 and S5.2). For the DM, use of a conformational *⟩ (R), was sufficient to fit the observed tunneling potential, V⟨iconf *⟩ (R) contains two unknown kinetics. However, because V⟨iconf *⟩ , there is a range of possibilities that parameters, R⟨i0 *⟩ and f⟨iconf can be found as delineated in Table 2. On the other hand, for

Figure 2. Active site of WT SLO from the crystal structure72 (PDB: 3PZW).

the WT there are very strong forces needed to position P⟨iDA*⟩(R) below the van der Waals contact distance, at ∼2.8 Å, which is necessary to account for the observed rates. We therefore used

carboxylate oxygen atom of C-terminal Ile839 only 2.5 Å away from the catalytically active iron- bound oxygen acceptor atom. Thus, the active site hydroxyl proton is almost certainly engaged in a separate strong hydrogen bond to the Ile839 carboxylate oxygen atom, which might facilitate the tunneling reaction by helping to activate the oxygen atom acceptor78 (i.e., |q| can be increased as the H-bond strengthens). Depending on the orientation constraints placed on His 499, the calculation returned values between q≌ −0.8e and q≌ −1.0e for the net charge on iron bound oxygen atom as calculated using a grid-based method (SI section S2). When δ ∼ 0.6, a similar charge (−0.9e) develops on the oxygen atom of the H3CH---[OH]δ− donor−acceptor test system and this was

⟨i*⟩ (R), which contains the eq 2 and the potential, Vcompress *⟩ , as the magnitude of the total compressive force, F⟨icompress only unknown R-dependent parameter. Under this condition a unique distribution emerges from the analysis and determines the magnitude of the compressive force (within the error bounds of the experiments).

The first two moments of P⟨iDA*⟩(R) are listed in Table 2 and they can be used for comparison with the usual harmonic approach. The first moment of P⟨iDA*⟩(R) generates the average distance, R, between the carbon donor of the substrate and the

Table 2. Rate Expression Parameters for WT and DM SLO Tunneling Kineticsa parameters bare potential CH---[OH] δ

WT

DM

0.6

0.6

δ−

protein force/potential

moments of P⟨iDA*⟩ (R) ⟨R⟩ (Å) σeff (Å) absolute scaling ΔH‡s + ΔH⟨i*⟩ (kcal/mol)b ⎡ ΔS⟨i*⟩ ⎤ |VRP|2 exp⎣⎢ k ⎦⎥ ([cm−1]2)c 5 B

*⟩ F⟨a compress (N/m) 1.25−1.35

2.85−2.83 0.091−0.088

0.74−0.84

*⟩ V⟨b conf (R) ⟨b*⟩ fconf (N/m) 1−25 *⟩ R⟨b (Å) 4.4−2.8 0

3.78−3.03 0.46−0.089

5.0−7.6 6.4−1.4 × 104

1.5−1.1

a The parameter ranges for the WT a* state conformation reflect the experimental error, while for the DM b*-state they reflect the possibility of multiple distributions consistent with experiment. The parameter ranges are ordered consistently throughout the table. bThe approximation ΔH‡s = (ΔH0 + λs)2/4λs is used to find the environmental reorganization energy, λs, where ΔG0 ≅ ΔH0 has been estimated using experimental data1,25 to be −5.4 kcal/mol. For WT we assume ΔH⟨a*⟩ ≅ 0 so that λs = 11.1−11.7 kcal/mol. The value of λs for DM may be slightly different and it depends on the value of ΔH⟨b*⟩ that is chosen. However, variation of λs has a very small effect. If we take ΔH⟨b*⟩ ≅ 0 for DM, the range of λs is 5.4−6.3 ⎡ ΔS⟨i*⟩ ⎤ |VRP|2 kcal/mol Taking ΔH‡s (DM) ≅ ΔH‡s (WT), as done previously,22 yields λs ≅ 3.3 kcal/mol . cThe values of 5 exp⎢⎣ k ⎥⎦ were selected to B reproduce the absolute magnitude of the experimental rates1,22 at 30 °C.

F


Article


Figure 3. Panel (a) shows an example probability distribution for WT SLO that is consistent with the tunneling kinetics1,22 and ENDOR results.31 The active WT a* -state (green) is shown with a ∼ 1% population weighting. Panels (b)−(d) examine the interconversion between the a-state (normalized distribution in orange) and the tunneling-active a* -state (shown as a normalized distribution in green). Panel (b) shows a situation where a strong conformational force, ∼ 1.7 nN, compresses the average D−A distance below van der Waals contact. The minimum protein conformational force ∼1.3 nN is shown in panel (c). Panel (d) presents a scenario where a strong electric field, ER ∼ 60−100 MV/cm (for q ≅ −1.0 *⟩ to −0.6e), helps to generate an overall compressive force, F⟨a compress ≅ 1.3 nN that shifts the electronic repulsion inward allowing for more reasonable conformational forces (≲300 pN). The a-state is also shown with ⟨R⟩ ≅3.1 Å near the van der Waals contact distance determined by ENDOR measurements.

minimum protein compressive force was still on the ∼1 nN scale (see section S5.2 of SI). The large magnitude of these forces indicates that strong protein electric fields are needed to define the tunneling-active protein conformation of WT SLO. For the DM, on the other hand, such strong compressive forces are not necessarily required. Although various possibilities exist to explain the DM tunneling kinetics, one limit involves a broad D−A distribution, centered significantly above the van der Waals contact distance, where only small conformational forces are required. At the other limit, the D−A distribution for DM is narrow and strong electric fields are again required. In order to determine the ranges for the DM in Table 2, we set an upper limit for the reduced force constant of the protein

active site oxygen acceptor. The second moment is the reduced variance of the distribution, which we denote as σ2eff (for the harmonic model this is denoted as σ2R). Near room temperature, the variance can be used to define an effective force constant, kT feff ≡ B2 , and this is helpful when comparing the anharmonic σ eff

distribution to models that use a simple harmonic D−A “vibrational” potential and its associated frequency.40 In the absence of unusually strong (≳300 pN) D−A compressive forces, the protein conformational potential *⟩ (R), can be used to dress the D−A interaction. energy, V⟨iconf However, when very strong compressive forces are required, *⟩ the electric field potential, V⟨ifield (R), should be applied. For the WT SLO tunneling kinetics, we were unsuccessful in finding a protein conformational force that fell below what we consider to be a generous upper limit for such forces (∼300 pN). In fact, the minimum protein conformational force was found to be ∼1.3 nN when δ = 0.6 and ∼1 nN when δ = 0.8. Moreover, *⟩ (R), were examined, they when the associated potentials, V⟨iconf were found to be nearly linear across the width of the D−A *⟩ distribution (e.g., see Figure 3b,c), so that F⟨iconf is approximately constant. Because forces in the nN range exceed, by nearly an order of magnitude, the conformational forces observed in proteins,57−61 the electric field potential was also included as prescribed by eq 2. Thus, even though it derives from two ⟨i*⟩ becomes the single unknown parameter sources, Fcompress needed to establish the tunneling-active D−A distribution for WT SLO. As noted above, when an extremely activated oxygen charge was created by moving the FeOH proton to Ile839, the

*⟩ conformational potential, f⟨b conf ≤ 25 N/m, which is typical of a strong hydrogen bond. More details, including direct fits to the tunneling rates can be found in section S5.2 of the SI. We note that, when fitting the SLO DM kinetic data, the use of the anharmonic D−A potential generally allows for the possibility of D−A distributions that have larger R and σeff compared with models using a harmonic D−A potential.22,30 Furthermore, the values of R and σeff can be significantly larger for the DM than for the WT as seen by the upper limit in Table 2 where ⟨R⟩ = 3.78 Å and σeff = 0.46 Å and such effects have been seen in the single mutants of SLO,30,36,76 as well as in other enzymes.24 An alternative limit for the DM can also be seen in Table 2, where ⟨R⟩ = 3.03 Å and σeff = 0.089 Å. Analysis of the SLO tunneling kinetics using harmonic D−A potentials22,30 focused on this limit for DM, which suggests little change in flexibility between DM and the WT because of their similar D−A distribution widths.

G


Article


Waals contact “doorway” a-state conformation and the inactive b-state conformation. The population weighting for the a*-state conformation is chosen to be ≲1% to demonstrate how the tunneling-active state with ⟨R⟩ ∼ 2.8 Å can remain unresolved in the wing of the majority a-state population. Based on estimates discussed in section S7 of the SI, the population of the a*-state could actually be significantly less than 1%. The nature of the tunneling-active a*-state is explored in more detail in panels b−d of Figure 3. The possibility that the a*-state is created by a nonelectrostatic protein conformational force alone leads to a very large value for the protein conformational compressive force. An example, using only *⟩ V⟨a conf and Vbare for δ = 0.6, is given in Figure 3b, where the

Proposed D−A Distributions and the ENDOR Observations. The parameters listed in Table 2 indicate D−A equilibrium distance of ⟨R⟩ = 2.84 Å for the WT and a range of ⟨R⟩ = 3.03−3.78 Å for the DM. Recent MD calculations31 yield C11---O distances of ⟨R⟩ ∼ 3.1 Å and ⟨R⟩ ∼ 3.8 Å for the two ENDOR-visible global conformations of WT and ⟨R⟩ ∼ 3.9 Å for the DM. Importantly, the ENDOR line shapes for these global conformations indicate31 that the corresponding widths of the D−A distributions are very similar for DM and WT. As can be seen in Table 2, the distribution width for WT is σeff ∼ 0.09 Å, while for the DM with ⟨R⟩ ≌ 3.8 Å (which is very close to the ENDOR distance), the width (σeff ∼ 0.46 Å) is significantly increased. Thus, in order to remain consistent with the ENDOR measurements,31 a tunneling-inactive DM conformation with a distribution width similar to WT and ⟨R⟩ = 3.9 Å must exist. We label this the “b-state”, following the ENDOR notation.31 We postulate that the majority of the DM population is trapped (or stuck) in this state and its relatively narrow distribution width is associated with restricted relative motion between the metal and the C11 carbon atom of the substrate. We further suggest that this tunneling-inactive global conformation is in equilibrium with a tunneling-active global conformation, or “b*state”, which in one limit can sample a much wider range of D− A distances due to additional conformational flexibility. For example, one possibility is that the tunneling-active b*-state conformation for DM has a distribution where ⟨R⟩ = 3.8 Å and σeff= 0.46 Å, as found in Table 2. We note that the presence of a small population of the disordered b*-state conformation will add only a broad background to the ENDOR line shape, making direct detection of this conformation unlikely. On the other hand, if the smaller limiting value of ⟨R⟩ ∼ 3.0 Å is selected22,30 for the b*-state conformation, the probability distribution narrows to σeff = 0.089 Å (Table 2). Because of its narrower width, choosing this limit for the b*-state conformation means that it must have an even smaller population in order to escape ENDOR detection.31 Such a narrow distribution also requires a large compressive electric field force, similar to what is found for WT (see below). The WT tunneling kinetics have been discussed in the context of two ENDOR-detected states with differing average D−A distances.31 The “a” state is located at ⟨R⟩ = 3.1 Å and the “b” state is located at ⟨R⟩ = 3.8 Å. The b-state of WT is very similar to the b-state of DM and both are inactive because the D−A distances are too large to account for the tunneling kinetics of either sample.31 On the other hand, the a-state has an average D−A distance that is close to what is expected for a van der Waals contact pair between D and A, but it is still too large to quantitatively account for the WT tunneling kinetic observables. Thus, we treat the a-state as a “doorway” state that allows passage to the WT tunneling-active conformation that we denote here as the a*-state with ⟨R⟩ ≌ 2.84 Å, falling significantly below the 3.1−3.3 Å van der Waals contact distance. The a*-state distribution in Table 2 is similar to what was found in a previous analysis using a harmonic D−A potential.30 It is evidently populated with relatively low probability via a conformational search31 because it is not detected as an independent state in the ENDOR analysis (e.g., it is hidden in the background). In Figure 3a we display an example D−A probability distribution for the WT of SLO that is consistent with both the available kinetic data and the ENDOR line shapes. The a*-state conformation is visualized in green, along with the van der

*⟩ ⟨a*⟩ conformational force, f⟨a − RDA conf (R0 0 ), is found to be ∼1.7 nN. However, as shown in Figure 3c, we also searched for the minimum possible magnitude of the conformational force and it was found to be ∼1.3 nN. Because this is still too large to be assigned to a protein conformational force, we suggest that the tunneling-active a*-state conformation must involve a protein electric field32 acting on the CH---Oδ− subsystem. This situation *⟩ can be simulated by dressing Vbare(R) with V⟨a field (R) as noted in eq 2 (and in section S5.3 of SI) so that the electronic repulsion is suppressed and the D−A equilibrium distance is reduced (see insert of Figure 1). The field-strength, ER ∼ 60−100 MV/cm shown in panel 3d corresponds to an oxygen acceptor charge in the range −0.6e to −1.0e. These fields are very similar to those found in recent Stark shift experiments using infrared vibrational probes in enzyme active sites.32 We conclude that the tunneling-active a*-state conformation of WT SLO involves protein induced electric field enhancement of the D−A binding interaction. However, it is likely that protein conformational forces act to maintain the a-state van der Waals “doorway” configuration. Thus, it appears that active site oxygen charge, polarizing electric fields, and protein compressive forces all play important roles in the PCET tunneling mechanism of SLO. The broad D−A distribution limit found in Table 2, for the tunneling-active b*-conformation of the DM, is displayed in Figure 4. The normalized tunneling-inactive (b-state) and tunneling-active (b*-state) distributions are shown as red and green curves in Figure 4a,b, respectively. In Figure 4c, an example is given that uses an arbitrary weighting factor for the tunneling-active b*-state conformation. The value of p⟨b*⟩ is chosen to be small enough (p⟨b*⟩ ≲ 0.1) that the population of the broad DA distribution falls below the detection threshold of the ENDOR measurements.31 The dominant inactive b-state DM conformation (p⟨b*⟩≌0.9) is located at ⟨R⟩ ≌ 3.9 Å and has a width that is similar to the WT a-state and b-state conformations as required by the ENDOR observations.

■

DISCUSSION Analytical Models for PCET in Enzymes. The analysis presented here helps to improve the analytical nonadiabatic tunneling model27,28,38,43 and can be used by experimentalists when interpreting kinetic measurements of PCET reactions. A more accurate treatment of proton tunneling will result when the potential includes the electronic repulsion between donor and acceptor atoms because most of the tunneling amplitude is found in the short distance tail of the D−A distribution. The D−A repulsion acts oppositely to anharmonicity in the proton coordinate and these two effects can sometimes nearly cancel as discussed in more detail in section S6 of the SI. The “bare” D− H


Article


for WT. This is not based on a random choice of scaling factors. Rather it is a condition that is imposed by the absolute rate measurements. Note that this observation is independent of The DM distribution limits on ⟨R⟩ and σeff, as can be seen from the full range of possibilities given in Table 2 (see also Tables S5.1 and S5.2 in SI). When the tails of the D−A distributions for DM and WT are compared, there is only a ∼ 0.2 Å difference in the region where 99% of proton tunneling takes place. Most of the tunneling occurs in the tail region between ∼2.5−2.7 Å for WT (see Figure S5.1 in SI), while in DM most of the tunneling occurs between ∼2.6−3.0 Å. The electron tunneling distance is generally much larger (e.g., ∼ 6.5 Å via a through-bond route79 from the pi electrons on adjacent carbon atoms (C10 or C12) through the C11H proton donor interacting with the active site FeO acceptor). Thus, the distance of the electron tunneling path should be essentially immune to the small sub-Å differences in the proton tunneling distance found for WT and DM. If anything, one might expect that the electron tunneling path for DM is slightly larger than that of WT. However, the calculation of |VRP|2 is complicated by the concerted nature of the proton and electron tunneling, which makes their relative time scales important.80 Thus, many of the lessons learned about the coupling associated with long-range electron tunneling in proteins81−83 may not apply under these circumstances.80 |V |2 ⎡ ΔS⟨i*⟩ ⎤ When the magnitudes of 5RP exp⎣⎢ k ⎦⎥ given in Table 2 are B combined with independent calculations80 of VRP ∼ 1600 cm−1 they yield estimates for the population factors, p⟨i*⟩. In the case of WT the value of p⟨a*⟩ ∼ 10−6−10−7 (under the assumption ΔH⟨a*⟩ ∼ 0). This validates the proposition that the tunnelingactive a*-state for WT SLO is not observable using the ENDOR technique. In fact, the very small value indicated for the tunneling-active population suggests that it might also be difficult to directly characterize a rare structural fluctuation into the a*-state using molecular modeling simulations. A more detailed discussion of the extracted population factors for the WT and DM can be found in S7 of the SI. Assuming that VRP(DM) ≈ VRP(WT), based on similar electron tunneling distances, the difference in the relative |V |2 ⎡ ΔS⟨i*⟩ ⎤ magnitude of the 5RP exp⎢⎣ k ⎥⎦ terms can be assumed to arise B from the entropic factors that are important in determining the population of the tunneling-active states. In the case of DM, we note that the b*-state, associated with the broad (σeff ≌ 0.46 Å) distribution, is located at ⟨R⟩ ≌ Å 3.78 Å, which is very close to the b-state ENDOR distance.31 This could be a coincidence, but it suggests that the broad distribution limit for the b*-state should be considered as a candidate for the tunneling-active state of the DM. Moreover, when the broad DM b*-state distribution (σeff ≌0.46 Å) is compared to the b-state distribution detected by ENDOR (which has a narrow width similar to the a-state, σeff ≌ 0.09 Å), it suggests that the b*-state has more disorder and a significantly larger entropy; ΔS⟨b*⟩ ≈ S⟨b*⟩ − S⟨b⟩ > 0. This choice for the b*-state allows sampling of a wider range of tunneling distances where the short-distance tail of this broad distribution enables the (very slow) proton tunneling reaction observed for the DM. On the other hand, the WT enzyme apparently evolved to undergo an entropically controlled search 69 where it interconverts to a rare tunneling-active D−A conformation (a*-state) with aligned electric fields that allow ⟨R⟩ to shift to

Figure 4. Panel (a) shows a normalized probability distribution for the tunneling-inactive b-state associated with the DM of SLO. The b-state is ENDOR-visible31 and has a width and average D−A distance similar to the WT b-state. Panel (b) shows one of the possible tunnelingactive b*-state distributions, consistent with the DM tunneling kinetics.22 Panel (c) shows the full probability distribution, following eq 3b, where the b*-state is weighted by p⟨b*⟩ ≅ 0.1 and the b-state by p⟨b⟩ ≅ 0.9. All distributions use δ = 0.6. The b-state uses R⟨i⟩ 0 = 4.0 Å and f⟨i⟩ conf = 20 N/m, while the b*-state uses the parameters from Table 2 for ⟨R⟩ = 3.78 Å.

A anharmonic potential is found using quantum chemical calculations, which can then be “dressed” with classical potentials that reflect the protein conformational forces or its internal electric fields. Magnitude of the Prefactors Associated with the DM and WT Tunneling Rates. When Table 2 is examined across the range of acceptable P⟨iDA*⟩ distributions, it becomes clear that ⎡ ΔS * ⎤ |V |2 the factors 5RP exp⎢⎣ k⟨i ⟩ ⎥⎦ for WT and DM must differ by many B orders-of-magnitude in order to reproduce the absolute magnitude of the tunneling rates, while at the same time fitting the other kinetic observables. Specifically, the value for |VRP|2 ⎡ ΔS * ⎤ exp⎣⎢ k⟨i ⟩ ⎦⎥ for the DM must be ∼104-105 times larger than 5 B I


Article

The Journal of Physical Chemistry B ∼2.8 Å, below the van der Waals contact distance of the associated “doorway” a-state. Thus, the WT tunneling reaction takes place only after accessing a lower entropy conformation with ΔS⟨a*⟩ ≈ S⟨a*⟩ − S⟨a⟩ < 0. The opposing signs of the entropic population factors for the WT and DM tunnelingactive states appear naturally within the scenario of a broad distribution for the DM tunneling-active conformation and this helps to explain the disparity in their prefactors. In addition, the choice of the more disordered b*-state for the tunneling active conformation of DM can also help to explain the unusually large value of the effective Arrhenius slope, EHa ≌10 kcal/mol, which is found for the DM in Table 1. The large EHa is consistent with enthalpic trapping of the tunneling-inactive b-state conformation. For example, the substrate in DM could be bound (or “stuck”) at a distance from the active site that is too large to allow tunneling. A larger entropy for the tunneling-active b*-state requires that it also have a higher enthalpy than the b-state because it must have a small population (i.e., not detectable by ENDOR). Thus, there is an enthalpic term, ΔH⟨b*⟩ ≈ H⟨b*⟩ − H⟨b⟩, that will combine with ΔH‡s to increase EHa for the DM. The quantity ΔH⟨b*⟩ could be interpreted as the substrate binding energy for an inactive binding conformation. The other limit in Table 2, with ⟨R⟩ ∼ 3.0 Å and σeff ≌0.089 Å, is an alternative choice for the b*-state of DM. This limit generates a distribution that is similar to what was found in a prior study22,30 using the standard harmonic D−A model and it is discussed further in section S7 of the SI. Electric fields and the Tunneling-Active Conformation of WT SLO. A key point is that the “doorway” a-state conformation of the WT enzyme, found near the van der Waals contact D−A distance ⟨R⟩ = 3.1 Å, becomes tunneling-active by making a rare (p⟨a*⟩ ≪ 1) conformational transition to the a*state, which resides significantly below the van der Waals contact distance.37 We first considered the possibility that protein conformational forces were involved in the formation of the a*-state (see Figure 3b,c). However, the minimum conformational force necessary to compress the D−A distribution to ⟨R⟩ ∼ 2.8 Å against the electronic repulsion forces is ∼1.3 nN, which is much larger than the forces (∼10− 100 pN) normally associated with molecular motors and protein unfolding.60,61 The effects of increased charge on the oxygen atom acceptor and increased anharmonicity along the proton coordinate37 were also explored, but the conformational force needed to establish the a*-state remained in the ∼ nN range. The only apparent explanation is that the a*-state conformation involves oriented electric fields33 that increase the attractive interaction between the substrate CH bond and the active site oxygen atom. In this simple model, such effects are quantified as a function of the field strength (see insert of Figures 1 and 3d) so that the need for a large protein conformational force is eliminated. The upper limit range for protein unfolding forces is on the order of hundreds of pN57−59 so that, depending on the magnitude chosen for Fconf ≲300 pN, we estimate Ffield ∼ 1 nN as can be inferred from Figure 3d. The internal protein electric fields associated with such forces are ∼60−100 MV/cm, depending on the charge associated with the acceptor oxygen atom. Such fields are measurable and they fall in the range of recent experiments using Stark shift spectroscopy in enzymes.32,62,63 It should also be noted that, for classical transitions, the effect of the electric field is usually considered to reduce the reaction barrier;63 however, in the

context of a tunneling reaction mechanism, a more important role for the electric field is to reduce the effective “width” of the barrier.

■

CONCLUSIONS

An often-used1,27,30,34,37,41 assumption in proton tunneling calculations involves a harmonic D−A potential and a Gaussian distribution of D−A distances. The anharmonic potential presented here utilizes DFT-based quantum calculations to model the DH---A interaction, including the steep repulsive wall that prevents the D−A distribution from accessing overly short distances. The conformational and electric field forces of the surrounding protein can also be explicitly included by introducing only a few unknown parameters. The compressive force needed to quantitatively explain the PCET reaction of WT SLO is found to fall in the nN range. A compressive force of this magnitude is much larger than the 10−100 pN forces observed for molecular motors or conformational unfolding.57−61 This indicates that internal protein electric fields are involved in preparing the tunneling-active a*-conformation needed for the PCET reaction. In contrast to the WT, a range of tunneling-active D−A distributions are possible for the DM. If we choose the possibility where the DM samples a broad range of D−A distances, the associated entropic factors can help to account for the 104−105 increase of the DM prefactor relative to the WT, as well as its larger effective Arrhenius slope. Finally, as discussed in more detail in S7 of the SI, the measured absolute rates, along with independent calculations of VRP, suggest that the population factors, p⟨i*⟩, for the tunneling active states of both WT and DM are very small (∼10−5−10−7). Generally, anharmonic effects are most apparent in the short distance “tail” of the D−A distribution (or wave function) where almost all of the tunneling occurs. Such effects are therefore quite important to characterize when attempting to extract physical parameters, such as the position and width of the D−A distance distribution. A more accurate determination of these parameters improves our understanding of the tunneling process and of the catalytic process as a whole. Currently, a harmonic model for the D−A probability distribution is widely used by experimentalists. The simple model presented here includes D−A electronic repulsion and accomplishes this without the need for introducing more unknown parameters. Moreover, when the total compressive force greatly exceeds the maximum-allowed protein conformational force, a single parameter, along with the quantumcalculated bare potential, is sufficient to define the D−A potential and its temperature dependent distribution. The treatment of tunneling in other enzymatic systems, which also exhibit coupled proton and electron transfer,2−6,14 may benefit from an analysis based on the anharmonic model for the D−A interaction. It is likely that this model can be improved by using more sophisticated potentials. However, the basic approach presented here offers a simple method for treating the anharmonic electronic repulsion and interaction forces associated with active site oxygen activity, as well as the electric field and conformational forces associated with the surrounding protein. As such, this simple model should be a useful aid to experimentalists when interpreting the kinetics of chemical reactions that involve tunneling. J


Article


■

(13) Yosca, T. H.; Rittle, J.; Krest, C. M.; Onderko, E. L.; Silakov, A.; Calixto, J. C.; Behan, R. K.; Green, M. T. Iron(Iv)Hydroxide Pk(a) and the Role of Thiolate Ligation in C-H Bond Activation by Cytochrome P450. Science 2013, 342, 825−829. (14) Rittle, J.; Green, M. T. Cytochrome P450 Compound I: Capture, Characterization, and C-H Bond Activation Kinetics. Science 2010, 330, 933−937. (15) Groves, J. T. Key Elements of the Chemistry of Cytochrome-P450 - the Oxygen Rebound Mechanism. J. Chem. Educ. 1985, 62, 928− 931. (16) Klinman, J. P.; Kohen, A. Hydrogen Tunneling Links Protein Dynamics to Enzyme Catalysis. Annu. Rev. Biochem. 2013, 82, 471− 496. (17) Layfield, J. P.; Hammes-Schiffer, S. Hydrogen Tunneling in Enzymes and Biomimetic Models. Chem. Rev. 2014, 114, 3466−3494. (18) Truhlar, D. G. Tunneling in Enzymatic and Nonenzymatic Hydrogen Transfer Reactions. J. Phys. Org. Chem. 2010, 23, 660−676. (19) Klinman, J. P. An Integrated Model for Enzyme Catalysis Emerges from Studies of Hydrogen Tunneling. Chem. Phys. Lett. 2009, 471, 179−193. (20) Davydov, R.; Perera, R.; Jin, S. X.; Yang, T. C.; Bryson, T. A.; Sono, M.; Dawson, J. H.; Hoffman, B. M. Substrate Modulation of the Properties and Reactivity of the Oxy-Ferrous and Hydroperoxo-Ferric Intermediates of Cytochrome P450cam as Shown by CryoreductionEpr/Endor Spectroscopy. J. Am. Chem. Soc. 2005, 127, 1403−1413. (21) Davydov, R.; Chemerisov, S.; Werst, D. E.; Rajh, T.; Matsui, T.; Ikeda-Saito, M.; Hoffman, B. M. Proton Transfer at Helium Temperatures During Dioxygen Activation by Heme Monooxygenases. J. Am. Chem. Soc. 2004, 126, 15960−15961. (22) Hu, S.; Sharma, S. C.; Scouras, A. D.; Soudackov, A.; Carr, C. A. M.; Hammes-Schiffer, S.; Alber, T.; Klinman, J. P. Extremely Elevated Room-Temperature Kinetic Isotope Effects Quantify the Critical Role of Barrier Width in Enzymatic C-H Activation. J. Am. Chem. Soc. 2014, 136, 8157−8160. (23) Bhabha, G.; Lee, J.; Ekiert, D. C.; Gam, J.; Wilson, I. A.; Dyson, H. J.; Benkovic, S. J.; Wright, P. E. A Dynamic Knockout Reveals That Conformational Fluctuations Influence the Chemical Step of Enzyme Catalysis. Science 2011, 332, 234−238. (24) Wang, L.; Goodey, N. M.; Benkovic, S. J.; Kohen, A. Coordinated Effects of Distal Mutations on Environmentally Coupled Tunneling in Dihydrofolate Reductase. Proc. Natl. Acad. Sci. U. S. A. 2006, 103, 15753−15758. (25) Hatcher, E.; Soudackov, A. V.; Hammes-Schiffer, S. ProtonCoupled Electron Transfer in Soybean Lipoxygenase. J. Am. Chem. Soc. 2004, 126, 5763−5775. (26) Olsson, M. H. M.; Siegbahn, P. E. M.; Warshel, A. Simulations of the Large Kinetic Isotope Effect and the Temperature Dependence of the Hydrogen Atom Transfer in Lipoxygenase. J. Am. Chem. Soc. 2004, 126, 2820−2828. (27) Kuznetsov, A. M.; Ulstrup, J. Proton and Hydrogen Atom Tunnelling in Hydrolytic and Redox Enzyme Catalysis. Can. J. Chem. 1999, 77, 1085−1096. (28) Hammes-Schiffer, S.; Stuchebrukhov, A. A. Theory of Coupled Electron and Proton Transfer Reactions. Chem. Rev. 2010, 110, 6939− 6960. (29) Cukier, R. I.; Nocera, D. G. Proton-Coupled Electron Transfer. Annu. Rev. Phys. Chem. 1998, 49, 337−369. (30) Soudackov, A. V.; Hammes-Schiffer, S. Proton-Coupled Electron Transfer Reactions: Analytical Rate Constants and Case Study of Kinetic Isotope Effects in Lipoxygenase. Faraday Discuss. 2016, 195, 171−189. (31) Horitani, M.; Offenbacher, A. R.; Carr, C. A. M.; Yu, T.; Hoeke, V.; Cutsail, G. E.; Hammes-Schiffer, S.; Klinman, J. P.; Hoffman, B. M. 13 C Endor Spectroscopy of Lipoxygenase-Substrate Complexes Reveals the Structural Basis for C-H Activation by Tunneling. J. Am. Chem. Soc. 2017, 139, 1984−1997. (32) Fried, S. D.; Boxer, S. G. Measuring Electric Fields and Noncovalent Interactions Using the Vibrational Stark Effect. Acc. Chem. Res. 2015, 48, 998−1006.

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpcb.7b05570. A description of DFT calculations, rate expressions, example D−A interaction potentials, extraction of fitting parameters, anharmonic cancellation effects, and discussion of the active-state populations (PDF)

■

AUTHOR INFORMATION

Corresponding Author

*Tel: 617-373-2918; E-mail: [email protected]. ORCID

Paul M. Champion: 0000-0001-6958-2836 Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS This work is supported by NSF CHE-1243948. We also thank Yihan Shao for assistance with the electric field calculations using Q-Chem.

■

REFERENCES

(1) Knapp, M. J.; Rickert, K.; Klinman, J. P. Temperature-Dependent Isotope Effects in Soybean Lipoxygenase-1: Correlating Hydrogen Tunneling with Protein Dynamics. J. Am. Chem. Soc. 2002, 124, 3865− 3874. (2) Groves, J. T. Enzymatic C-H Bond Activation Using Push to Get Pull. Nat. Chem. 2014, 6, 89−91. (3) Sikorski, R. S.; Wang, L.; Markham, K. A.; Rajagopalan, P. T. R.; Benkovic, S. J.; Kohen, A. Tunneling and Coupled Motion in the Escherichia Coli Dihydrofolate Reductase Catalysis. J. Am. Chem. Soc. 2004, 126, 4778−4779. (4) Whittaker, M. M.; Whittaker, J. W. Catalytic Reaction Profile for Alcohol Oxidation by Galactose Oxidase. Biochemistry 2001, 40, 7140−7148. (5) Reinhardt, L. A.; Svedruzic, D.; Chang, C. H.; Cleland, W. W.; Richards, N. G. J. Heavy Atom Isotope Effects on the Reaction Catalyzed by the Oxalate Decarboxylase from Bacillus Subtilis. J. Am. Chem. Soc. 2003, 125, 1244−1252. (6) Su, Q. J.; Klinman, J. P. Probing the Mechanism of Proton Coupled Electron Transfer to Dioxygen: The Oxidative Half-Reaction of Bovine Serum Amine Oxidase. Biochemistry 1998, 37, 12513− 12525. (7) Saraste, M. Oxidative Phosphorylation at the Fin De Siecle. Science 1999, 283, 1488−1493. (8) Belevich, I.; Verkhovsky, M. I.; Wikstrom, M. Proton-Coupled Electron Transfer Drives the Proton Pump of Cytochrome C Oxidase. Nature 2006, 440, 829−832. (9) Okamura, M. Y.; Paddock, M. L.; Graige, M. S.; Feher, G. Proton and Electron Transfer in Bacterial Reaction Centers. Biochim. Biophys. Acta, Bioenerg. 2000, 1458, 148−163. (10) Schlichting, I.; Berendzen, J.; Chu, K.; Stock, A. M.; Maves, S. A.; Benson, D. E.; Sweet, R. M.; Ringe, D.; Petsko, G. A.; Sligar, S. G. The Catalytic Pathway of Cytochrome P450cam at Atomic Resolution. Science 2000, 287, 1615−1622. (11) Madrona, Y.; Hollingsworth, S. A.; Khan, B.; Poulos, T. L. P450cin Active Site Water: Implications for Substrate Binding and Solvent Accessibility. Biochemistry 2013, 52, 5039−5050. (12) Imai, M.; Shimada, H.; Watanabe, Y.; Matsushima-Hibiya, Y.; Makino, R.; Koga, H.; Horiuchi, T.; Ishimura, Y. Uncoupling of the Cytochrome P-450cam Monooxygenase Reaction by a Single Mutation, Threonine-252 to Alanine or Valine: A Possible Role of the Hydroxy Amino Acid in Oxygen Activation. Proc. Natl. Acad. Sci. U. S. A. 1989, 86, 7823−7827. K


Article

The Journal of Physical Chemistry B (33) Warshel, A.; Sharma, P. K.; Kato, M.; Xiang, Y.; Liu, H. B.; Olsson, M. H. M. Electrostatic Basis for Enzyme Catalysis. Chem. Rev. 2006, 106, 3210−3235. (34) Borgis, D.; Hynes, J. T. Dynamic Theory of Proton Tunneling Transfer Rates in Solution - General Formulation. Chem. Phys. 1993, 170, 315−346. (35) Mincer, J. S.; Schwartz, S. D. Rate-Promoting Vibrations and Coupled Hydrogen-Electron Transfer Reactions in the Condensed Phase: A Model for Enzymatic Catalysis. J. Chem. Phys. 2004, 120, 7755−7760. (36) Meyer, M. P.; Klinman, J. P. Modeling Temperature Dependent Kinetic Isotope Effects for Hydrogen Transfer in a Series of Soybean Lipoxygenase Mutants: The Effect of Anharmonicity Upon Transfer Distance. Chem. Phys. 2005, 319, 283−296. (37) Siebrand, W.; Smedarchina, Z. Temperature Dependence of Kinetic Isotope Effects for Enzymatic Carbon-Hydrogen Bond Cleavage. J. Phys. Chem. B 2004, 108, 4185−4195. (38) Benabbas, A.; Salna, B.; Sage, J. T.; Champion, P. M. Deep Proton Tunneling in the Electronically Adiabatic and Non-Adiabatic Limits: Comparison of the Quantum and Classical Treatment of Donor-Acceptor Motion in a Protein Environment. J. Chem. Phys. 2015, 142, 114101. (39) Salna, B.; Benabbas, A.; Sage, J. T.; Van Thor, J. J.; Champion, P. M. Wide-Dynamic-Range Kinetic Investigations of Deep Proton Tunnelling in Proteins. Nat. Chem. 2016, 8, 874−880. (40) Soudackov, A. V.; Hammes-Schiffer, S. Nonadiabatic Rate Constants for Proton Transfer and Proton-Coupled Electron Transfer Reactions in Solution: Effects of Quadratic Term in the Vibronic Coupling Expansion. J. Chem. Phys. 2015, 143, 194101. (41) Kiefer, P. M.; Hynes, J. T. Kinetic Isotope Effects for Nonadiabatic Proton Transfer Reactions in a Polar Environment. 1. Interpretation of Tunneling Kinetic Isotopic Effects. J. Phys. Chem. A 2004, 108, 11793−11808. (42) Salna, B.; Benabbas, A.; Champion, P. M. Proton-Coupled Electron Transfer and the ″Linear Approximation″ for Coupling to the Donor-Acceptor Distance Fluctuations. J. Phys. Chem. A 2017, 121, 2199. (43) Kiefer, P. M.; Hynes, J. T. Kinetic Isotope Effects for Nonadiabatic Proton Transfer Reactions in a Polar Environment. 2. Comparison with an Electronically Diabatic Description. J. Phys. Chem. A 2004, 108, 11809−11818. (44) Kiefer, P. M.; Hynes, J. T. Adiabatic and Nonadiabatic Proton Transfer Rate Constants in Solution. Solid State Ionics 2004, 168, 219− 224. (45) Tang, J.; Lee, M. T.; Lin, S. H. Effects of the Duschinsky ModeMixing Mechanism on Temperature Dependence of Electron Transfer Processes. J. Chem. Phys. 2003, 119, 7188−7196. (46) Soudackov, A.; Hatcher, E.; Hammes-Schiffer, S. Quantum and Dynamical Effects of Proton Donor-Acceptor Vibrational Motion in Nonadiabatic Proton-Coupled Electron Transfer Reactions. J. Chem. Phys. 2005, 122, 014505. (47) Hatcher, E.; Soudackov, A.; Hammes-Schiffer, S. Comparison of Dynamical Aspects of Nonadiabatic Electron, Proton, and ProtonCoupled Electron Transfer Reactions. Chem. Phys. 2005, 319, 93−100. (48) Hay, S.; Scrutton, N. S. Good Vibrations in Enzyme-Catalysed Reactions. Nat. Chem. 2012, 4, 161−168. (49) Hammes-Schiffer, S.; Hatcher, E.; Ishikita, H.; Skone, J. H.; Soudackov, A. V. Theoretical Studies of Proton-Coupled Electron Transfer: Models and Concepts Relevant to Bioenergetics. Coord. Chem. Rev. 2008, 252, 384−394. (50) Cha, Y.; Murray, C. J.; Klinman, J. P. Hydrogen Tunneling in Enzyme Reactions. Science 1989, 243, 1325−1330. (51) Kohen, A.; Cannio, R.; Bartolucci, S.; Klinman, J. P. Enzyme Dynamics and Hydrogen Tunnelling in a Thermophilic Alcohol Dehydrogenase. Nature 1999, 399, 496−499. (52) Masgrau, L.; Roujeinikova, A.; Johannissen, L. O.; Hothi, P.; Basran, J.; Ranaghan, K. E.; Mulholland, A. J.; Sutcliffe, M. J.; Scrutton, N. S.; Leys, D. Atomic Description of an Enzyme Reaction Dominated by Proton Tunneling. Science 2006, 312, 237−241.

(53) Edwards, S. J.; Soudackov, A. V.; Hammes-Schiffer, S. Analysis of Kinetic Isotope Effects for Proton-Coupled Electron Transfer Reactions. J. Phys. Chem. A 2009, 113, 2117−2126. (54) Austin, A.; Petersson, G. A.; Frisch, M. J.; Dobek, F. J.; Scalmani, G.; Throssell, K. A Density Functional with Spherical Atom Dispersion Terms. J. Chem. Theory Comput. 2012, 8, 4989−5007. (55) Dunning, T. H. Gaussian-Basis Sets for Use in Correlated Molecular Calculations. I. The Atoms Boron through Neon and Hydrogen. J. Chem. Phys. 1989, 90, 1007−1023. (56) Woon, D. E.; Dunning, T. H. Gaussian-Basis Sets for Use in Correlated Molecular Calculations. 3. The Atoms Aluminum through Argon. J. Chem. Phys. 1993, 98, 1358−1371. (57) Popa, I.; Berkovich, R.; Alegre-Cebollada, J.; Badilla, C. L.; Rivas-Pardo, J. A.; Taniguchi, Y.; Kawakami, M.; Fernandez, J. M. Nanomechanics of Halotag Tethers. J. Am. Chem. Soc. 2013, 135, 12762−12771. (58) Jagannathan, B.; Marqusee, S. Protein Folding and Unfolding under Force. Biopolymers 2013, 99, 860−869. (59) Rief, M.; Gautel, M.; Oesterhelt, F.; Fernandez, J. M.; Gaub, H. E. Reversible Unfolding of Individual Titin Immunoglobulin Domains by Afm. Science 1997, 276, 1109−1112. (60) Rickgauer, J. P.; Fuller, D. N.; Grimes, S.; Jardine, P. J.; Anderson, D. L.; Smith, D. E. Portal Motor Velocity and Internal Force Resisting Viral DNA Packaging in Bacteriophage Φ29. Biophys. J. 2008, 94, 159−167. (61) Andreasson, J. O. L.; Shastry, S.; Hancock, W. O.; Block, S. M. The Mechanochemical Cycle of Mammalian Kinesin-2 Kif3a/B under Load. Curr. Biol. 2015, 25, 1166−1175. (62) Liu, C. T.; Layfield, J. P.; Stewart, R. J.; French, J. B.; Hanoian, P.; Asbury, J. B.; Hammes-Schiffer, S.; Benkovic, S. J. Probing the Electrostatics of Active Site Microenvironments Along the Catalytic Cycle for Escherichia Coli Dihydrofolate Reductase. J. Am. Chem. Soc. 2014, 136, 10349−10360. (63) Fried, S. D.; Bagchi, S.; Boxer, S. G. Extreme Electric Fields Power Catalysis in the Active Site of Ketosteroid Isomerase. Science 2014, 346, 1510−1514. (64) Glickman, M. H.; Klinman, J. P. Nature of Rate-Limiting Steps in the Soybean Lipoxygenase-1 Reaction. Biochemistry 1995, 34, 14077−14092. (65) Srajer, V.; Reinisch, L.; Champion, P. M. Protein Fluctuations, Distributed Coupling, and the Binding of Ligands to Heme-Proteins. J. Am. Chem. Soc. 1988, 110, 6656−6670. (66) Ye, X.; Ionascu, D.; Gruia, F.; Yu, A.; Benabbas, A.; Champion, P. M. Temperature-Dependent Heme Kinetics with Nonexponential Binding and Barrier Relaxation in the Absence of Protein Conformational Substates. Proc. Natl. Acad. Sci. U. S. A. 2007, 104, 14682− 14687. (67) Benabbas, A.; Karunakaran, V.; Youn, H.; Poulos, T. L.; Champion, P. M. Effect of DNA Binding on Geminate Co Recombination Kinetics in Co-Sensing Transcription Factor Cooa. J. Biol. Chem. 2012, 287, 21729−21740. (68) Truhlar, D. G.; Kohen, A. Convex Arrhenius Plots and Their Interpretation. Proc. Natl. Acad. Sci. U. S. A. 2001, 98, 848−851. (69) Nagel, Z. D.; Dong, M.; Bahnson, B. J.; Klinman, J. P. Impaired Protein Conformational Landscapes as Revealed in Anomalous Arrhenius Prefactors. Proc. Natl. Acad. Sci. U. S. A. 2011, 108, 10520−10525. (70) Weintraub, O.; Bixon, M. Entropy of Reaction and Number of Active-Modes in the Theory of Electron-Transfer Reactions. J. Phys. Chem. 1994, 98, 3407−3411. (71) Hopfield, J. J. Electron Transfer between Biological Molecules by Thermally Activated Tunneling. Proc. Natl. Acad. Sci. U. S. A. 1974, 71, 3640−3644. (72) Sutin, N. Theory of Electron-Transfer Reactions-Insights and Hindsights. Prog. Inorgan. Chem. 1983, 30, 441−498. (73) Jortner, J. Temperature Dependent Activation Energy for Electron Transfer between Biological Molecules. J. Chem. Phys. 1976, 64, 4860. L


Article

The Journal of Physical Chemistry B (74) Hanggi, P.; Talkner, P.; Borkovec, M. Rate-Reaction Theory: Fifty Years after Kramers. Rev. Mod. Phys. 1990, 62, 251−341. (75) Edwards, S. J.; Soudackov, A. V.; Hammes-Schiffer, S. Driving Force Dependence of Rates for Nonadiabatic Proton and ProtonCoupled Electron Transfer: Conditions for Inverted Region Behavior. J. Phys. Chem. B 2009, 113, 14545−14548. (76) Meyer, M. P.; Tomchick, D. R.; Klinman, J. P. Enzyme Structure and Dynamics Affect Hydrogen Tunneling: The Impact of a Remote Side Chain (I553) in Soybean Lipoxygenase-1. Proc. Natl. Acad. Sci. U. S. A. 2008, 105, 1146−1151. (77) Chruszcz, M.; Wlodawer, A.; Minor, W. Determination of Protein Structures - a Series of Fortunate Events. Biophys. J. 2008, 95, 1−9. (78) Phatak, P.; Sumner, I.; Iyengar, S. S. Gauging the Flexibility of the Active Site in Soybean Lipoxygenase-1 (Slo-1) through an AtomCentered Density Matrix Propagation (Admp) Treatment That Facilitates the Sampling of Rare Events. J. Phys. Chem. B 2012, 116, 10145−10164. (79) Lehnert, N.; Solomon, E. I. Density-Functional Investigation on the Mechanism of H-Atom Abstraction by Lipoxygenase. J. Biol. Inorg. Chem. 2003, 8, 294−305. (80) Harshan, A. K.; Yu, T.; Soudackov, A. V.; Hammes-Schiffer, S. Dependence of Vibronic Coupling on Molecular Geometry and Environment: Bridging Hydrogen Atom Transfer and Electron-Proton Transfer. J. Am. Chem. Soc. 2015, 137, 13545−13555. (81) Fox, L. S.; Kozik, M.; Winkler, J. R.; Gray, H. B. Gaussian FreeEnergy Dependence of Electron-Transfer Rates in Iridium Complexes. Science 1990, 247, 1069−1071. (82) Winkler, J. R.; Gray, H. B. Long-Range Electron Tunneling. J. Am. Chem. Soc. 2014, 136, 2930−2939. (83) Winkler, J. R.; Gray, H. B. Electron Flow through Metalloproteins. Chem. Rev. 2014, 114, 3369−3380.

M


Tunneling Kinetics and Nonadiabatic Proton-Coupled Electron

Recommend Documents