Theoretical and Computational Chemistry

---

• Proton and hydride transfer reactions in enzymes
• Proton-coupled electron transfer reactions
• Mixed quantum/classical molecular dynamics methods
• Multistate continuum theory
• Nuclear-electronic orbital (NEO) method
  
  

---

Proton and hydride transfer reactions in enzymes

 

We have developed a hybrid approach for simulating proton and hydride transfer reactions in enzymes.^41,46,48  This hybrid approach includes electronic and nuclear quantum effects, as well as the motion of the entire solvated enzyme.  The methodology provides detailed mechanistic information at the molecular level and allows the calculation of rates and kinetic isotope effects.  This hybrid approach also enables us to investigate the relation between enzyme motion and activity.⁵⁵

 

1. Liver alcohol dehydrogenase (LADH)^34,37,41,48 

a)     We performed electronic structure calculations at various levels of theory on a 148-atom model of the active site and classical molecular dynamics simulations on the entire solvated LADH dimer.³⁴  These calculations support the hypothesis that alcohol deprotonation occurs prior to the hydride transfer step and that the alcohol deprotonation facilitates the hydride transfer by lowering the barrier for hydride transfer.  The calculations indicate that the electrostatic interaction between the substrate alkoxide oxygen and the zinc counterion in the active site lowers the barrier to hydride transfer.

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

b)     We studied the nuclear quantum effects in LADH by combining electronic structure methods with the calculation of hydrogen vibrational wavefunctions along various reaction paths.³⁷  The results indicate that the hydride transfer is adiabatic and hydrogen tunneling does not play a critical role along the minimum energy path.  In contrast, nonadiabatic effects and hydrogen tunneling were shown to be important along the more relevant straight-line reaction paths.  The secondary hydrogens were found to be significantly coupled to the transferring hydride near the transition state.  Moreover, hydrogen tunneling in LADH was strongly impacted by the puckering of the NAD+ ring and the distance between the donor and acceptor carbons.

c)     We applied the hybrid quantum/classical molecular dynamics approach to hydride transfer in LADH.^41,48  The transmission coefficient was calculated to be nearly unity, implying that dynamical barrier recrossings are not dominant for this reaction.  The calculated deuterium and tritium kinetic isotope effects for the overall rate were in agreement with experimental results.  The simulations provided evidence of hydrogen tunneling in the direction along the donor-acceptor axis.  An analysis of the geometrical parameters during the equilibrium and nonequilibrium simulations elucidated the relation between specific enzyme motions and enzyme activity.  The simulations indicated that the experimentally observed effect of mutating Val-203 on the enzyme activity is due to the alteration of the equilibrium free energy barrier rather than nonequilibrium dynamical factors.  The promoting motion of Val-203 was characterized in terms of steric interactions involving Thr-178 and the coenzyme.

 

 

 

 

 

 

 

 

 

 

 

 

  LADH movies:

LADH-productive
LADH-recrossing
LADH-unproductive

  In order see the movies you must enable your browser to play AVI movies.

2. Dihydrofolate reductase (DHFR)^50,51,56  

a)     We applied the hybrid approach to hydride transfer in DHFR.⁵¹  The calculated deuterium kinetic isotope effect agreed with the experimental value.  The simulations provided evidence of hydrogen tunneling in the direction along the donor-acceptor axis.  The transmission coefficient was calculated to be 0.8, indicating dynamical barrier recrossings.  Nonadiabatic transitions among the vibrational states were observed but did not strongly affect the transmission coefficient.

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

b)     In collaboration with the Benkovic group, we investigated the relation of motion and activity with three approaches: genomic analysis for sequence conservation, kinetics of multiple mutations, and analysis of mixed quantum/classical molecular dynamics simulations of hydride transfer in DHFR.⁵⁰  Our studies suggested the existence of a network of coupled promoting motions in DHFR, where promoting motions are defined as systematic changes in thermally averaged equilibrium properties as the reaction evolves along the collective reaction coordinate.  These promoting motions represent the conformational changes that occur during the reaction and correspond to the reorganization of the environment for the charge transfer process.  Such equilibrium motions along the collective reaction coordinate influence the activation free energy barrier and occur on the timescale of the catalyzed chemical reaction (i.e., the millisecond timescale for DHFR).  The coupled promoting motions were found to extend throughout the DHFR enzyme.

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

c)     We also studied a mutant DHFR enzyme.⁵⁶  Although residue 121 is on the exterior of the enzyme and is more than 10 Angstroms from the active site, experimental studies showed that the mutation of Gly-121 to valine reduces the rate of hydride transfer by a factor of 163.  Our simulations suggested that this mutation interrupts the network of coupled promoting motions and therefore increases the free energy barrier by an amount consistent with the experimentally observed rate reduction.  The calculated transmission coefficients were found to be comparable for the wild type and mutant enzymes.

  DHFR movies:

