Why Do Empirical Valence Bond Simulations Yield Accurate Arrhenius Plots?

Computer simulations of the temperature dependence of enzyme reactions using the empirical valence bond (EVB) method have proven to give very accurate results in terms of the thermodynamic activation parameters. Here, we analyze the reasons for why such simulations are able to correctly capture activation enthalpies and entropies and how sensitive these quantities are to parametrization of the reactive potential energy function. We examine first the solution reference reaction for the enzyme ketosteroid isomerase, which corresponds to the acetate catalyzed deprotonation of the steroid in water. The experimentally determined activation parameters for this reaction turn out to be remarkably well reproduced by the calculations. By modifying the EVB potential so that the activation and reaction free energies become significantly shifted, we show that the activation entropy is basically invariant to such changes and that ΔS⧧ is instead determined by the specific mixture of the underlying force fields in the transition state region. The coefficients of this mixture do not change appreciably when the EVB potential is modified within reasonable limits, and hence, the estimate of ΔS⧧ becomes very robust. This is further verified by examining a more complex concerted hydride and proton transfer reaction in the enzyme hydroxybutyrate dehydrogenase.


■ INTRODUCTION
−7 In 2008, we devised a procedure analogous to experimental Arrhenius plots to obtain activation free energies as a function of temperature from molecular dynamics (MD) simulations in combination with the EVB method. 8With high enough precision for the calculated activation free energies, computational Arrhenius plots of either ΔG ⧧ /T vs 1/T or ΔG ⧧ vs T can thus be constructed, which allows the values of ΔH ⧧ and ΔS ⧧ to be reliably extracted by linear regression.Today, this precision can typically be pushed to 0.1−0.2kcal/mol, with many replicate MD/EVB free energy simulations, and the accuracy of the predicted ΔH ⧧ and ΔS ⧧ values is rather determined by how well the Arrhenius plot fits a straight line, than by the precision of the ΔG ⧧ values.For a number of examined enzyme reactions, we found that the calculated thermodynamic activation parameters are in remarkably good agreement with experimental measurements.Some recent examples of this are cytidine deaminase, 9 EF-Tu catalyzed GTP hydrolysis on the ribosome, 10 α-amylase, 4 hydroxybutyrate dehydrogenase, 5 lactate dehydrogenase, 6 and chorismate mutase. 7he most recent case of chorismate mutase is particularly instructive, since it is a unimolecular reaction in the enzyme (and in solution), with no enzymic groups participating in the actual chemistry.In that case, the EVB potential was parametrized on density functional theory (DFT) calculations on the uncatalyzed reaction in solution, with a few explicit water molecules together with a continuum solvent model. 7recisely since the reaction is unimolecular, no thermodynamic corrections to the free energy profile in water, associated with bringing the reactants together in the same solvent cage, are required.The perfectly linear Arrhenius plots obtained from MD/EVB simulations of both the water reaction and the catalyzed reactions by enzymes from two different species turn out to give very accurate values of both activation enthalpies and entropies compared to experimental results. 7he question we want to address in this work is what is the underlying reason for this methodology giving such good estimates of the thermodynamic activation parameters?Here, it should, of course, be kept in mind that any enzyme EVB model has the value of ΔG ⧧ somehow parametrized but, notably, not ΔH ⧧ and ΔS ⧧ .This parametrization was initially most often performed against experimentally derived data for an uncatalyzed reference reaction in water. 1,2More lately, and for reference reactions lacking experimental data, the water reaction has been parametrized based on quantum chemical calculations with a suitable representation of the surrounding solvent. 3,4,7In both of these cases, it should be noted that no parametrization against enzyme data has been done, and the catalytic effect of "moving" the reference reaction into an enzyme is entirely a predicted result of the MD/EVB free energy calculations.An alternative that we have used recently for cases where the reference reaction would be rather complex is to parametrize the EVB potential energy surface directly on DFT calculations for the enzyme reaction.This can be done using either large cluster models of the active site with a continuum solvent model 9 or using QM/MM calculations with explicit solvation. 5,11In order for such approaches to be viable, it is, of course, necessary that the DFT calculations give energetics that reasonably well match experimental results.
Among the examples given above, we have representatives of all four variants of EVB parametrization (solution experiments, 10 solution DFT, 3,4,7 enzyme cluster DFT, 9 and enzyme QM/MM DFT) 5,6 and the predicted values of ΔH ⧧ and ΔS ⧧ are in all cases remarkably good.

