Empirical Valence Bond Method Research Articles

Toward transferable empirical valence bonds: Making classical force fields reactive.

The empirical valence bond technique allows classical force fields to model reactive processes. However, parametrization from experimental data or quantum mechanical calculations is required for each reaction present in the simulation. We show that the parameters present in the empirical valence bond method can be predicted using a neural network model and the SMILES strings describing a reaction. This removes the need for quantum calculations in the parametrization of the empirical valence bond technique. In doing so, we have taken the first steps toward defining a new procedure for enabling reactive atomistic simulations. This procedure would allow researchers to use existing classical force fields for reactive simulations, without performing additional quantum mechanical calculations.

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.

Examination of the performance of degraded Nafion membrane with graphene oxide

The degradation of Nafion membrane easily occurs after a long run of proton exchange membrane fuel cell. Graphene-oxide (GO) is regarded as an ideal material to form composite membrane. In this work, a multistate empirical valence bond (MS-EVB) method and molecular dynamics method are employed to investigate the effect of GO on transport and solvation characteristics of protons in the degraded hydration membrane. The hydrophilic cluster and its connectivity are analyzed. The interaction between GO, Nafion membrane and water molecules is discussed. The influence of GO on the thermal conduction of Nafion membrane is also evaluated. The results demonstrate that with the addition of GO, the water cluster connectivity of Nafion membrane becomes better. The distribution of hydronium ions and membranes in the water is more uniform, which increases the proton conductivity. The larger GO may not interact well with water molecules, but promotes the thermal conductivity of degraded Nafion membrane more obviously.

Effect of Environmental Factors on the Catalytic Activity of Intramembrane Serine Protease.

The cleavage of protein inside cell membranes regulates pathological pathways and is a subject of major interest. Thus, the nature of the coupling between the physical environment and the function of such proteins has recently attracted significant experimental and theoretical efforts. However, it is difficult to determine the nature of this coupling uniquely by experimental and theoretical studies unless one can separate the chemical and the environmental factors. This work describes calculations of the activation barriers of the intramembrane rhomboid protease in neutral and charged lipid bilayers and in detergent micelle, trying to explore the environmental effect. The calculations of the chemical barrier are done using the empirical valence bond (EVB) method. Additionally, the renormalization method captures the energetics and dynamical effects of the conformational change. The simulations indicate that the physical environment around the rhomboid protease is not a major factor in changing the chemical catalysis and that the conformational and substrate dynamics do not exhibit long-time coupling. General issues about the action of membrane-embedded enzymes are also considered.

In silico simulation of the excited state proton transfer reaction of 2-(2-furyl)-3-hydroxychromone (FHC) in solution by empirical valence bond (EVB) method in conjunction with classical molecular dynamics

2-(2-Furyl)-3-hydroxychromone (FHC) is a dual emitter fluorescent probe exhibiting excited state intramolecular proton transfer (ESIPT) reaction in solution and inserted in biological relevant macromolecules. Its two fluorescence bands are attributed to normal form (N*) for the short wavelength and the long-wavelength band to the tautomer form (T*), the product of the ESIPT reaction. To rationalise the mechanisms that lead to the occurrence of the ESIPT reaction and predict the fluorescence spectra, the excited and ground state dynamic of this compound was followed in acetonitrile, ethanol, methanol, N-methylformamide and formamide, by the classical molecular dynamic simulation in conjunction with the empirical valence bond (EVB) model to describe bond dissociation process. A set of transferable parameters were used to fit the EVB potential energy curve to TDDFT data. Using these parameters, the occurrence time of the ESIPT reaction was predicted to be 7, 84, 113, 98, 207ps, respectively in the considered solvents, in excellent agreement with the experimentally reported values. These values are comparatively higher than those determined for the parent flavone dye (less than 100 fs), giving FHC a higher potential for sensor development. By examining structural parameters such as the distance from the proton to acceptor, the intermolecular H-bond between FHC and the solvent molecules and proton trajectories, we found that ESIPT occured in acetonitrile and alcohols as a combination of low frequency high amplitude donor-acceptor bending motion and charge redistribution in the excited state, that were modulated by the solvent H-bond donor potency. Differently in amides ESIPT occured via tunnelling mechanism due to the close energy of N* and T*. Fluorescence spectra calculated from potential energies of the EVB model along the reaction coordinates were qualitatively in agreement with the recorded spectra.

Exploring the Mechanism of Covalent Inhibition: Simulating the Binding Free Energy of α-Ketoamide Inhibitors of the Main Protease of SARS-CoV-2.

