Surface-Driven ElectronLocalization and Defect Heterogeneityin Ceria

Introduction: Biological macromolecules, such as proteins or nucleic acids, are (still) molecules and thus they follow the same chemical rules that any simple molecule follows, even if their size generally renders accurate studies unhelpful. However, in the context of drug discovery, a detailed analysis of ligand association is required for understanding or predicting their interactions and hybrid quantum mechanics/molecular mechanics (QM/MM) computations are relevant tools to help elucidate this process.Areas covered: In this review, the authors explore the use of QM/MM for drug discovery. After a brief description of the molecular mechanics (MM) technique, the authors describe the subtractive and additive techniques for QM/MM computations. The authors then present several application cases in topics involved in drug discovery.Expert opinion: QM/MM have been widely employed during the last decades to study chemical processes such as enzyme–inhibitor interactions. However, despite the enthusiasm around this area, plain MM simulations may be more meaningful than QM/MM. To obtain reliable results, the authors suggest fixing several keystone parameters according to the underlying chemistry of each studied system.

Butyrylcholinesterase (BChE) is considered as an efficient stoichiometric antidote against organophosphorus (OP) poisons. Recently we utilized combination of calculations and ultrahigh-throughput screening (uHTS) to select BChE variants capable of catalytic destruction of OP pesticide paraoxon. The purpose of this study was to elucidate the molecular mechanism underlying enzymatic hydrolysis of paraoxon by BChE variants using hybrid quantum mechanical/molecular mechanical (QM/MM) calculations. Detailed analysis of accomplished QM/MM runs revealed that histidine residues introduced into the acyl-binding loop are always located in close proximity with aspartate residue at position 70. Histidine residue acts as general base thus leading to attacking water molecule activation and subsequent SN2 inline hydrolysis resulting in BChE reactivation. This combination resembles canonical catalytic triad found in active centers of various proteases. Carboxyl group activates histidine residue by altering its pKa, which in turn promotes the activation of water molecule in terms of its nucleophilicity. Observed re-protonation of catalytic serine residue at position 198 from histidine residue at position 438 recovers initial configuration of the enzyme’s active center, facilitating next catalytic cycle. We therefore suggest that utilization of uHTS platform in combination with deciphering of molecular mechanisms by QM/MM calculations may significantly improve our knowledge of enzyme function, propose new strategies for enzyme design and open new horizons in generation of catalytic bioscavengers against OP poisons.

Protein splicing is a post-translational process in which a biologically inactive protein is activated by the release of a segment denoted as an intein. The process involves four steps. In the third, the scission of the intein takes place after the cyclization of the last amino acid of the segment, an asparagine. Little is known about the chemical reaction necessary for this cyclization. Experiments demonstrate that two histidines (the penultimate amino acid of the intein, and a histidine located 10 amino acids upstream) are relevant in the cyclization of the asparagine. We have investigated the mechanism and determinants of reaction in the GyrA intein focusing on the requirements for asparagine activation for its cyclization. First, the influence that the protonation states of these two histidines have on the orientation of the asparagine side chain is investigated by means of molecular dynamics simulation. Molecular dynamics simulations using the CHARMM27 force field were carried out on the three possible protonation states for each of these two histidines. The results indicate that the only protonation state in which the conformation of the system is suitable for cyclization is when the penultimate histidine is fully protonated (positively charged), and the upstream histidine is in the His(ε) neutral tautomeric form. The free energy profile for the reaction in which the asparagine is activated by a proton transfer to the upstream histidine is presented, computed by hybrid quantum mechanics/molecular mechanics (QM/MM) umbrella sampling molecular dynamics at the SCCDFTB/CHARMM27 level of theory. The calculated free energy barrier for the reaction is 19.0 kcal mol(-1). B3LYP/6-31+G(d) QM/MM single-point calculations give a qualitatively a similar energy profile, although with somewhat higher energy barriers, in good agreement with the value derived from experiment of 25 kcal mol(-1) at 60 °C. QM/MM molecular dynamics simulations of the reactant, activated reactant and intermediate states highlight the importance of the Arg181-Val182-Asp183 segment in catalysing the reaction. Overall, the results indicate that nucleophilic activation of the asparagine for its cyclization by the upstream histidine acting as the base is a plausible mechanism for the C-terminal cleavage in protein splicing.