■ THEORETICAL BACKGROUND
EVB Methodology.To understand the origin of the thermodynamic activation parameters in the EVB model, it is useful to first review its basic features.An elementary chemical reaction step is usually treated as a two-state problem where the system is then represented by a 2 × 2 Hamiltonian matrix Here, c is the eigenvector and the ground-state energy (E g ) is obtained as the lowest eigenvalue of the secular equation pertaining to the above Hamiltonian: where the diagonal matrix elements H 11 and H 22 are classical force fields representing the reactant and product states of the reaction step.The off-diagonal element H 12 = H 21 represents the coupling or mixing of the two states and is often taken as a constant term, although it can also be given a distance or energy gap dependent form. 1,2,12As said, H 11 and H 22 are standard force fields of the type where the bonded terms (U bonded ) consist of bond (Morse potentials), angle and torsion terms, while the nonbonded terms (U nonbonded ) denote electrostatic and van der Waals (Lennard-Jones) interactions as usual.There is, however, something missing here.Consider for example the simple proton transfer between two water molecules, H 2 O + H 2 O ⇌ OH − + H 3 O + .If we keep the two molecules on each side of the equation at their equilibrium geometries and at infinite separation from each other (meaning no nonbonded interactions between them), then the force field energies would be U 1 (R) = U 2 (R) = 0.This is clearly not correct since there is an absolute energy difference between the left-and right-hand sides of the equation under such a condition.That is, there is a real difference in the heat of formation between 2H 2 O and OH − + H 3 O + that is not accounted for by standard force fields.Hence, we must include such a constant term, and this is, e.g., done by modifying U 2 (R) to U 2 (R) + Δα, where Δα is formally the absolute energy difference between the reactant and product when the reacting groups are at their equilibrium geometries and noninteracting.
We now thus have two unknown parameters in our EVB model that are not given by any force field, namely, H 12 and Δα.This is also the beauty of the method, that these two parameters can be used to calibrate the energetics of the EVB model so that it reproduces whatever data are used to parametrize the model.In practice, it is the activation free energy ΔG ⧧ and the reaction free energy ΔG 0 that are parametrized by parameters H 12 and Δα.To clarify the concepts a bit further, let us consider an even simpler test case, the transfer of an electronic charge between a sodium atom and a sodium ion at fixed distance in water (Na + Na + ⇌ Na + + Na), with H 12 = 0 and Δα = 0.This is thus just a hypothetical transformation without consideration of electron tunnelling or anything like that.We use the free energy perturbation (FEP) method to map the free energy change upon moving from reactant to product, via a mapping potential expressed as a linear combination of the two end-states: The free energy of the system for any value of the mapping parameter λ m relative to the reactant state (λ = 0) is given by summation of Zwanzig's formula: 13 where λ has now been discretized into a number of points λ n (or windows) and the MD ensemble average at each such point estimates the free energy change associated with moving to the next point.The free energy change as a function of λ going from 0 to 1 is shown in Figure 1a together with the free energy change on the ground-state potential given by eq 2. This ground-state free energy profile is obtained with the umbrella sampling formula from the free energy on the mapping potential ΔG(λ m ) according to where the generalized reaction coordinate is the energy gap occurring in eq 2, X s = ΔU = U 1 − U 2 .Since each value of λ m only samples a certain range of the energy gap, the overall free energy profile is obtained by piecing together these ranges according to 12 G X p G X ( ) ( , ) where the sum runs over all the sampling windows (λ m ) that contribute to the bin X s and p m is the normalized statistical weight of the m th sampling window in that bin.This gives the solid black ground-state free energy curve shown in Figure 1a, with a free energy barrier of ΔG ‡ = 37 kcal/mol that is entirely due to solvent reorganization.Likewise, the free energy functions corresponding to the pure reactant and product states can be obtained by the equivalent of eq 6 where E g is replaced by either where i = 1 and 2 and pieced together again as before These diabatic free energy curves are shown with blue and red circles in Figure 1a.As can be seen, ΔG 1 and ΔG 2 intersect at precisely 37 kcal/mol, and this is also the value of ΔG(λ m ) at λ = 0.5.Now, what happens if we introduce a coupling between the two states in terms of a nonzero value of H 12 ?The result of introducing such a coupling for our Na + Na + ⇌ Na + + Na test case with H 12 = 20 kcal/mol is shown in Figure 1b.As can be seen, neither ΔG(λ m ) nor the diabatic ΔG i 's have changed, but the adiabatic ground-state free energy function now has a lower barrier of ΔG ‡ = 20 kcal/mol due to the coupling between reactant and product states, while ΔG 0 is still zero.If we instead keep H 12 = 0 and change our Δα from zero to 20 kcal/ mol, all that happens to our free energy curves is that the product state is lifted by 20 kcal/mol as can be seen from ΔG(λ m ) and ΔG 2 in Figure 1c, and the reaction free energy thus becomes ΔG 0 = 20 kcal/mol.As the diabatic free energy curve of the product state becomes shifted upward, the intersection between ΔG 1 and ΔG 2 will also be up-shifted, and the barrier on ΔG g now becomes 47 kcal/mol.Finally, to illustrate how the EVB potential can be parametrized to any target values of ΔG ⧧ and ΔG 0 , Figure 1d shows the effect of setting the parameters to H 12 = 20 kcal/mol and Δα = 20 kcal/ mol.Now we get ΔG ‡ = 31 kcal/mol and ΔG 0 = 19 kcal/mol on the ground-state free energy curve.The calculation of the ground-state free energy profile according to eqs 6 and 7 can be further illustrated graphically, as shown in Figure 2, where it can be seen how different values of the mapping parameter λ will sample different regions of the energy gap ΔU.It should be noted here that all the dynamics used for sampling the system are run on the artificial potential U map (λ).However, with a parametrized EVB model at hand, it is straightforward to also run the dynamics on the actual ground-state potential.This may be of interest for studying the true dynamic properties of the reacting system, for example, for evaluating the transmission coefficient in transition state theory.The forces acting on the particles of the system will then depend on the energy gap and the off-diagonal matrix element according to the derivatives: where x i is any given Cartesian coordinate, and H 12 is assumed to be a constant for simplicity (with a distance or energy gap dependent H 12 , the expression becomes slightly more complicated).
Computational Arrhenius Plots from MD/EVB Simulations.As an example of how the enthalpy and entropy components of the activation free energy are extracted from MD/EVB simulations, let us again consider the simple Na + Na + ⇌ Na + + Na test case with H 12 = 0 and Δα = 0. We thus repeat the free energy calculations at a set of different temperatures, e.g., in the range 280−320 K in steps of 10 K, and plot ΔG ⧧ /T vs 1/T or, alternatively, ΔG ⧧ vs T. Depending on the magnitude of ΔH ⧧ relative to ΔS ⧧ the linear regression will be more accurate with a higher R 2 value for either of the two types of plots.That is, for a reaction with a high enthalpy component, the ΔG ⧧ /T vs 1/T plot will be the most accurate since the enthalpy then corresponds to the slope of the plot, as opposed to the ΔG ⧧ vs T variant where the entropy becomes the slope.This is illustrated for our test case in Figure 3, where it can be seen that ΔH ‡ = 40.7 kcal/mol and ΔS ‡ = 0.0134 kcal/mol/K, corresponding to TΔS ‡ = +4.02kcal/mol at 300 K. Hence, the free energy barrier for transfer of an electronic charge from Na to Na + is completely dominated by the enthalpy barrier, while the entropy contribution is favorable at the transition state.This can be understood in terms of the strong electrostatic interactions between water and a solvated ion that partly have to break in the transition state region, where the positive charge becomes delocalized between the two Na atoms.As a solvated alkali ion also has an unfavorable entropy contribution to the solvation free energy of about TΔS sol = −6 kcal/mol, 14 the delocalization of the charge reduces this entropy penalty, due to weaker ion−water interactions, and the contribution to the reaction barrier (−TΔS ‡ < 0) thus becomes favorable.
Below we will examine the performance of this approach for calculating thermodynamic activation parameters for a number of test cases, including reactions both in solution and in eznymes.We are particularly interested in assessing the influence of EVB parameters H 12 and Δα on the enthalpy− entropy balance of the activation free energy.Our analysis demonstrates that the EVB parameters, together with the target value of ΔG ⧧ , largely dictate the activation enthalpy, while the entropic contribution is rather intrinsic to the force field.That is, the entropy change upon moving from the reactants to the particular mixture of U 1 and U 2 representing the transition state is found to be almost invariant, with respect to the EVB parameters.