DHFR-productive
DHFR-recrossing
DHFR-unproductive

  In order see the movies you must enable your browser to play AVI movies.

---

Proton-coupled electron transfer reactions

 

Proton-coupled electron transfer (PCET) reactions play a critical role in a variety of chemical and biological processes.  We have developed a multistate continuum theory for proton-coupled electron transfer (PCET) that includes both electronic and nuclear quantum effects.^30,35,40,43  In this theory, the solute is represented by a multistate valence bond model, the solvent is described by a dielectric continuum, and the transferring hydrogen nuclei are represented as quantum mechanical wavefunctions.  A PCET reaction involving the transfer of one electron and one proton is described in terms of four diabatic states. The free energy surfaces are calculated as functions of two collective solvent coordinates corresponding to proton transfer (PT) and electron transfer (ET), respectively.  Rate expressions for PCET were derived in various limits.³⁵  Calculation of the rates requires the gas phase valence bond matrix elements and the reorganization energies.  The gas phase valence bond matrix elements are represented by molecular mechanical terms fit to electronic structure calculations or experimental data.  The inner-sphere (solute) reorganization energy matrix elements may be calculated from the equilibrium force constants and bond lengths.  The outer-sphere (solvent) reorganization energy matrix elements are calculated with an electrostatic dielectric continuum model.  This theory was applied to a series of model PCET systems to predict the dependence of the rates, mechanisms, and kinetic isotope effects (KIEs) on the solute and solvent properties.³⁹  Our applications to chemically and biologically important systems are described below.^{33,44,45,52,53}  We have also developed the methodology for mixed quantum/classical molecular dynamics simulations with explicit solvent for PCET reactions.^{15,18,19,27,47} 

 

1. Iron bi-imidazoline complexes⁴⁴

 

 

 

 

 

 

 

 

 

 

 

 

 

A comparative experimental study of single ET and PCET reactions in iron bi-imidazoline complexes indicated that the rates of ET and PCET are similar.  Previously this result was explained in the context of adiabatic Marcus theory, and the PCET reaction was viewed as a hydrogen atom transfer involving negligible solute charge rearrangement, leading to zero solvent reorganization energy.  The similarity of the ET and PCET rates was thought to be due to the compensation of the larger solvent reorganization energy for ET by a larger solute reorganization energy for PCET.  The kinetic isotope effect (the ratio of the rate for hydrogen to the rate for deuterium) for PCET was measured to be a moderate value of 2.3.  Our calculations, which were based on nonadiabatic rate expressions for ET and PCET, provided an alternative explanation for the experimental results. In our calculations, the inner-sphere reorganization involving the Fe-N bonds was assumed to be the same for both ET and PCET.  The solvent reorganization energy for PCET was found to be substantial and was »1-3 kcal/mol lower than the solvent reorganization energy for ET.  The overall coupling for PCET was found to be smaller than the coupling for ET due to averaging over the reactant and product proton vibrational wavefunctions (i.e., multiplying by the vibrational overlap factor).  The calculations indicated that the similarity of the rates for ET and PCET was due mainly to the compensation of the larger solvent reorganization energy for ET by the smaller coupling for PCET.  The moderate KIE was determined to arise from the relatively large overlap factor and the significant contributions from excited vibronic states. 

 

2. Ruthenium polypyridyl complexes⁵²

 

 

 

 

 

 

 

 

 

 

 

 

 

An experimental study of PCET in ruthenium polypyridyl complexes revealed that the ComB rate is nearly one order of magnitude larger than the CompA rate, and the CompA KIE of 16.1 is larger than the CompB KIE of 11.4.  Density functional theory calculations illustrated that the steric crowding near the oxygen proton acceptor is significantly greater for CompA than for CompB.  Consistent with this observation, our calculations implied that the proton donor-acceptor distance is larger for CompA than for CompB, leading to a larger overlap between the reactant and product proton vibrational wavefunctions for CompB than for CompA.  The rate for CompB is larger than the rate for CompA because the rate increases as this overlap factor increases.  The KIE for CompB is smaller than the KIE for CompA because the KIE decreases as this overlap factor increases.  Both of these KIEs are greater than the KIE for the iron bi-imidazoline complexes because the vibrational overlap factor is smaller for the ruthenium systems.

 

3. Osmium complex to benzoquinone

 

 

 

 

 

 

 

 

 

