Multiscale electrostatic embedding simulations for modeling structure and dynamics of molecules in solution: A tutorial review

The main concepts and important technical details of electrostatic embedding quantum mechanics/molecular mechanics (QM/MM) simulations are explained and illustrated with the intent of assisting newcomers in performing and gauging the accuracy of such simulations, focused on smaller molecules in solution. Beginners are advised on how to increase the reliability and accuracy of the simulations through benchmarking. Central considerations on methodologies for QM/MM Molecular Dynamics (MD) simulations are presented, alongside technical fundamentals regarding the construction and manipulation of simulation systems using the python-based Atomic Simulation Environment (ASE). A worked example of QM/MM Born–Oppenheimer MD is included, and a flowchart summarizing the most salient decisions and tasks within the methodology is presented.     

An extensible interface for QM/MM molecular dynamics simulations with AMBER

We present an extensible interface between the AMBER molecular dynamics (MD) software package and electronic structure software packages for quantum mechanical (QM) and mixed QM and classical molecular mechanical (MM) MD simulations within both mechanical and electronic embedding schemes. With this interface, ab initio wave function theory and density functional theory methods, as available in the supported electronic structure software packages, become available for QM/MM MD simulations with AMBER. The interface has been written in a modular fashion that allows straight forward extensions to support additional QM software packages and can easily be ported to other MD software. Data exchange between the MD and QM software is implemented by means of files and system calls or the message passing interface standard. Based on extensive tests, default settings for the supported QM packages are provided such that energy is conserved for typical QM/MM MD simulations in the microcanonical ensemble. Results for the free energy of binding of calcium ions to aspartate in aqueous solution comparing semiempirical and density functional Hamiltonians are shown to demonstrate features of this interface. © 2013 Wiley Periodicals, Inc.     

Dynamic structures of phosphodiesterase-5 active site by combined molecular dynamics simulations and hybrid quantum mechanical/molecular mechanical calculations

Various quantum mechanical/molecular mechanical (QM/MM) geometry optimizations starting from an x-ray crystal structure and from the snapshot structures of constrained molecular dynamics (MD) simulations have been performed to characterize two dynamically stable active site structures of phosphodiesterase-5 (PDE5) in solution. The only difference between the two PDE5 structures exists in the catalytic, second bridging ligand (BL2) which is HO- or H2O. It has been shown that, whereas BL2 (i.e. HO-) in the PDE5(BL2 = HO-) structure can really bridge the two positively charged metal ions (Zn2+ and Mg2+), BL2 (i.e. H2O) in the PDE5(BL2 = H2O) structure can only coordinate Mg2+. It has been demonstrated that the results of the QM/MM geometry optimizations are remarkably affected by the solvent water molecules, the dynamics of the protein environment, and the electronic embedding charges of the MM region in the QM part of the QMM/MM calculation. The PDE5(BL2 = H2O) geometries optimized by using the QM/MM method in different ways show strong couplings between these important factors. It is interesting to note that the PDE5(BL2 = HO-) and PDE5(BL2 = H2O) geometries determined by the QM/MM calculations neglecting these three factors are all consistent with the corresponding geometries determined by the QM/MM calculations that account for all of these three factors. These results suggest the overall effects of these three important factors on the optimized geometries can roughly cancel out. However, the QM/MM calculations that only account for some of these factors could lead to considerably different geometries. These results might be useful also in guiding future QM/MM geometry optimizations on other enzymes.     

Combined QM(DFT)/MM molecular dynamics simulations of the deamination of cytosine by yeast cytosine deaminase (yCD)