Simulations of the Water Reference Reaction for Ketosteroid Isomerase.
As an example of a relatively simple solution reaction, we choose the acetate catalyzed proton abstraction from 5-androstene-3,17-dione (Figure 4a).This is the solution reference for the first step of the ketosteroid isomerase reaction, where an aspartate residue abstracts the 4β proton from the steroid substrate.Houck and Pollack have determined the activation parameters for this reaction in water both in the forward and backward directions. 15,16We thus simulate this reaction in a 40 Å diameter water sphere with the MD/EVB methodology described above and the standard OPLS-AA/M force field, 17 where the experimentally determined values of ΔG ⧧ = 22.0 and ΔG 0 = 10.9 kcal/mol at 298 K are used as calibration target values for the EVB model.The parametrization against experimental data also means that any nuclear (proton) tunnelling effects are implicitly included in the effective EVB potential. 18ere, however, we encounter the typical problem of simulating bimolecular reactions in solution where the experimental standard state is taken as 1 M.In microscopic terms, reaction rates with this standard state also include the probability of an encounter between the donor and acceptor atoms in addition to the reaction rate for the contact complex.If we set our origin at one of the reacting partners, a 1 M concentration of the other reactant around it means that it is free to move within 1661 Å 3 , or within a radius of about 7.3 Å.Hence, for the encounter probability to be accurately treated, it would be necessary for the MD simulation at the reactant state to correctly sample the diffusion within this 1661 Å 3 sphere.Since such a calculation of what is mostly entropy is known to converge very slowly, 19 one instead usually employs a formal correction to the 1 M reaction rates.That is, a distance restraint that keeps the donor and acceptor in contact with each other is used, and the formal concentration of the  The experimental activation free energies and entropies have been corrected to represent a 55 M standard state, as described in the text.Error bars for the calculated free energies denote the standard error of the mean (s.e.m.) at 300 K from 60 replicate simulations, while those for enthalpies and entropies are the asymptotic standard errors from linear regression.otherwise diffusing partner is then considered to be the same as the surrounding solvent.For a reaction in water, this means that the concentration of the diffusing reactant has changed from 1 to 55 M, implying that the measured activation free energy should be corrected by −RT ln(55) = −2.39kcal/mol at room temperature.With this translational entropy correction, our target value for ΔG ⧧ becomes 19.6 kcal/mol at a 55 M standard state, while ΔG 0 is unaffected since the equilibrium constant between reactant and product does not depend on the standard state.A typical calculated reaction free energy profile at 300 K is shown in Figure 4b.
The results from these calculations of acetate catalyzed deprotonation of 5-androstene-3,17-dione at five different temperatures (290, 295, 300, 305, and 310 K) are shown in the Arrhenius plots of Figure 4c and a typical MD snapshot from the transition state region is shown in Figure 4d.At each temperature, 60 independent reaction free energy profiles were calculated, and each such profile corresponds to 1.05 ns of MD simulation, where 21 distinct λ windows (eq 4) were used.As can be seen, the Arrhenius plots are highly linear with R 2 values above 0.93, and the extracted thermodynamic activation parameters are summarized in Table 1.It is immediately apparent that the enthalpies and entropies for both forward and backward reactions as well as for the equilibrium constant are reproduced by the MD/EVB simulations to within the combined simulation and experimental errors.Particularly, the fact that the entropy terms become correct clearly supports the approximation above, where the change from a 1 to a 55 M standard state is considered solely as a translational entropy effect.That the free energies are exactly reproduced is obvious in view of the EVB parametrization, but the agreement with regard to the enthalpy and entropy components is quite remarkable since only their sum is subject to fitting.
To shed some light on the origin of the enthalpy and entropy terms, it is useful to explore the effect of modifying the EVB potential by varying H 12 and Δα.That is, by changing these parameters, but using exactly the same underlying MD trajectories, we can move the diabatic free energy curves and the transition state on the ground-state potential up or down, as shown in Figure 1.Starting from our standard EVB parametrization for the proton transfer reaction, with H 12 = 52.6 and Δα = 99.55 kcal/mol, we add or subtract 10 kcal/mol to each of these parameters and recalculate the Arrhenius plots for the forward reaction.The results of this exercise are summarized in Table 2, where it can be seen that ΔG ⧧ goes up when Δα is increased and when H 12 is decreased and down when these parameters are changed in the opposite directions.This thus follows the logic of Figure 1.However, what is interesting here is that it is ΔH ⧧ that closely follows the shifts in ΔG ⧧ while TΔS ⧧ is almost invariant with our modifications of the EVB potential.Hence, the TΔS ⧧ term must be intrinsically determined by the underlying force fields U 1 (R) and U 2 (R), OPLS-AA/M in our case. 17ince the TΔS ⧧ term is found to be invariant with respect to the EVB parametrization for our proton transfer reaction, this must mean that the geometry and charge distribution of the transition state region does not change appreciably when we modify H 12 and Δα as above.This would imply that the state coefficients, c 1 2 + c 2 2 = 1, of the ground-state eigenvector at the transition state do not change very much when we modify the EVB potential.These coefficients are given by ( ) and the average values of c 1 2 and c 2 2 for our standard parametrization are 0.45 and 0.55 at the transition state, respectively.As an example, when we change Δα to 89.55 kcal/mol (third entry in Table 2), these coefficients change only to 0.47 and 0.53.This means that the transition state region will largely be sampled by the same λ windows, and hence, its entropy relative to the reactant state is determined by the corresponding mixture of U 1 and U 2 .This can be further verified by simply calculating the Arrhenius plot corresponding to the free energy difference between ΔG(λ m ) at λ = 0.55 and λ = 0.10.This free energy difference thus has nothing to do with EVB but only represents the change of mapping potential from the reactant minimum in Figure 4b (λ = 0.10) to the mixture of U 1 and U 2 that roughly corresponds to the transition state (λ = 0.55).While the enthalpy change between these two mapping points is much larger than ΔH ⧧ , since there is no H 12 coupling involved in U map (λ) (Figure 1), it turns out that the corresponding TΔS(λ = 0.55) − TΔS(λ = 0.10) = −2.6 kcal/ mol, which is very close to the value of TΔS ⧧ obtained in Table 2. Hence, we can conclude that the entropy term at the transition state is largely determined simply by a mixture of the underlying state potentials.
The above conclusion can be further illustrated for our proton transfer reaction by comparing the TΔS term along the mapping potential (eq 4) with that on the ground state along the reaction coordinate.In the latter case, we therefore bin the value of c 2 2 along the reaction coordinate ΔU to get data points that more closely correspond to the λ values than the ΔU coordinate does (Figure 5a).For each such bin, we then calculate an Arrhenius plot from the five temperatures, both for ΔG(λ) and for ΔG g (c 2 2 ), and the result is shown in Figure 5b in terms of TΔS(λ), TΔS g (c 2 2), and the ground-state enthalpy term ΔH g (c 22 ).The curves in Figure 5b are plotted from the reactant minimum at λ = c 2 2 = 0.10 to the product minimum at λ = c 2 2 = 0.90.It can be seen that TΔS g (c 2 2 ) closely follows TΔS(λ) with a value at the transition state (λ = c 2 2 = 0.55) of −2.5 kcal/mol (at 300 K).Moreover, the entropy term decreases monotonically along the path from reactant to product for this reaction, while the corresponding enthalpy term ΔH g (c 2 2) reaches a maximum at the transition state and then decreases toward the product.
Simulations of Hydride Transfer in a Short Chain Dehydrogenase.In order to address the generality of the above findings, particularly the insensitivity of TΔS ⧧ to the EVB parametrization, we also examine the energetics of an enzyme reaction.Here, we will consider the reduction of 3oxovalerate to (R)-3-hydroxyvalerate by the psychrophilic hydroxybutyrate dehydrogenase enzyme (HBDH) from Psychrobacter arcticus.This reaction proceeds by concerted hydride and proton transfer to the substrate from the cofactor NADH and Tyr161, respectively (Figure 6a). 11,20It is an example with a relatively large experimentally determined entropy penalty for the tetrameric enzyme of TΔS ⧧ = −6.1 kcal/mol at 283 K. 20 We earlier modeled it via MD/EVB simulations using both the Q 21 and GROMACS 22 software for different oligomeric states of the enzyme. 5,23The tetramer calculations gave an entropy penalty of TΔS calcd ⧧ = −6.4kcal/ mol, in excellent agreement with experimental results, and interestingly, the thermodynamic activation parameters were found to be more or less insensitive to the oligomeric state. 5urther, the calculations using Q and GROMACS for the HBDH monomer gave virtually identical results (TΔS calcd ⧧ = −5.5 and TΔS calcd ⧧ = −5.2kcal/mol, respectively, at 283 K), 5,23 and here we will use the latter set of calculations as an example.
The EVB model for HBDH was parametrized previously on average minimum energy paths from 10 replicate QM/MM minimum energy paths with final energies calculated at the M062X/ma-def2-TZVPP level, including thermal corrections at 283 K. 11 These calculations turned out to yield very accurate results compared to experimental kinetic data, and the predicted values of the activation and reaction free energies were ΔG ⧧ = 15.9 and ΔG 0 = +3.9kcal/mol (ΔG expt ⧧ = 16 kcal/ mol). 11The EVB calibration from GROMACS simulations of the monomer reaction in a periodic dodecahedral box at 283 K gave parameter values of H 12 = 168.54and Δα = 96.70kcal/ mol and the resulting free energy profile is shown in Figure 6b. 23The corresponding Arrhenius plot obtained from these MD/EVB simulations at 273, 283, 293, 303, and 313 K is shown in Figure 6c, and it predicts values of ΔH ⧧ = 10.5 and TΔS ⧧ = kcal/mol at 283 K.
As above, we now examine the effect of changing the EVB parametrization by adding 10 kcal/mol to either H 12 or Δα.Addition of this energy quantity to Δα basically lifts the product state of the reaction, and as expected, we now get an increase of both ΔG ⧧ and ΔG 0 to 19.1 and 9.9 kcal/mol, respectively (Figure 6b).The activation enthalpy increases from ΔH ⧧ = 10.5 to ΔH ⧧ = 13.7 kcal/mol and, hence, perfectly matches the corresponding shift in ΔG ⧧ (Figure 6c).There is thus no significant change of the entropy penalty, TΔS ⧧ = −5.3kcal/mol (at 283 K).Likewise, if we add 10 kcal/ mol to H 12 , the main effect is that the free energy barrier goes down due to the increased coupling between reactant and product states, and we now get ΔG ⧧ = 14.0 and ΔG 0 = 4.1 kcal/mol (Figure 6b).Again, this shift is basically only enthalpic, and the activation parameters now become ΔH ⧧ = 8.9 and TΔS ⧧ = −5.0kcal/mol (Figure 6c).Hence, as with the case of our proton transfer reaction in solution, we can ) along the ground-state reaction path from the reactant to the product minimum, together with the entropy term along the mapping potential (eq 4).Note that the corresponding enthalpy on the mapping potential (not shown) is considerably higher since it does not involve any coupling between the states.Data for all three curves are plotted relative to the reactant minimum.conclude that the activation entropy penalty is not really determined by the EVB parameters, but by the degree of mixing of the two states in the transition state region.In this case, we also have c 1 2 and c 2 2 both ∼0.5 (c 1 2 = 0.48, c 2 2 = 0.52), and the entropy penalty is thus mainly determined by this particular mixing of the reactant and product force fields.As before, the entropy term on the ground state basically follows that on the mapping potential, and both give entropy penalties of TΔS ≈ − 5 kcal/mol at λ = c 2 2 = 0.52 (at 283 K).