Experimental studies of PCET reactions from a series of osmium complexes to benzoquinone identified unusually high KIE’s of up to 408.  Our theoretical calculations illustrated that these colossal KIEs arise from the relatively small overlap between the reactant and product hydrogen vibrational wavefunctions.    The KIE increases as the vibrational overlap decreases and as the contribution of transitions between the lowest energy reactant and product vibronic states increases. The trends in the KIEs for a series of complexes were found to be determined by a balance among several factors, including the X-H frequencies and proton transfer distances for the different proton donors (X=N, P, S), as well as the solvent reorganization energies and reaction free energies for the different complexes.  These characteristics of the osmium systems influence the overlaps between the reactant and product hydrogen vibrational wavefunctions and the relative contributions of the excited vibronic states, which in turn impact the KIE.

 

 

4. Amidinium-carboxylate salt bridges^33,45

 

 

 

 

 

 

 

 

 

In PCET through amidinium-carboxylate salt bridges, the ET reaction is coupled to the motion of two protons at the proton transfer interface.  In this case, the reaction was described in terms of eight valence bond states to include all possible charge transfer states, two hydrogen nuclei were treated quantum mechanically, and the free energy surfaces depended on three solvent coordinates corresponding to the electron and two proton transfer reactions.  Experimental studies of photoinduced PCET in analogous systems revealed that the rate for the donor-(amidinium-carboxylate)-acceptor system is substantially slower than the rate for the switched interface donor-(carboxylate-amidinium)-acceptor system.  The calculations illustrated that this difference in rates is due mainly to the opposite dipole moments at the proton transfer interfaces for the two systems, leading to an endothermic reaction for the donor-(amidinium-carboxylate)-acceptor system and an exothermic reaction for the switched interface system. 

 

 

5. DNA-acrylamide complexes⁵³

 

 

 

 

 

 

 

 

 

Experiments imply that PCET may occur in DNA-acrylamide complexes.  The influence of neighboring DNA base pairs was determined theoretically by studying both the solvated thymine-acrylamide complex and solvated DNA-acrylamide models.  The calculations indicated that the final product corresponds to single ET for the solvated thymine-acrylamide complex but to a net PCET reaction for the solvated DNA-acrylamide complex.  This difference is due to a decrease in solvent accessibility in the presence of DNA, which alters the relative free energies of the ET and PCET product states.  Thus, the balance between ET and PCET in the DNA-acrylamide system is highly sensitive to the solvation properties of the system.

 

6. Tyrosine oxidation in Photosystem II model

 

 

 

 

 

 

 

 

 

 

 

 

The pH- and temperature-dependence of the rates for a model representing tyrosine oxidation in Photosystem II were measured experimentally.  The mechanism was determined to be PCET at pH<10 when the tyrosine is initially protonated and single ET for pH>10 when the tyrosine is initially deprotonated.  The PCET rate was found to increase monotonically with pH, whereas the single ET rate was found to be independent of pH and two orders of magnitude faster than the PCET rate.  Initially the slower rate for PCET was attributed to a larger reorganization energy.  Our calculations reproduced the experimentally observed trends and provided an alternative explanation for the difference in rates.  The pH dependence for the PCET reaction resulted from the decrease in the reaction free energies with pH.  The calculated solvent reorganization energy for the overall PCET reaction was only slightly larger than that for single ET.  The calculations indicated that the larger rate for single ET arises from the averaging of the coupling for PCET over the reactant and product hydrogen vibrational wavefunctions (i.e., the vibrational overlap factor).

---

Mixed quantum/classical molecular dynamics methods

 

Proton, hydride, and proton-coupled electron transfer reactions play an essential role in a wide variety of chemical and biological processes.  Kinetic isotope effect experiments indicate that nuclear quantum effects such as hydrogen tunneling are significant in these types of reactions.  We have been developing mixed quantum/classical molecular dynamics methodology for the simulation of proton, hydride, and proton-coupled electron transfer reactions in solution and proteins.^24,26

 

1. A hybrid approach was developed for including electronic and nuclear quantum effects in molecular dynamics simulations of hydrogen transfer reactions in enzymes.^41,46,48  The electronic quantum effects are incorporated with an empirical valence bond potential.  The nuclear quantum effects of the transferring hydrogen are included with a mixed quantum/classical molecular dynamics method in which the hydrogen nucleus is represented as a multidimensional vibrational wavefunction.  The free energy profiles are obtained as functions of a collective reaction coordinate.  A perturbation formula was derived to incorporate the vibrationally adiabatic nuclear quantum effects into the free energy profiles.  The dynamical effects are studied with the molecular dynamics with quantum transitions (MDQT) surface hopping method,^10-12,26 which incorporates nonadiabatic transitions among the adiabatic hydrogen vibrational states.  The MDQT method is combined with a reactive flux approach to calculate the transmission coefficient and to investigate the real-time dynamics of reactive trajectories.  This hybrid approach has been applied to several enzyme reactions.^{41,48,50,51,56}
2. The MDQT method was compared to exact quantum dynamical and mean field calculations for model single proton, double proton, and proton-coupled electron transfer reactions.^14,15,19,29 
3. An improved treatment of classically forbidden nonadiabatic transitions was implemented.^29,32  
4. The Fourier grid Hamiltonian multiconfigurational self-consistent-field (FGH-MCSCF) and multidimensional grid generation methods^38,42 were developed for the calculation of multidimensional hydrogen vibrational wavefunctions.  These methods were utilized in applications to hydride transfer in enzymes.^{41,48,50,51,56}
5. The multiconfigurational MDQT (MC-MDQT) method was developed for the simulation of multiple proton transfer reactions, where several hydrogen nuclei must be treated quantum mechanically.^13,17  The MC-MDQT method was applied to proton wires (i.e., proton transport along hydrogen-bonded water chains).^{17,20,21,28,31}  

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