The development of reliable ways of predicting the binding free energies of covalent inhibitors is a challenge for computer-aided drug design. Such development is important, for example, in the fight against the SARS-CoV-2 virus, in which covalent inhibitors can provide a promising tool for blocking Mpro, the main protease of the virus. This work develops a reliable and practical protocol for evaluating the binding free energy of covalent inhibitors. Our protocol presents a major advance over other approaches that do not consider the chemical contribution of the binding free energy. Our strategy combines the empirical valence bond method for evaluating the reaction energy profile and the PDLD/S-LRA/β method for evaluating the noncovalent part of the binding process. This protocol has been used in the calculations of the binding free energy of an α-ketoamide inhibitor of Mpro. Encouragingly, our approach reproduces the observed binding free energy. Our study of covalent inhibitors of cysteine proteases indicates that in the choice of an effective warhead it is crucial to focus on the exothermicity of the point on the free energy surface of a peptide cleavage that connects the acylation and deacylation steps. Overall, we believe that our approach should provide a powerful and effective method for in silico design of covalent drugs.

Molecular Origins of the Barriers to Proton Transport in Acidic Aqueous Solutions.

The self-consistent iterative multistate empirical valence bond (SCI-MS-EVB) method is used to analyze the structure, thermodynamics, and dynamics of hydrochloric acid solutions. The reorientation time scales of irreversible proton transport are elucidated by simulating 0.43, 0.85, 1.68, and 3.26 M HCl solutions at 270, 285, 300, 315, and 330 K. The results indicate increased counterion pairing with increasing concentration, which manifests itself via a reduced hydronium oxygen-chloride (O*-Cl) structuring in the radial distribution functions. Increasing ionic concentration also reduces the diffusion of the hydrated excess protons, principally by reducing the contribution of the Grotthuss proton hopping (shuttling) mechanism to the overall diffusion process. In agreement with prior experimental findings, a decrease in the activation energy of reorientation time scales was also observed, which is explicitly explained by using activated rate theory and an energy-entropy decomposition of the state-averaged radial distribution functions. These results provide atomistic verification of suggestions from recent two-dimensional infrared spectroscopy experiments that chloride anions (as opposed to hydrated excess protons) create entropic barriers to proton transport.

Reactive Molecular Dynamics at Constant Pressure via Nonreactive Force Fields: Extending the Empirical Valence Bond Method to the Isothermal-Isobaric Ensemble

The Empirical Valence Bond (EVB) method offers a suitable framework to obtain reactive potentials through the coupling of nonreactive force fields. In this formalism, most of the implemented coupling terms are built using functional forms that depend on spatial coordinates, while parameters are fitted against reference data to model the change of chemistry between the participating nonreactive states. In this work, we demonstrate that the use of such coupling terms precludes the computation of the stress tensor for condensed phase systems and prevents the possibility to carry out EVB molecular dynamics in the isothermal-isobaric (NPT) ensemble. Alternatively, we make use of coupling terms that depend on the energy gaps, defined as the energy differences between the participating nonreactive force fields, and derive a general expression for the EVB stress tensor suitable for computation. Implementation of this new methodology is tested for a model of a single reactive malonaldehyde solvated in nonreactive water. Mass densities and probability distributions for the values of the energy gaps computed in the NPT ensemble reveal a negligible role of the reactive potential in the limit of low concentrated solutions, thus corroborating for the first time the validity of approximations based on the canonical NVT ensemble, customarily adopted for EVB simulations. The presented formalism also aims to contribute to future implementations and extensions of the EVB method to research the limit of highly concentrated solutions.

Oxide Relay: An Efficient Mechanism for Catalytic Water Oxidation at Hydrophobic Electrode Surfaces.

In order to combine the advantages of molecular catalysts with the stability of solid-state catalysts, hybrid systems with catalysts immobilized on carbon nanotubes are prominent candidates. Here we explore our recent mechanistic proposal for Ru(tda)(py)2, the oxide relay mechanism, in a hybrid system from an experimental study. It reacts with the same efficiency but with increased stability compared to the homogeneous molecular catalyst. We used the empirical valence bond method and molecular dynamics with enhanced sampling approaches to investigate the two key steps in the mechanism: the intramolecular O-O bond formation and the OH- nucleophilic attack. The results on these calculations show that the oxide relay mechanism remains unaltered in the new environment. We believe that the principles should apply to other oxide containing dangling groups and to other metal centers, opening new possibilities of future developments on hybrid molecular catalyst-based water splitting devices.

Exploring the Proteolysis Mechanism of the Proteasomes.