■ DISCUSSION
Herein, we have addressed the question of why EVB simulations are able to produce such an accurate description of the temperature dependence of chemical in enzymes and solution in terms of computational Arrhenius plots.Obviously, if the activation free energy is correctly predicted at a given temperature, this just means that the sum of ΔH ⧧ and − TΔS ⧧ is correctly predicted.However, since these components of the free energy are never individually involved in the EVB parametrization, it appears somewhat surprising that they are also almost always correctly predicted.
We have shown here that the answer to this lies in the use of an accurate classical force field for the reactant and product states, in the case of a chemical reaction step that can be represented by a two-state EVB Hamiltonian.Moreover, for enzyme reactions that are parametrized on a reference reaction in solution and not on large DFT cluster or QM/MM enzyme models, the force field description must also capture the effect of changing the surroundings of the reacting groups from water to the enzyme environment.A case in point here is the recent MD/EVB study of chorismate mutase where, first, the uncatalyzed solution reaction at 298 K was parametrized on DFT calculations (in terms of ΔG ⧧ and ΔG 0 ) with a mixed explicit/implicit solvent representation. 7MD/EVB simulations of the temperature dependence then gave predictions of ΔH ⧧ and −TΔS ⧧ within 0.4 kcal/mol of the experimental values for the solution reaction.Upon transfer of the solution EVB model to two different enzyme variants, simulations of the catalyzed reactions also gave excellent agreement with the experimentally derived thermodynamic parameters and, particularly so, for the TΔS ⧧ terms that were within 0.7 kcal/mol of the measured values. 7ur analysis herein shows that the entropy term is basically determined by the mixture of U 1 and U 2 that defines the transition state region.This mixture is typically near 50/50 proportions for reactions that are not extremely endothermic or exothermic (for an example of more complex three-state EVB simulations, see ref 24).The acetate catalyzed abstraction of a proton from androstenedione is an example of this.The entropy change on the mapping potential U map (λ) between the λ values corresponding to the transition state and reactant minimum is thus almost the same as the entropy change on the EVB ground-state potential.Even for more complex chemical reactions that involve more bonds being simultaneously broken and formed, the correspondence between the entropy on the two potentials is close, and TΔS g remains essentially insensitive to the EVB parametrization.This is illustrated by the concerted HBDH reaction above, which involves four bonds being simultaneously broken or formed.Hence, we can conclude that an accurate force field, itself, will ensure a good estimate of the activation entropy penalty, irrespective of how the EVB model is parametrized.This means that if the EVB model gives an accurate prediction of ΔG ⧧ , then it will also give a good estimate of ΔH ⧧ (since ΔG ⧧ = ΔH ⧧ − TΔS ⧧ ).
It may further be useful to point out here that, besides the H 12 and Δα parameters, the EVB model also offers the possibility of more detailed geometric parametrization of the transition state.That is, a softer exponential repulsion term (U rep = Ae −βr ij ) instead of the Lennard-Jones 1/r 12 is usually used between atoms involved in bond formation and breaking, and this repulsion can be tuned to obtain desired distances. 23oreover, coupling terms between the Morse potentials describing changes in bonding are also usually introduced with respect to bond angles, torsions, and improper dihedrals, so that such terms vanish as bonds dissociate. 25Taken together, the parametrization of ΔG ⧧ in terms of H 12 and Δα from suitable target data, possibly with additional fine-tuning of the transition state geometry, makes the EVB model very powerful for describing chemical reactions both in solution and in enzymes.What makes the method capable of correctly also describing the temperature dependence of reactions turns out to be the fact that the underlying force fields are accurate enough to properly estimate the entropy changes along the reaction path.