6. The MDQT method was extended to treat the donor-acceptor vibrational motion as well as the hydrogen motion quantum mechanically for hydrogen transfer reactions.  This approach was applied to a model representing intramolecular proton transfer within a phenol-amine complex in liquid methyl chloride. 
7. Mixed quantum/classical molecular dynamics methods were developed for simulating proton-coupled electron transfer reactions.^{15,18,19,27,47}

 

---

Multistate continuum theory

 

We have developed a multistate continuum theory for multiple charge transfer reactions such as proton-coupled electron transfer and multiple proton transfer reactions.^30,40,43  The solute is represented by a multistate valence bond model, the solvent is described by a dielectric continuum, and the transferring hydrogen nuclei are represented by quantum mechanical wavefunctions.  This theory provides free energy surfaces that depend on a set of scalar solvent coordinates corresponding to the individual charge transfer reactions.  For processes involving significant inner-sphere (i.e., solute) reorganization, the effects of solute intramolecular vibrations can be included.  Rate expressions for proton-coupled electron transfer reactions were derived in various limits.³⁵  Calculation of the rates requires the gas phase valence bond matrix elements and the reorganization energies.  The gas phase valence bond matrix elements are represented by molecular mechanical terms fit to electronic structure calculations or experimental data.  The inner-sphere (solute) reorganization energy matrix elements can be calculated from the equilibrium force constants and bond lengths.  The outer-sphere (solvent) reorganization energy matrix elements are calculated with an electrostatic dielectric continuum model.  Currently this theory is being extended to include the proton donor-acceptor vibrational motion.  This theoretical formulation has been applied to a wide range of proton-coupled electron transfer reactions.^{33,39,44,45,52,53}

 

 

 

 

 

 

 

 

 

 

 

 

 

 

---

Nuclear-electronic orbital (NEO) method

 

We have developed the nuclear-electronic orbital (NEO) method for the incorporation of nuclear quantum effects in electronic structure calculations.^54,57  In the NEO approach, specified nuclei are treated quantum mechanically on the same level as the electrons, and mixed nuclear-electronic wavefunctions are calculated variationally with molecular orbital methods.  Both electronic and nuclear molecular orbitals are expressed as linear combinations of Gaussian basis functions.  The variational method is utilized to minimize the energy with respect to all molecular orbitals, as well as to optimize the centers of the nuclear basis functions.  Significant correlation effects are included using a multiconfigurational self-consistent-field (MCSCF) approach.  Analytic gradients allow for geometry optimizations and for the generation of minimum energy paths and dynamic reaction paths.  For hydrogen transfer reactions, the transferring hydrogen nuclei, as well as all electrons, are treated quantum mechanically to include nuclear quantum effects such as zero point energy and hydrogen tunneling.  In this case, the reaction coordinate depends explicitly on only the classical nuclei, and the imaginary mode at the transition state corresponds to heavy-atom motion that drives the charge transfer reaction via reorganization of the environment.  The advantages of the NEO approach are that nuclear quantum effects are incorporated during the electronic structure calculation, the Born-Oppenheimer separation of electrons and nuclei is avoided, excited vibrational-electronic states may be calculated, and its accuracy may be improved systematically.

 

1. We derived and implemented the NEO-HF (Hartree-Fock), NEO-CI (configuration interaction), and NEO-MCSCF (multiconfigurational self-consistent-field) methods.⁵⁴  The NEO approach was tested with applications to small molecules.  The bond lengths and vibrational energy splittings for the H₂ and HF molecules are in agreement with experimental results.  The tunnel splitting for malonaldehyde is also reasonable.
2. We have developed the methodology for a vibrational analysis within the NEO framework.⁵⁷  This approach has been applied to a series of molecules, including HCN, the protonated water dimer, triazene, and a symmetric S_N2 reaction.

 

 

 

 

 

 

 

 

3. Currently the NEO method is being applied to proton transfer in the condensation steps required for the synthesis of POSS (polyhedral silsesquioxanes) and to proton transfer in ionic liquids.