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.     

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.     

Importance of accurate charges in molecular docking: quantum mechanical/molecular mechanical (QM/MM) approach

The extent to which accuracy of electric charges plays a role in protein-ligand docking is investigated through development of a docking algorithm, which incorporates quantum mechanical/molecular mechanical (QM/MM) calculations. In this algorithm, fixed charges of ligands obtained from force field parameterization are replaced by QM/MM calculations in the protein environment, treating only the ligands as the quantum region. The algorithm is tested on a set of 40 cocrystallized structures taken from the Protein Data Bank (PDB) and provides strong evidence that use of nonfixed charges is important. An algorithm, dubbed "Survival of the Fittest" (SOF) algorithm, is implemented to incorporate QM/MM charge calculations without any prior knowledge of native structures of the complexes. Using an iterative protocol, this algorithm is able in many cases to converge to a nativelike structure in systems where redocking of the ligand using a standard fixed charge force field exhibits nontrivial errors. The results demonstrate that polarization effects can play a significant role in determining the structures of protein-ligand complexes, and provide a promising start towards the development of more accurate docking methods for lead optimization applications.     

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     

QM/MM Methods for Biomolecular Systems

Two are better than one : Quantum mechanics/molecular mechanics (QM/MM) methods are the state‐of‐the‐art computational technique for treating reactive and other “electronic” processes in biomolecular systems. This Review presents the general methodological aspects of the QM/MM approach, its use within optimization and simulation techniques, and its areas of application, always with a biomolecular focus.

     

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     

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.     

应用QM与ABEEM/MM结合的方法研究NAMI-A的结构性质

应用量子力学(QM)与ABEEM浮动电荷力场(ABEEM/MM)相结合的方法研究了抗癌药物NAMI-A在水溶液中的结构性质. 所有的结构优化都是在DFT的B3LYP方法下采用6-31G(d,p)和LanL2DZ基组完成的, 没有加入任何限制性条件. 结果表明, 优化得到的NAMI-A构型受不同环境及方法的影响均有变化. 与气相中得到的构型相比, QM/MM迭代优化得到构型要比PCM的构型变化更明显. QM/MM (ABEEM/MM)迭代优化得到的NAMI-A构型比QM/MM (OPLS-AA)的变化要小. 总之, 溶剂通过极化效应对NAMI-A结构、电荷分布及径向分布函数等性质均有影响, 客观地处理极化效应才能正确地反映QM区的性质.     

量子力学和分子力学组合方法及其应用

QM/MM组合方法在研究凝聚态中的化学反应及生物大分子的结构和活性之间的关系等方面已取得重要进展。这一方法的要点在于将大体系配分成几部分,根据需要对不同部分进行不同级别的处理,因此既利用了量子力学的精确性,又利用了分子力学的高效性。对QM/MM组合理论及其一些最新进展作一简单介绍,并以最近进行了几个工作为例说明QM、MM组合方法的应用。     

On the effect of a variation of the force field, spatial boundary condition and size of the QM region in QM/MM MD simulations

During the past years, the use of combined quantum-classical, QM/MM, methods for the study of complex biomolecular processes, such as enzymatic reactions and photocycles, has increased considerably. The quality of the results obtained from QM/MM calculations is largely dependent on five aspects to be considered when setting up a molecular model: the QM Hamiltonian, the MM Hamiltonian or force field, the boundary and coupling between the QM and MM regions, the size of the QM region and the boundary condition for the MM region. In this study, we systematically investigate the influence of a variation of the molecular mechanics force field and the size of the QM region in QM/MM MD simulations on properties of the photoactive part of the blue light photoreceptor protein AppA. For comparison, we additionally performed classical MD simulations and studied the effect of a variation of the type of spatial boundary condition. The classical boundary conditions and the force field used in a QM/MM MD simulation are shown to have non-neglegible effects upon the structural and energetic properties of the protein which makes it advisable to minimize computational artifacts in QM/MM MD simulations by application of periodic boundary conditions and a thermodynamically calibrated force field. A comparison of the structural and energetic properties of MD simulations starting from two alternative, different X-ray structures for the blue light utilizing flavin protein in its dark state indicates a slight preference of the two force fields used for the so-called Anderson structure over the Jung structure.     