The proteasome is a key protease in the eukaryotic cells which is responsible for various important cellular processes such as the control of the cell cycle, immune responses, protein homeostasis, inflammation, apoptosis, and the response to proteotoxic stress. Acting as a major molecular machine for protein degradation, proteasome first identifies damaged or obsolete regulatory proteins by attaching ubiquitin chains and subsequently utilizes conserved pore loops of the heterohexameric ring of AAA+ (ATPases associated with diverse cellular activities) to pull and mechanically unfold and translocate the misfolded protein to the active site for proteolysis. A detailed knowledge of the reaction mechanism for this proteasomal proteolysis is of central importance, both for fundamental understanding and for drug discovery. The present study investigates the mechanism of the proteolysis by the proteasome with full consideration of the protein's flexibility and its impact on the reaction free energy. Major attention is paid to the role of the protein electrostatics in determining the activation barriers. The reaction mechanism is studied by considering a small artificial fluorogenic peptide substrate (Suc-LLVY-AMC) and evaluating the activation barriers and reaction free energies for the acylation and deacylation steps, by using the empirical valence bond method. Our results shed light on the proteolysis mechanism and thus should be important for further studies of the proteasome action.

Evolution Conserves the Network of Coupled Residues in Dihydrofolate Reductase.

Understanding protein motions and their role in enzymatic reactions is an important and timely topic in enzymology. Protein motions that are involved in the chemical step of catalysis are particularly intriguing but difficult to identify. A global network of coupled residues in Escherichia coli dihydrofolate reductase (E. coli DHFR), which assists in catalyzing the chemical step, has previously been demonstrated through quantum mechanical/molecular mechanical and molecular dynamics simulations as well as bioinformatic analyses. A few specific residues (M42, G121, F125, and I14) were shown to function synergistically with measurements of single-turnover rates and the temperature dependence of intrinsic kinetic isotope effects (KIEsint) of site-directed mutants. This study hypothesizes that the global network of residues involved in the chemical step is evolutionarily conserved and probes homologous residues of the potential global network in human DHFR through measurements of the temperature dependence of KIEsint and computer simulations based on the empirical valence bond method. We study mutants M53W and S145V. Both of these remote residues are homologous to network residues in E. coli DHFR. Non-additive isotope effects on activation energy are observed between M53 and S145, indicating their synergistic effect on the chemical step in human DHFR, which suggests that both of these residues are part of a network affecting the chemical step in enzyme catalysis. This finding supports the hypothesis that human and E. coli DHFR share similar networks, consistent with evolutionary preservation of such networks.

Alkyl hydrogen atom abstraction reactions of the CN radical with ethanol.

We present a study of the abstraction of alkyl hydrogen atoms from the β and α positions of ethanol by the CN radical in solution using the Empirical Valence Bond (EVB) method. We have built separate 2 × 2 EVB models for the Hβ and Hα reactions, where the atom transfer is parameterized using ab initio calculations. The intra- and intermolecular potentials of the reactant and product molecules were modelled with the General AMBER Force Field, with some modifications. We have carried out the dynamics in water and chloroform, which are solvents of contrasting polarity. We have computed the potential of mean force for both abstractions in each of the solvents. They are found to have a small and early barrier along the reaction coordinate with a large energy release. Analyzing the solvent structure around the reaction system, we have found two solvents to have little effect on either reaction. Simulating the dynamics from the transition state, we also study the fate of the energies in the HCN vibrational modes. The HCN molecule is born vibrationally hot in the CH stretch in both reactions and additionally in the HCN bends for the Hα abstraction reaction. In the early stage of the dynamics, we find that the CN stretch mode gains energy at the expense of the energy in CH stretch mode.

Cold Adaptation of Triosephosphate Isomerase.

The main problem for enzymes from psychrophilic species, which need to work near the freezing point of liquid water, is the exponential decay of reaction rates as the temperature is decreased. Cold-adapted enzymes have solved this problem by shifting the activation enthalpy-entropy balance for the catalyzed reaction compared to those of their mesophilic orthologs. To understand the structural basis of this universal feature, it is necessary to examine pairs of such orthologous enzymes, with known three-dimensional structures, at the microscopic level. Here, we use molecular dynamics free energy calculations in combination with the empirical valence bond method to evaluate the temperature dependence of the activation free energy for differently adapted triosephosphate isomerases. The results show that the enzyme from the psychrophilic bacterium Vibrio marinus indeed displays the characteristic shift in enthalpy-entropy balance, compared to that of the yeast ortholog. The origin of this effect is found to be located in a few surface-exposed protein loops that show differential mobilities in the two enzymes. Key mutations render these loops more mobile in the cold-adapted triosephosphate isomerase, which explains both the reduced activation enthalpy contribution from the protein surface and the lower thermostability.