Thiols are widely present in biological systems, most notably as the side chain of cysteine amino acids in proteins. Thiols can be deprotonated to form a thiolate which affords a diverse range of enzymatic activity and modes for chemical modification of proteins. Parameters for modeling thiolates using molecular mechanical force fields have not yet been validated, in part due to the lack of structural data on thiolate solvation. Here, the CHARMM36 and Amber models for thiolates in aqueous solutions are assessed using free energy perturbation and hybrid quantum mechanics/molecular mechanics (QM/MM) molecular dynamics (MD) simulations. The hydration structure of methylthiolate was calculated from 1 ns of QM/MM MD (PBE0-D3/def2-TZVP//TIP3P), which shows that the water-S- distances are approximately 2 Å with a coordination number near 6. The CHARMM thiolate parameters predict a thiolate S radius close to the QM/MM value and predict a hydration Gibbs energy of -329.2 kJ/mol, close to the experimental value of -318 kJ/mol. The cysteine thiolate model in the Amber force field underestimates the thiolate radius by 0.2 Å and overestimates the thiolate hydration energy by 119 kJ/mol because it uses the same Lennard-Jones parameters for thiolates as for thiols. A recent Drude polarizable model for methylthiolate with optimized thiolate parameters also performs well. SAPT2+ [Symmetry Adapted Perturbation Theory (SAPT)] analysis indicates that exchange repulsion is larger for the methylthiolate, consistent with it having a more diffuse electron density distribution in comparison with the parent thiol. These data demonstrate that it is important to define distinct non-bonded parameters for the protonated/deprotonated states of amino acid side chains in molecular mechanical force fields.

Characterization of a Catalytic Ligand Bridging Metal Ions in Phosphodiesterases 4 and 5 by Molecular Dynamics Simulations and Hybrid Quantum Mechanical/Molecular Mechanical Calculations

Hybrid quantum mechanics/molecular mechanics (QM/MM) methods were employed to study electronic effects of 5,5'-substituents on the enantioselectivity of (salen)Mn catalysts. These methods were shown to effectively capture electronic effects of substituents if the substituents were treated with a simple electrostatic coupling using point charges. Correlations between calculated properties of (salen)Mn=O indicative of its reactivity and the modified Hammett constant σ + para were applied to study electronic framework effects on enantioselectivity for (salen)Mn heterogenized in a metal—organic framework (MOF). The results suggest that these effects are too small to account for all of the decrease in enantioselectivity observed experimentally and that steric constraints of the framework likely play a role. The reactivity correlations also enabled fast ligand screening for the prediction of highly enantioselective (salen)Mn catalysts that could be used in MOFs. A modified version of the DREIDING force field was also developed to allow QM/MM calculations of (salen)Mn species.

MurD (UDP-N-acetylmuramoyl-L-alanine:D-glutamate ligase), a three-domain bacterial protein, catalyses a highly specific incorporation of D-glutamate to the cytoplasmic intermediate UDP-N-acetyl-muramoyl-L-alanine (UMA) utilizing ATP hydrolysis to ADP and P(i). This reaction is part of a biosynthetic path yielding bacterial peptidoglycan. On the basis of structural studies of MurD complexes, a stepwise catalytic mechanism was proposed that commences with a formation of the acyl-phosphate intermediate, followed by a nucleophilic attack of D-glutamate that, through the formation of a tetrahedral reaction intermediate and subsequent phosphate dissociation, affords the final product, UDP-N-acetyl-muramoyl-L-alanine-D-glutamate (UMAG). A hybrid quantum mechanical/molecular mechanical (QM/MM) molecular modeling approach was utilized, combining the B3LYP QM level of theory with empirical force field simulations to evaluate three possible reaction pathways leading to tetrahedral intermediate formation. Geometries of the starting structures based on crystallographic experimental data and tetrahedral intermediates were carefully examined together with a role of crucial amino acids and water molecules. The replica path method was used to generate the reaction pathways between the starting structures and the corresponding tetrahedral reaction intermediates, offering direct comparisons with a sequential kinetic mechanism and the available structural data for this enzyme. The acquired knowledge represents new and valuable information to assist in the ongoing efforts leading toward novel inhibitors of MurD as potential antibacterial drugs.

We present classical molecular dynamics (MD), Born-Oppenheimer molecular dynamics (BOMD), and hybrid quantum mechanics/molecular mechanics (QM/MM) data. MD was performed using the GPU accelerated pmemd module of the AMBER14MD package. BOMD was performed using CP2K version 2.6. The reaction rates in BOMD were accelerated using the Metadynamics method. QM/MM was performed using ONIOM in the Gaussian09 suite of programs. Relevant input files for BOMD and QM/MM are available.

Pathogen resistance to carbapenem antibiotics compromises effective treatments of superbug infections. One major source of carbapenem resistance is the bacterial production of carbapenemases which effectively hydrolyze carbapenem drugs. In this computational study, the deacylation reaction of imipenem (IPM) by GES-5 carbapenemases (GES) is modeled to unravel the mechanistic factors that facilitate carbapenem resistance. Hybrid quantum mechanical/molecular mechanical (QM/MM) calculations are applied to sample the GES/IPM deacylation barriers on the minimum energy pathways (MEPs). In light of the recent emergence of graph-based deep-learning techniques, we construct graph representations of the GES/IPM active site. An edge-conditioned graph convolutional neural network (ECGCNN) is trained on the acyl-enzyme conformational graphs to learn the underlying correlations between the GES/IPM conformations and the deacylation barriers. A perturbative approach is proposed to interpret the latent representations from the graph-learning (GL) model and extract essential mechanistic understanding with atomistic detail. In general, our study combining QM/MM MEPs calculations and GL models explains mechanistic landscapes underlying the IPM resistance driven by GES carbapenemases. We also demonstrate that GL methods could effectively assist the post-analysis of QM/MM calculations whose data span high dimensionality and large sample-size.

In this work, we apply quantum mechanics/molecular mechanics (QM/MM) approach to predict core-electron binding energies and chemical shifts of polymers, obtainable via X-ray photoelectron spectroscopy(XPS), using polymethyl methacrylate as a demonstration example. The results indicate that standard parametrizations of the quantum part (basis sets, level of correlation) and the molecular mechanics parts (decomposed charges, polarizabilities, and capping technique) are sufficient for the QM/MM model to be predictive for XPS of polymers. It is found that the polymer environment produces contributions to the XPS binding energies that are close to monotonous with the number of monomer units, totally amounting to approximately an eV decrease in binding energies. In most of the cases, the order of the shifts is maintained, and even the relative size of the differential shifts is largely preserved. The coupling of the internal core-hole relaxation to the polymer environment is found to be weak in each case, amounting only to one or two tenths of an eV. The main polymeric effect is actually well estimated already at the frozen orbital level of theory, which in turn implies a substantial computational simplification. These conclusions are best represented by the cases where the ionized monomer and its immediate surrounding are treated quantum mechanically. If the QM region includes only a single monomer, a couple of anomalies are spotted, which are referred to the QM/MM interface itself and to the neglect of a possible charge transfer.

ConspectusThe atmosphere is a key part of the earth system comprising myriad chemical species in all basic forms of matter. Ubiquitous nano- and microscopic aerosol particles and cloud droplets suspended in the air play crucial roles in earth's climate and the formation of air pollution. Surfaces are a prominent part of aerosols and droplets, due to the high surface area to bulk volume ratios, but very little is known about their specific properties. Many atmospheric compounds are surface-active, leading to enhanced surface concentrations in aqueous solutions. Their distribution between the surface and bulk may determine heterogeneous chemistry and many other properties of aerosol and cloud droplets, but has not been directly observed.We used X-ray photoelectron spectroscopy (XPS) to obtain direct molecular-level information on the surface composition and structure of aqueous solutions of surface-active organics as model systems for atmospheric aerosol and cloud droplets. XPS is a vacuum-based technique enabled for volatile aqueous organic samples by the application of a high-speed liquid microjet. In combination with brilliant synchrotron X-rays, the chemical specificity of XPS allows distinction between elements in different chemical states and positions within molecular structures. We used core-level C 1s and N 1s signals to identify the alkyl and hydrophilic groups of atmospheric carboxylic acids, alkyl-amines, and their conjugate acids and bases. From this, we infer changes in the orientation of surface-adsorbed species and quantify their relative abundances in the surface. XPS-derived surface enrichments of the organics follow trends expected from their surface activities and we observed a preferential orientation at the surface with the hydrophobic alkyl chains pointing increasingly outward from the solution at higher concentrations. This provides a first direct experimental observation of well-established concepts of surface adsorption and confirms the soundness of the method.We mapped relative abundances of conjugate acid-base pairs in the aqueous solution surfaces from the respective intensities of distinctive XPS signals. For each pair, the protonation equilibrium was significantly shifted toward the neutral form in the surface, compared to the bulk solution, across the full pH range. This represents an apparent shift of the pKa in the surface, which may be toward either higher or lower pH, depending on whether the acid or base form of the pair is the neutral species. The surface shifts are broadly consistent with the relative differences in surface enrichment of the individual acid and base conjugates in binary aqueous solutions, with additional contributions from nonideal interactions in the surface. In aqueous mixtures of surface-active carboxylate anions with ammonium salts at near-neutral pH, we found that the conjugate carboxylic acids were further strongly enhanced. This occurs as the coadsorption of weakly basic carboxylate anions and weakly acidic ammonium cations forms ion-pair surface layers with strongly enhanced local abundances, increasing the probability of net proton transfer according to Le Chatelier's principle. The effect is stronger when the evaporation of ammonia from the surface further contributes to irreversibly perturb the protonation equilibrium, leaving a surplus of carboxylic acid. These surface-specific effects may profoundly influence atmospheric chemistry mediated by aqueous aerosols and cloud droplets but are currently not taken into account in atmospheric models.

Enoate reductases catalyze the reduction of activated C═C bonds with high enantioselectivity. The oxidative half-reaction, which involves the addition of a hydride and a proton to opposite faces of the C═C bond, has been studied for the first time by hybrid quantum mechanics/molecular mechanics (QM/MM). The reduction of 2-cyclohexen-1-one by YqjM from Bacillus subtilis was selected as the model system. Two-dimensional QM/MM (B3LYP-D/OPLS2005) reaction pathways suggest that the hydride and proton are added as distinct steps, with the former step preceding the latter. Furthermore, we present interesting insights into the reactivity of this enzyme, including the weak binding of the substrate in the active site, the role of the two active site histidine residues for polarization of the substrate C═O bond, structural details of the transition states to hydride and proton transfer, and the role of Tyr196 as proton donor. The information presented here will be useful for the future design of enantioselective YqjM mutants for other substrates.

Quantum mechanical/molecular mechanical (QM/MM) geometry optimizations of the X-ray crystal structures of PDE10-AMP (PDB code 2OUN ) and PDE10-GMP (PDB code 2OUQ ) complexes have been performed to characterize the state of the AMP and GMP products, respectively. Results show that only one phosphate oxygen atom (O1) is protonated for both AMP and GMP product complexes. In addition, QM/MM calculations have resolved the orientation of the amide group of Gln726 in PDE10-GMP which was in conflict with the assignment of the guanine group of GMP in the X-ray crystal structure. Calculations reveal that the amide oxygen and nitrogen atom of Gln726 are rotated 180 degrees, resulting in two strong hydrogen bonds formed between the amide group of Gln726 and the guanine group of GMP. Binding free energy calculations for both QM/MM-optimized structures confirm the new conformational assignment of Gln726 in PDE10-GMP. The calculated binding free energy of the rotated structure is approximately 22 kcal/mol lower than the X-ray crystal assignment. The lower energy is mainly derived from the formation of two hydrogen bonds between the amide group of Gln726 and the guanine group of GMP. This implies that the orientation of the amide oxygen and nitrogen atoms in PDE10-AMP is different from PDE10-GMP. Finally, our results help to understand why PDE10 can hydrolyze both cAMP and cGMP.

Photoinduced processes play a crucial role in a multitude of important molecular phenomena. Accurately modeling these processes in an environment other than a vacuum requires a detailed description of the electronic states involved as well as how energy flows are coupled to the surroundings. Nonadiabatic effects must also be included in order to describe the exchange of energy between electronic and nuclear degrees of freedom correctly. In this work, we revisit the ring-opening reaction 1,3-cylohexadiene (CHD) in a solvent environment. Using our newly developed Interface for Non-Adiabatic Quantum mechanics/molecular mechanics in Solvent (INAQS) we trace the evolution of the reaction via hybrid quantum mechanics/molecular mechanics (QM/MM) surface hopping with a focus on the solvent's participation in the nonadiabatic relaxation process and the long-time approach to equilibrium. We explicitly include the MM solvent contribution to the nonadiabatic coupling vector─enabling an accurate approach to equilibrium at long times─and find that in highly multidimensional systems gradients can have little or nothing to do with the nonadiabatic couplings.

The interaction of hydrogen (H) with dislocations in tungsten (W) must be understood in order to model the mechanical response of future plasma-facing materials for fusion applications. Here, hybrid quantum mechanics/molecular mechanics (QM/MM) simulations are employed to study the ⟨111⟩ screw dislocation glide in W in the presence of H, using the virtual work principle to obtain energy barriers for dislocation glide, H segregation, and pipe diffusion. We provide a convincing validation of the QM/MM approach against full DFT energy-based methods. This is possible because the compact core and relatively weak elastic fields of ⟨111⟩ screw dislocations allow them to be contained in periodic DFT supercells. We also show that H segregation stabilizes the split-core structure while leaving the Peierls barrier almost unchanged. Furthermore, we find an energy barrier of less than 0.05 eV for pipe diffusion of H along dislocation cores. Our quantum-accurate calculations provide important reference data for the construction of larger-scale material models.