Simulation of Charge Transfer between Na and Na + .
The Na atom and Na + ion, described by standard parameters, 26 were enclosed in a TIP3P water sphere of diameter 40 Å and kept fixed at a 6 Å separation, by restraining each ion to its initial position (by 100 kcal/mol/Å 2 ).For this system, 10 replicate MD/FEP simulations (Δt = 1 fs) were carried out at each temperature (280, 290, 300, 310, and 320 K), starting from different initial velocities.Each replica involved 1.25 ns of heating and equilibration at the target temperature, followed by a 1.05 ns FEP calculation where the charge was moved from Na + to Na in 21 discrete λ steps (eq 4).
Acetate Catalyzed Deprotonation of 5-Andro-3,17dione.The reactants and products were described by the OPLS-AA/M force field, 17 and parameters were generated with Schrodinger's ffld_server. 27The reacting molecules (acetate and steroid) were enclosed in a water sphere of diameter 40 Å, and the steroid position and orientation were weakly restrained by applying harmonic positional restraints (5 kcal/mol/Å 2 ) to three carbon atoms of the ring system.To account for the 55 M standard state, a 10 kcal/mol/Å 2 distance restraint was applied between the acetate acceptor oxygen and the steroid donor carbon, which keeps these atoms in contact during the simulations.A soft exponential repulsion term (U rep = Ae −βr ij ) was applied between the donor and acceptor atoms with A = 12 100 kcal/mol and β = 2.5 Å −1 .At each temperature (290, 295, 300, 305, and 310 K), 60 replicate MD/EVB free energy calculations were carried out with a 1 fs time step.Each of these replicas involved 200 ps of initial equilibration, followed by 1.05 ns of data collection during the FEP transformation, which involved 21 discrete λ windows.This gives error bars (s.e.m.) for the calculated activation and reaction free energies of less than 0.1 kcal/mol at all temperatures.The MD/EVB simulations (also for the Na + Na + system) were performed with the Q software package, 21 where water molecules at the sphere boundary were subjected to radial and polarization restraints following the SCAAS model. 21,28Water−water Lennard-Jones interactions were truncated beyond 10 Å and long-range electrostatic inter-