A comparison of sodium and hydrogen halides at the air-water interface.

New molecular models, parameterized to ab initio calculations, were developed to describe HBr and HI at the air-water interface. These were used to compare how the air-water interface influenced dissociation of NaX and HX, with X being Cl, Br, or I, and also their propensity for the interface. The polarizable multistate empirical valence bond method, which explicitly describes proton sharing, was used to model HX. Results showed that the air-water interface suppressed HX dissociation from a contact ion pair to a solvent separated to a greater degree than NaX dissociation. Furthermore, HX had a greater propensity for the interface than NaX, which was a consequence of the hydronium ion having a greatest interfacial activity of all species studied. As a consequence of this, the average configuration of dissociated HX, while in both contact ion and solvent separated ion pairs near the air-water interface, is with the dissociated hydrogen oriented more towards the air than the X atom.

Thermodynamics of the Purine Nucleoside Phosphorylase Reaction Revealed by Computer Simulations

Enzymes are able to catalyze chemical reactions by reducing the activation free energy, yielding significant increases in the reaction rates. This can thermodynamically be accomplished by either reducing the activation enthalpy or increasing the activation entropy. The effect of remote mutations on the thermodynamic activation parameters of human purine nucleoside phosphorylase is examined using extensive molecular dynamics and free energy simulations. More than 2700 independent reaction free energy profiles for six different temperatures have been calculated to obtain high-precision computational Arrhenius plots. On the basis of these, the activation enthalpies and entropies were computed from linear regression of the plots with ΔG⧧ as a function of 1/T, and the obtained thermodynamic activation parameters are in very good agreement with those from experiments. The Arrhenius plots immediately show that the 6-oxopurines (INO and GUO) have identical slopes, whereas the 6-aminopurine (ADO) has a significantly different slope, indicating that the substrate specificity is related to the difference in thermodynamic activation parameters. Furthermore, the calculations show that the human PNP specificity for 6-oxopurines over 6-aminopurines originates from significant differences in electrostatic preorganization. The effect of the remote double mutation, K22E and H104R (E:R), has also been examined, as it alters human PNP toward the bovine PNP. These residues are situated on the protein surface, 28-35 Å from the active site, and the mutation alters the enthalpy-entropy balance with little effect on the catalytic rates. It is thus quite remarkable that the empirical valence bond method can reproduce the enthalpies and entropies induced by these long-range mutations.

Potentials of mean force in acidic proton transfer reactions in constrained geometries

Free energy barriers associated with the transfer of an excess proton in water and related to the potentials of mean force in proton transfer episodes have been computed in a wide range of thermodynamic states, from low-density amorphous ices to high-temperature liquids under the critical point for unconstrained and constrained systems. The latter were represented by set-ups placed inside hydrophobic graphene slabs at the nanometric scale allocating a few water layers, namely one or two in the narrowest case. Water–proton and carbon–proton forces were modelled with a Multi-State Empirical Valence Bond method. As a general trend, a competition between the effects of confinement and temperature is observed on the local hydrogen-bonded structures around the lone proton and, consequently, on the mean force exerted by its environment on the water molecule carrying the proton. Free energy barriers estimated from the computed potentials of mean force tend to rise with the combined effect of increasing temperatures and the packing effect due to a larger extent of hydrophobic confinement. The main reason observed for such enhancement of the free energy barriers was the breaking of the second coordination shell around the lone proton.

Molecular Dynamics Study on Proton Transport in Supported Nafion Ionomer Thin Films on Lennard-Jones Walls