Extensive combined quantum mechanical (B3LYP/6‐31G*) and molecular mechanical (QM/MM) molecular dynamics simulations have been performed to elucidate the hydrolytic deamination mechanism of cytosine to uracil catalyzed by the yeast cytosine deaminase (yCD). Though cytosine has no direct binding to the zinc center, it reacts with the water molecule coordinated to zinc, and the adjacent conserved Glu64 serves as a general acid/base to shuttle protons from water to cytosine. The overall reaction consists of several proton‐transfer processes and nucleophilic attacks. A tetrahedral intermediate adduct of cytosine and water binding to zinc is identified and similar to the crystal structure of yCD with the inhibitor 2‐pyrimidinone. The rate‐determining step with the barrier of 18.0 kcal/mol in the whole catalytic cycle occurs in the process of uracil departure where the proton transfer from water to Glu64 and nucleophilic attack of the resulting hydroxide anion to C2 of the uracil ring occurs synchronously. © 2016 Wiley Periodicals, Inc.     

Modeling of phytochrome absorption spectra

Phytochromes constitute one of the six well‐characterized families of photosensory proteins in Nature. From the viewpoint of computational modeling, however, phytochromes have been the subject of much fewer studies than most other families of photosensory proteins, which is likely a consequence of relevant high‐resolution structural data becoming available only in recent years. In this work, hybrid quantum mechanics/molecular mechanics (QM/MM) methods are used to calculate UV‐vis absorption spectra of Deinococcus radiodurans bacteriophytochrome. We investigate how the choice of QM/MM methodology affects the resulting spectra and demonstrate that QM/MM methods can reproduce the experimental absorption maxima of both the Q and Soret bands with an accuracy of about 0.15 eV. Furthermore, we assess how the protein environment influences the intrinsic absorption of the bilin chromophore, with particular focus on the Q band underlying the primary photochemistry of phytochromes. © 2013 Wiley Periodicals, Inc.     

Revisiting the structure and dynamics of hydrated Cd2+ in aqueous solutions: Insights from the RI-SCS-MP2/MM molecular dynamics simulation

The spin component scale MP2/molecular mechanics molecular dynamics simulation investigated the hydration shell formation and hydrated Cd²⁺ dynamics in the water environment. At the first hydration shell, six water molecules with 2.27 Å for the average distance between water and Cd²⁺. Dynamical properties were analyzed by computing the water molecule's mean residence time (MRT) in its first and second hydration shells. The MRT of each shell was determined to be 31.8 and 1.92 ps, suggesting the strong influence of Cd²⁺ in the first hydration shell. The second shell was labile, with an average number of water molecules being 18. Despite the strong interaction between Cd²⁺ and water molecules in the first shell, the influence of ions in the second hydration shell remained weak.     

Flooding in GROMACS: accelerated barrier crossings in molecular dynamics

The major bottleneck of today's atomistic molecular dynamics (MD) simulations is that because of the enormous computational effort involved, only processes at nanoseconds to microseconds time scales or faster can be studied directly. Unfortunately, apart from a few exceptions, relevant processes, such as chemical reactions or many large scale conformational transitions in proteins, occur at slower time scales and therefore are currently far out of reach for conventional MD. The flooding technique addresses this problem by inclusion of a flooding potential into the force field. This flooding potential locally destabilizes the educt state and thereby significantly accelerates the escape from the initial energy well without affecting the reaction pathway. Here, we summarize the theory and method for the computational chemistry community and detail the implementation within the official version 3.3 of the freely available MD program package GROMACS. Two examples shall demonstrate the application of flooding to accelerate conformational transitions and chemical reactions. The second example was carried out within a QM/MM framework.     

Molecular dynamics simulations of the detoxification of paraoxon catalyzed by phosphotriesterase