Journal of Chemical Theory and Computation
actions beyond this cutoff were treated with the local reaction field multipole expansion method. 29imulation the HBDH Reaction.The HBDH calculations presented are taken from ref 23, where details of the simulations are given.Briefly, they were performed using the GROMACS 2022 software 22 with the OPLS-AA/L force field 30 and utilized the crystallographic structure with PDB code 6ZZP 11 and the starting point.The protonation states of the ionizable residues were calculated with PROPKA, 31 and the system was enclosed in a dodecahedral box with its nearest wall at 20 Å from the edge of the protein.The box encapsulates 23 710 water molecules and nine Na + ions, for a total number of 75 191 atoms.The topologies for the MD/EVB simulations were generated using the gmxtools presented in ref 23 since, for each FEP window, a separate topology corresponding to the appropriate λ needs to be provided. 23The system was initially equilibrated in the reactant configuration in the NPT ensemble for 70 ns at 283 K. Following an additional short equilibration at the different target temperatures, 20 replicate FEP calculations were carried out with 51 λ windows for a total simulation time of 20.04 ns at each temperature.The Arrhenius plots were obtained from simulations at 273, 283, 293, 303, and 313 K.During the FEP simulations, coordinates were recorded every 10 time steps and were later used to recalculate the energies corresponding to the main diagonal of the EVB matrix.This was accomplished using the rerun option of the mdrun tool of GROMACS.The energies for each of the two valence states were then collected with the energy tool of GROMACS and analyzed with a modified version 23 of the Qfep tool of the Q software.