The catalyst layer, one of the components of membrane electrode assembly in polymer electrolyte fuel cells, was modeled by molecular dynamics methods, and proton transport in ionomer thin films was simulated to elucidate the effect of surface wetting interaction on the ionomer thin film morphology and proton transport property. Ionomer thin film properties are considered to be affected by the surface property of supported carbon due to their nano-scale thickness (4-10 nm thick). One of the surface properties is surface wetting interaction basically attributed to the physical or chemical property of surfaces.In this study, modeled Nafion thin films with 4 nm and 10 nm thick were generated on the two ideal LJ smooth surfaces which have different strength of interactions to focus on the effect of the surface wetting interaction on the ionomer thin film morphology and proton transport. Protons exist as hydronium ions in ionomer thin films for chemical stability and are transported by the two transport mechanisms which are the vehicle mechanism (just move as hydronium ions) and the Grotthuss mechanism called proton hopping. It is common that proton transport model is restricted to only the vehicle mechanism when using classical molecular dynamics method, however, our proton transport model includes the Grotthuss mechanism by using aTS-EVB method (anharmonic two-state empirical valence bond method) in order to make proton transport model more reproducible.From the analysis in the ionomer thin films with 4 nm thick, surface wetting interaction was found to affect the orientation of sulfonic group near the surface. The orientation is well-ordered on the higher hydrophobic surface, and the sulfonic groups point to the gas phase direction, while it is randomly distributed on the lower hydrophobic surface, which make differences in the water cluster connectivity and size in ionomer thin films between the two surfaces. Proton transport property was evaluated by the self-diffusion coefficient. The proton self-diffusion coefficient in ionomer thin films was indicated to be improved on higher hydrophobic surface probably because higher hydrophobic surface makes ideal water cluster structures for proton transport in ionomer thin films.As compared to the ionomers with the different film thickness, the proton self-diffusion coefficient in ionomer thin films on higher hydrophobic surface was indicated to decrease in the increase of the film thickness, while that on lower hydrophobic surface was indicated to increase slightly. This result is expected to come from the confinement and interface effects which can be remarkable when the thickness decrease. The ionomer morphology on higher hydrophobic surface has advantages for proton transport probably thanks to the effects as observed in the ionomer with 4 nm thick, while that on lower hydrophobic surfaces has disadvantage. Since the effects can be weaker with the increase in the film thickness, the proton self-diffusion coefficient in ionomer thin films on higher hydrophobic surface was found to decrease with the increase in the film thickness, while that on lower hydrophobic surface was found to increase slightly.

Multistate empirical valence bond study of temperature and confinement effects on proton transfer in water inside hydrophobic nanochannels.

Microscopic characteristics of an aqueous excess proton in a wide range of thermodynamic states, from low density amorphous ices (down to 100 K) to high temperature liquids under the critical point (up to 600 K), placed inside hydrophobic graphene slabs at the nanometric scale (with interplate distances between 3.1 and 0.7 nm wide) have been analyzed by means of molecular dynamics simulations. Water-proton and carbon-proton forces were modeled with a multistate empirical valence bond method. Densities between 0.07 and 0.02 Å(-3) have been considered. As a general trend, we observed a competition between effects of confinement and temperature on structure and dynamical properties of the lone proton. Confinement has strong influence on the local structure of the proton, whereas the main effect of temperature on proton properties is observed on its dynamics, with significant variation of proton transfer rates, proton diffusion coefficients, and characteristic frequencies of vibrational motions. Proton transfer is an activated process with energy barriers between 1 and 10 kJ/mol for both proton transfer and diffusion, depending of the temperature range considered and also on the interplate distance. Arrhenius-like behavior of the transfer rates and of proton diffusion are clearly observed for states above 100 K. Spectral densities of proton species indicated that in all states Zundel-like and Eigen-like complexes survive at some extent. © 2016 Wiley Periodicals, Inc.

Importance of hydrophobic traps for proton diffusion in lyotropic liquid crystals.

The diffusion of protons in self-assembled systems is potentially important for the design of efficient proton exchange membranes. In this work, we study proton dynamics in a low-water content, lamellar phase of a sodium-carboxylate gemini surfactant/water system using computer simulations. The hopping of protons via the Grotthuss mechanism is explicitly allowed through the multi-state empirical valence bond method. We find that the hydronium ion is trapped on the hydrophobic side of the surfactant-water interface, and proton diffusion then proceeds by hopping between surface sites. The importance of hydrophobic traps is surprising because one would expect the hydronium ions to be trapped at the charged headgroups. The physics illustrated in this system should be relevant to the proton dynamics in other amphiphilic membrane systems, whenever there exist exposed hydrophobic surface regions.

Effects of Water Structure on Proton Transport in Nafion Thin Films with Molecular Dynamics Simulations

Proton transport properties in Nafion thin films have been investigated using molecular dynamics simulations to characterize the nanoscopic flow phenomena in the electrode of polymer electrolyte fuel cell. Protons transfer in the ionomer by a combination of the Vehicle mechanism and the Grotthuss mechanism. In this study, the Grotthuss mechanism is described using the anharmonic two-state empirical valence bond method to evaluate quantitative proton transport within the framework of classical molecular dynamics simulations. The electrode consists of platinum-supported carbon and ionomers, and it is well known that the wettability of supported carbon depends on materials and operating environment. Our results suggested that the wettability affects morphology of the ionomer and proton transport.