Structure and dynamics of La(III) in aqueous solution – An ab initio QM/MM MD approach

Ab initio QM/MM MD simulations have allowed to clarify some of the ambiguities arising from various studies on the hydrated La(III) ion. Both nine- and ten-coordinated hydrates co-exist and interchange in a dissociative process on the nano- or even subnanosecond scale, and thus much faster than any other trivalent main group or transition metal ions. The weak ion–ligand bond (53 N/m) supplies a reasonable explanation for it. The simulation results for La(III) are also compared to those for the isoelectronic ions Cs(I) and Ba(II) obtained by the same ab initio MD procedure, leading to conclusions on the influence of central ion charge on structural and dynamic properties of hydrate complexes.     

Quantifying Electronic Effects in QM and QM/MM Biomolecular Modeling with the Fukui Function

Multi-scale quantum-mechanical/molecular-mechanical(QM/MM) and large-scale QM simulation provide valuable insight into enzyme mechanism and structure-property relationships. Analysis of the electron density afforded through these methods can enhance our understanding of how the enzyme environment modulates reactivity at the enzyme active site. From this perspective, tools from conceptual density functional theory to interrogate electron densities can provide added insight into enzyme function. We recently introduced the highly parallelizable Fukui shift analysis(FSA) method, which identifies how frontier states of an active site are altered by the presence of an additional QM residue to identify when QM treatment of a residue is essential as a result of quantum-mechanically affecting the behavior of the active site. We now demonstrate and analyze distance and residue dependence of Fukui function shifts in pairs of residues representing different non-covalent interactions. We also show how the interpretation of the Fukui function as a measure of relative nucleophilicity provides insight into enzymes that carry out S_N2 methyl transfer. The FSA method represents a promising approach for the systematic, unbiased determination of quantum mechanical effects in enzymes and for other complex systems that necessitate multi-scale modeling.     

The Jahn-Teller effect of the Cr(2+) ion in aqueous solution: Ab initio QM/MM molecular dynamics simulations

The hydration structure of Cr(2+) has been studied using molecular dynamics (MD) simulations including three-body corrections and combined ab initio quantum mechanical/molecular mechanical (QM/MM) MD simulations at the Hartree-Fock level. The structural properties are determined in terms of radial distribution functions, coordination numbers, and several angle distributions. The mean residence time was evaluated for describing ligand exchange processes in the second hydration shell. The Jahn-Teller distorted octahedral [Cr(H(2)O)(6)](2+) complex was pronounced in the QM/MM MD simulation. The first-shell distances of Cr(2+) are in the range of 1.9-2.8 A, which are slightly larger than those observed in the cases of Cu(2+) and Ti(3+). No first-shell water exchange occurred during the simulation time of 35 ps. Several water-exchange processes were observed in the second hydration shell with a mean residence time of 7.3 ps.     

QM/MM study of the taxadiene synthase mechanism

Combined quantum mechanics/molecular mechanics (QM/MM) calculations were used to investigate the reaction mechanism of taxadiene synthase (TXS). TXS catalyzes the cyclization of geranylgeranyl diphosphate (GGPP) to taxadiene (T) and four minor cyclic products. All these products originate from the deprotonation of carbocation intermediates. The reaction profiles for the conversion of GGPP to T as well as to minor products were calculated for different configurations of relevant TXS carbocation complexes. The QM region was treated at the M06-2X/TZVP level, while the CHARMM27 force field was used to describe the MM region. The QM/MM calculations suggest a reaction pathway for the conversion of GGPP to T, which slightly differs from previous proposals regarding the number of reaction steps and the conformation of the carbocations. The QM/MM results also indicate that the formation of minor products via water-assisted deprotonation of the carbocations is highly exothermic, by about −7 to −23 kcal/mol. Curiously, however, the computed barriers and reaction energies indicate that the formation of some of the minor products is more facile than the formation of T. Thus, the present QM/MM calculations provide detailed insights into possible reaction pathways and into the origin of the promiscuity of TXS, but they do not reproduce the product distribution observed experimentally. © 2019 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/MM and path integral simulations of kinetic isotope effects in the proton transfer reaction between nitroethane and acetate ion in water