Figure 1 .
Figure 1.Effect of the EVB parameters H 12 and Δα on the free energy profiles for transfer of a negative charge between Na and Na + in water at a 6 Å separation.The curves with purple squares denote the mapping free energy ΔG(λ m ).Red and blue circles are the diabatic (ΔG 1 and ΔG 2 ) free energy curves, and the solid black curves denote the ground-state free energy (ΔG g ).The following sets of H 12 and Δα parameters are used: (a) 0, 0; (b) 20, 0; (c) 0, 20; and (d) 20, 20.The quantities ΔG 1 , ΔG 2 , and ΔG g are functions of energy gap ΔU.

Figure 2 .
Figure 2. Illustration of how eqs 6 and 7 work in practice.Each value of λ (color coded form green to blue) contributes to several bins of the discretized energy gap ΔU (denoted by solid lines).The resulting free energy value of each bin is obtained by statistical weighting of the contributions from each λ.

Figure 3 .
Figure 3. Calculated Arrhenius plots of (a) ΔG ⧧ /T vs 1/T and (b) ΔG ⧧ vs T for the transfer of a negative charge between Na and Na + in water at 6 Å separation with H 12 = 0.

Figure 4 .
Figure 4. (a) Mechanism of the acetate catalyzed proton abstraction from the 5-androstene-3,17-dione steroid.(b) Free energy profile of the reaction obtained from MD/EVB simulations.(c) Calculated Arrhenius plots from the forward (k 1 ) and backward (k −1 ) reactions and for the equilibrium constant (K eq ).(d) MD snapshot from the transition state region of the reaction.

Figure 5 .
Figure 5. (a) Relationship between the energy gap and the EVB state coefficients corresponding to the reactant and product weights obtained from the MD/EVB simulations at 300 K. (b) Calculated thermodynamic parameters (ΔH g and TΔS g ) along the ground-state reaction path from the reactant to the product minimum, together with the entropy term along the mapping potential (eq 4).Note that the corresponding enthalpy on the mapping potential (not shown) is considerably higher since it does not involve any coupling between the states.Data for all three curves are plotted relative to the reactant minimum.

Figure 6 .
Figure 6.(a) View of the EVB transition state for concerted hydride and proton transfer in HBDH, where a typical MD snapshot at the top of the free energy barrier is shown.Calculated free energy profiles for the reaction with the standard H 12 and Δα parameters (168.54 and 96.70 kcal/ mol−black circles) 23 and with 10 kcal/mol added to either H 12 (red circles) or Δα (purple circles).(c) The corresponding Arrhenius plots from five different temperatures for the three cases of EVB parameters (color coding as in panel b).

Table 2 .
Effect of Modifying H 12 and Δα on the Activation Parameters for the Forward Proton Transfer Reaction at 300 K (kcal/mol) a aThe standard EVB parametrization corresponds to the first row in the table.