Combined QM(PM3)/MM molecular dynamics simulations together with QM(DFT)/MM optimizations for key configurations have been performed to elucidate the enzymatic catalysis mechanism on the detoxification of paraoxon by phosphotriesterase (PTE). In the simulations, the PM3 parameters for the phosphorous atom were reoptimized. The equilibrated configuration of the enzyme/substrate complex showed that paraoxon can strongly bind to the more solvent‐exposed metal ion Zn_β, but the free energy profile along the binding path demonstrated that the binding is thermodynamically unfavorable. This explains why the crystal structures of PTE with substrate analogues often exhibit long distances between the phosphoral oxygen and Zn_β. The subsequent S_N2 reaction plays the key role in the whole process, but controversies exist over the identity of the nucleophilic species, which could be either a hydroxide ion terminally coordinated to Zn_α or the μ‐hydroxo bridge between the α‐ and β‐metals. Our simulations supported the latter and showed that the rate‐limiting step is the distortion of the bound paraoxon to approach the bridging hydroxide. After this preparation step, the bridging hydroxide ion attacks the phosphorous center and replaces the diethyl phosphate with a low barrier. Thus, a plausible way to engineer PTE with enhanced catalytic activity is to stabilize the deformed paraoxon. Conformational analyses indicate that Trp131 is the closest residue to the phosphoryl oxygen, and mutations to Arg or Gln or even Lys, which can shorten the hydrogen bond distance with the phosphoryl oxygen, could potentially lead to a mutant with enhanced activity for the detoxification of organophosphates. © 2009 Wiley Periodicals, Inc. J Comput Chem, 2009     

Reaching multi‐nanosecond timescales in combined QM/MM molecular dynamics simulations through parallel horsetail sampling

We report an enhanced sampling technique that allows to reach the multi‐nanosecond timescale in quantum mechanics/molecular mechanics molecular dynamics simulations. The proposed technique, called horsetail sampling, is a specific type of multiple molecular dynamics approach exhibiting high parallel efficiency. It couples a main simulation with a large number of shorter trajectories launched on independent processors at periodic time intervals. The technique is applied to study hydrogen peroxide at the water liquid–vapor interface, a system of considerable atmospheric relevance. A total simulation time of a little more than 6 ns has been attained for a total CPU time of 5.1 years representing only about 20 days of wall‐clock time. The discussion of the results highlights the strong influence of the solvation effects at the interface on the structure and the electronic properties of the solute. © 2017 Wiley Periodicals, Inc.     

molUP: A VMD plugin to handle QM and ONIOM calculations using the gaussian software

The notable advances obtained by computational (bio)chemistry provided its widespread use in many areas of science, in particular, in the study of reaction mechanisms. These studies involve a huge number of complex calculations, which are often carried out using the Gaussian suite of programs. The preparation of input files and the analysis of the output files are not easy tasks and often involve laborious and complex steps. Taking this into account, we developed molUP: a VMD plugin that offers a complete set of tools that enhance the preparation of QM and ONIOM (QM/MM, QM/QM, and QM/QM/MM) calculations. The starting structures for these calculations can be imported from different chemical formats. A set of tools is available to help the user to examine or modify any geometry parameter. This includes the definition of layers in ONIOM calculations, choosing fixed atoms during geometry optimizations, the recalculation or adjustment of the atomic charges, performing SCANs or IRC calculations, etc. molUP also extracts the geometries from the output files as well as the energies of each of them. All of these tasks are performed in an interactive GUI that is extremely helpful for the user. MolUP was developed to be easy to handle by inexperienced users, but simultaneously to be a fast and flexible graphical interface to allow the advanced users to take full advantage of this plugin. The program is available, free of charges, for macOS, Linux, and Windows at the PortoBioComp page https://www.fc.up.pt/PortoBioComp/database/doku.php?id=molup . © 2018 Wiley Periodicals, Inc.     

Increasing the time step with mass scaling in Born‐Oppenheimer ab initio QM/MM molecular dynamics simulations

Born‐Oppenheimer ab initio QM/MM molecular dynamics simulation with umbrella sampling is a state‐of‐the‐art approach to calculate free energy profiles of chemical reactions in complex systems. To further improve its computational efficiency, a mass‐scaling method with the increased time step in MD simulations has been explored and tested. It is found that by increasing the hydrogen mass to 10 amu, a time step of 3 fs can be employed in ab initio QM/MM MD simulations. In all our three test cases, including two solution reactions and one enzyme reaction, the resulted reaction free energy profiles with 3 fs time step and mass scaling are found to be in excellent agreement with the corresponding simulation results using 1 fs time step and the normal mass. These results indicate that for Born‐Oppenheimer ab initio QM/MM molecular dynamics simulations with umbrella sampling, the mass‐scaling method can significantly reduce its computational cost while has little effect on the calculated free energy profiles. © 2009 Wiley Periodicals, Inc. J Comput Chem, 2009     