An integrated Feynman path integral-free energy perturbation and umbrella sampling (PI-FEP/UM) method has been used to investigate the kinetic isotope effects (KIEs) in the proton transfer reaction between nitroethane and acetate ion in water. In the present study, both nuclear and electronic quantum effects are explicitly treated for the reacting system. The nuclear quantum effects are represented by bisection sampling centroid path integral simulations, while the potential energy surface is described by a combined quantum mechanical and molecular mechanical (QM/MM) potential. The accuracy essential for computing KIEs is achieved by a FEP technique that transforms the mass of a light isotope into a heavy one, which is equivalent to the perturbation of the coordinates for the path integral quasiparticle in the bisection sampling scheme. The PI-FEP/UM method is applied to the proton abstraction of nitroethane by acetate ion in water through molecular dynamics simulations. The rule of the geometric mean and the Swain-Schaad exponents for various isotopic substitutions at the primary and secondary sites have been examined. The computed total deuterium KIEs are in accord with experiments. It is found that the mixed isotopic Swain-Schaad exponents are very close to the semiclassical limits, suggesting that tunneling effects do not significantly affect this property for the reaction between nitroethane and acetate ion in aqueous solution.     

Origin of the solvent effects on the barrier to amide isomerization from the combined QM/MM Monte Carlo simulations

The origin of the solvent effects on the free energy of activation for the isomerization of N,N-dimethylformamide in water, CHCl₃ and CCl₄ has been investigated through statistical mechanical simulations using the combined quantum mechanical and molecular mechanical AM1/OPLS potential. The differential solvations between the ground state and transition state in various solvents can be attributed to the differences in molecular dipole moments in solution, and to the solvent polarization effects. In polar solvents, DMF is polarized more favorably in the ground state than in the rotated conformers, leading to greater solvent contributions. The modest solvent effects in CCl₄ are a reflection of its much smaller dielectric constant.     

How Chorismatases Regulate Distinct Reaction Channels in a Single Conserved Active Pocket: Mechanistic Analysis with QM/MM (ONIOM) Investigations

The FkbO and Hyg5 subfamilies of chorismatases share the same active-site architectures, but perform distinct reaction mechanisms, that is, FkbO employs a hydrolysis reaction whereas Hyg5 proceeds through an intramolecular mechanism. Despite extensive research efforts, the detailed mechanism of the product selectivity in chorismatases need to be further unmasked. In this study, the effects of the A/G residue group (A244^FkbO/G240^Hyg5) and the V/Q residue group (V209^FkbO/Q201^Hyg5) on the catalytic mechanisms are investigated by employing molecular dynamics simulations and hybrid quantum mechanical/molecular mechanical (QM/MM) calculations of the two wild-type models (FkbO/CHO and Hyg5/CHO; CHO=chorismate) and four mutants models (A244G-FkbO/CHO and G240A-Hyg5/CHO; V209Q-FkbO/CHO and Q201V-Hyg5/CHO). Our results showed that the A/G residue group mentioned by previous works would cause changes in the binding states of the substrate and the orientation of the catalytic glutamate, but only these changes affect the product selectivity in chorismatases limitedly. Interestingly, the distal V/Q residue group, which determines the internal water self-regulating ability at the active site, has significant impact on the selectivity of the catalytic mechanisms. The V/Q residue group is suggested to be an important factor to control the catalytic activities in chorismatases. The results are consistent with biochemical and structural experiments, providing novel insight into the mechanism of product selectivity in chorismatases.