Interfacing the GROMOS (bio)molecular simulation software to quantum‐chemical program packages

The newly implemented quantum‐chemical/molecular‐mechanical (QM/MM) functionality of the Groningen molecular simulation (GROMOS) software for (bio)molecular simulation is described. The implementation scheme is based on direct coupling of the GROMOS C++ software to executables of the quantum‐chemical program packages MNDO and TURBOMOLE, allowing for an independent further development of these packages. The new functions are validated for different test systems using program and model testing techniques. The effect of truncating the QM/MM electrostatic interactions at various QM/MM cutoff radii is discussed and the application of semiempirical versus density‐functional Hamiltonians for a solute molecule in aqueous solution is compared. © 2012 Wiley Periodicals, Inc.     

Influence of hydrostatic pressure on water absorption of polyoxymethylene: experiment and molecular dynamics simulation

This work experimentally investigated the influence of hydrostatic pressure on the water absorption of polyoxymethylene, followed by analyzing its micro‐mechanism via molecular dynamics simulation. Tests results show that the polyoxymethylene water absorption decreases with the increase in hydrostatic pressure when the latter is within 0–3.0 MPa; it subsequently increases with the increase in hydrostatic pressure in the range of 3.0–5.0 MPa. Simulation of water molecules diffusion on polyoxymethylene surface shows that water molecules diffuse into polyoxymethylene surface during an equilibration run, and water molecule displacement value of maximum relative concentration gives almost the same characteristics: firstly decreases and then increases with the increase of hydrostatic pressure. Simulation of water molecule diffusion inside polyoxymethylene suggests the following: (i) water molecules vibrate from the interior to the edge of polyoxymethylene cell during the equilibration run, and (ii) water diffusion coefficient in polyoxymethylene gives trend of firstly decreasing and then increasing with the increase of hydrostatic pressure. Copyright © 2016 John Wiley & Sons, Ltd.     

Analytic first and second derivatives of the energy in the fragment molecular orbital method combined with molecular mechanics

Analytic first and second derivatives of the energy are developed for the fragment molecular orbital method interfaced with molecular mechanics in the electrostatic embedding scheme at the level of Hartree-Fock and density functional theory. The importance of the orbital response terms is demonstrated. The role of electrostatic embedding upon molecular vibrations is analyzed, comparing force field and quantum mechanical treatments for an ionic liquid and a solvated protein. The method is applied for 100 protein conformations sampled in molecular dynamics (MD) to take into account the complexity of a flexible protein structure in solution, and a good agreement with experimental data is obtained: Frequencies from an experimental infrared (IR) spectrum are reproduced within 17 cm⁻¹ .     

Structure and dynamics of hydrated NH(4) (+): an ab initio QM/MM molecular dynamics simulation

A combined ab initio quantum mechanical/molecular mechanical (QM/MM) molecular dynamics simulation has been performed to investigate solvation structure and dynamics of NH(4) (+) in water. The most interesting region, the sphere includes an ammonium ion and its first hydration shell, was treated at the Hartree-Fock level using DZV basis set, while the rest of the system was described by classical pair potentials. On the basis of detailed QM/MM simulation results, the solvation structure of NH(4) (+) is rather flexible, in which many water molecules are cooperatively involved in the solvation shell of the ion. Of particular interest, the QM/MM results show fast translation and rotation of NH(4) (+) in water. This phenomenon has resulted from multiple coordination, which drives the NH(4) (+) to translate and rotate quite freely within its surrounding water molecules. In addition, a "structure-breaking" behavior of the NH(4) (+) is well reflected by the detailed analysis on the water exchange process and the mean residence times of water molecules surrounding the ion.     

Response properties of embedded molecules through the polarizable embedding model

The polarizable embedding (PE) model is a fragment-based quantum-classical approach aimed at accurate inclusion of environment effects in quantum-mechanical response property calculations. The aim of this tutorial review is to give insight into the practical use of the PE model. Starting from a set of molecular structures and until you arrive at the final property, there are many crucial details to consider in order to obtain trustworthy results in an efficient manner. To lower the threshold for new users wanting to explore the use of the PE model, we describe and discuss important aspects related to its practical use. This includes directions on how to generate input files and how to run a calculation.     

Understanding the complete bioluminescence cycle from a multiscale computational perspective: A review

Bioluminescence (BL) is an amazing natural phenomenon whose visible light is produced by living organisms. BL phenomenon is quite pervasive and has been observed in 17 phyla of 4 kingdoms. This fascinating natural phenomenon has unceasingly attracted people’s curiosity from ancient era to today. For a very long time, we can only receive some sporadic and static information from experimental observations, the mechanism of most BL remains is unclear. How the chemical reaction of BL process is initiated? Where the energy for light emission comes from? How does the light emitter produce? What is the light emitter for a wild bioluminescent organism? How to regain luciferin for next bioluminescence when it is used up? The luciferin is utilized forthwith or stored and release for subsequent light emission? What factors affect the color and strength of a bioluminescence? How to artificially tune the bioluminescence for special application? Computational BL plays unreplaceable role in answering these mechanistic questions. In contrast with experimental BL, computational BL came very late. In the past two decades, computational BL has touched nearly all the bioluminescent systems with chemical bases via the method of multiscale simulation. In this review, the author firstly introduced the history, types and general chemical process of BL. Then, the computational scheme on BL was briefly epitomized. Using firefly BL as a paradigmatic case, the author summarized theoretical investigation on the six stages of general chemical process in a BL cycle: luciferin oxidation, peroxide thermolysis, light emission, luciferin regeneration, luciferin storage and luciferin release. At each stage, the available theoretical studies of other bioluminescent organisms are briefly introduced and compared with the firefly system. Basing on the mechanistic understanding, the author reviewed the up-to-date theoretical design on bioluminescent systems. Again, the firefly was mainly focused on, and the other possible systems were just briefly introduced. This review summarized the theoretical studies to date on BL and addressed the status, critical challenges and future prospects of computational BL.     

An averaged solvent electrostatic potential from molecular dynamics study of the anomeric equilibrium of D-xylose in aqueous solution

We applied the free-energy perturbation method together with the averaged solvent electrostatic potential from molecular dynamics (ASEP/MD) method to study the anomeric equilibrium of d-xylose in aqueous solution. The level of calculation, 6-311G++(2d,2p) basis set and density functional theory, permits one to explain the main characteristics of the anomeric equilibrium of d-xylopyranose: in vacuo, the anomeric effect predominates and the form is the stabler. In water, solvation leads to the form being the stabler. A comparison between the performances of the ASEP/MD and polarizable continuum models is also presented.Contribution to the Jacopo Tomasi Honorary Issue     

Protonation Equilibrium in the Active Site of the Photoactive Yellow Protein

The role and existence of low-barrier hydrogen bonds (LBHBs) in enzymatic and protein activity has been largely debated. An interesting case is that of the photoactive yellow protein (PYP). In this protein, two short HBs adjacent to the chromophore, p-coumaric acid (pCA), have been identified by X-ray and neutron diffraction experiments. However, there is a lack of agreement on the chemical nature of these H-bond interactions. Additionally, no consensus has been reached on the presence of LBHBs in the active site of the protein, despite various experimental and theoretical studies having been carried out to investigate this issue. In this work, we perform a computational study that combines classical and density functional theory (DFT)-based quantum mechanical/molecular mechanical (QM/MM) simulations to shed light onto this controversy. Furthermore, we aim to deepen our understanding of the chemical nature and dynamics of the protons involved in the two short hydrogen bonds that, in the dark state of PYP, connect pCA with the two binding pocket residues (E46 and Y42). Our results support the existence of a strong LBHB between pCA and E46, with the H fully delocalized and shared between both the carboxylic oxygen of E46 and the phenolic oxygen of pCA. Additionally, our findings suggest that the pCA interaction with Y42 can be suitably described as a typical short ionic H-bond of moderate strength that is fully localized on the phenolic oxygen of Y42.