Systematic QM/MM investigation of factors that affect the cytochrome P450-catalyzed hydrogen abstraction of camphor

The hydrogen abstraction reaction of camphor in cytochrome P450(cam) has been investigated in the native enzyme environment by combined quantum mechanical/molecular mechanical (QM/MM) calculations and in the gas phase by density functional calculations. This work has been motivated by contradictory published QM/MM results. In an attempt to pinpoint the origin of these discrepancies, we have systematically studied the factors that may affect the computed barriers, including the QM/MM setup, the optimization procedures, and the choice of QM region, basis set, and protonation states. It is found that the ChemShell and QSite programs used in the published QM/MM calculations yield similar results at given geometries, and that the discrepancies mainly arise from two technical issues (optimization protocols and initial system preparation) that need to be well controlled in QM/MM work. In the course of these systematic investigations, new mechanistic insights have been gained. The crystallographic water 903 placed near the oxo atom of Compound I lowers the hydrogen abstraction barrier by ca. 4 kcal/mol, and thus acts as a catalyst for this reaction. Spin density may appear at the A-propionate side chain of the heme if the carboxylate group is not properly screened, which might be expected to happen during protein dynamics, but not in static equilibrium situations. There is no clear correlation between the computed A-propionate spin density and the hydrogen abstraction barrier, and hence, no support for a previously proposed side-chain mediated transition state stabilization mechanism. Standard QM/MM optimizations yield an A-propionate environment close to the X-ray structure only for protonated Asp297, and not for deprotonated Asp297, but the computed barriers are similar in both cases. An X-ray like A-propionate environment can also be obtained when deprotonated Asp297 is included in the QM region and His355 is singly protonated, but this Compound II-type species with a closed-shell porphyrin ring has a higher hydrogen abstraction barrier and should thus not be mechanistically relevant.     

Quantum mechanical/molecular mechanical simulations of the Tl(III) ion in water

Classical molecular dynamics (MD) and combined quantum mechanical/molecular mechanical (QM/MM) MD simulations have been performed to investigate the structural and dynamical properties of the Tl(III) ion in water. A six-coordinate hydration structure with a maximum probability of the Tl-O distance at 2.21 A was observed, which is in good agreement with X-ray data. The librational and vibrational spectra of water molecules in the first hydration shell are blue-shifted compared with those of pure liquid water, and the Tl-O stretching force constant was evaluated as 148 Nm(-1). Both structural and dynamical properties show a distortion of the first solvation shell structure. The second shell ligands' mean residence time was determined as 12.8 ps. The Tl(III) ion can be classified as "structure forming" ion; the calculated hydration energy of -986 +/- 9 kcal mol agrees well with the experimental value of -986 kcal mol.     

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 mechanisms for compound I formation in the catalytic cycle of cytochrome P450cam

In the catalytic cycle of cytochrome P450cam, after molecular oxygen binds as a ligand to the heme iron atom to yield a ferrous dioxygen complex, there are fast proton transfers that lead to the formation of the active species, Compound I (Cpd I), which are not well understood because they occur so rapidly. In the present work, the conversion of the ferric hydroperoxo complex (Cpd 0) to Cpd I has been investigated by combined quantum-mechanical/molecular-mechanical (QM/MM) calculations. The residues Asp(251) and Glu(366) are considered as proton sources. In mechanism I, a proton is transported to the distal oxygen atom of the hydroperoxo group via a hydrogen bonding network to form protonated Cpd 0 (prot-Cpd0: FeOOH(2)), followed by heterolytic O-O bond cleavage that generates Cpd I and water. Although a local minimum is found for prot-Cpd0 in the Glu(366) channel, it is very high in energy (more than 20 kcal/mol above Cpd 0) and the barriers for its decay are only 3-4 kcal/mol (both toward Cpd 0 and Cpd I). In mechanism II, an initial O-O bond cleavage followed by a concomitant proton and electron transfer yields Cpd I and water. The rate-limiting step in mechanism II is O-O cleavage with a barrier of about 13-14 kcal/mol. According to the QM/MM calculations, the favored low-energy pathway to Cpd I is provided by mechanism II in the Asp(251) channel. Cpd 0 and Cpd I are of similar energies, with a slight preference for Cpd I.     

Classical and quantum mechanical/molecular mechanical molecular dynamics simulations of alanine dipeptide in water: comparisons with IR and vibrational circular dichroism spectra

We have implemented the combined quantum mechanical (QM)/molecular mechanical (MM) molecular dynamics (MD) simulations of alanine dipeptide in water along with the polarizable and nonpolarizable classical MD simulations with different models of water. For the QM/MM MD simulation, the alanine dipeptide is treated with the AM1 or PM3 approximations and the fluctuating solute dipole moment is calculated by the Mulliken population analysis. For the classical MD simulations, the solute is treated with the polarizable or nonpolarizable AMBER and polarizable CHARMM force fields and water is treated with the TIP3P, TIP4P, or TIP5P model. It is found that the relative populations of right-handed alpha-helix and extended beta and P(II) conformations in the simulation trajectory strongly depend on the simulation method. For the QM/MM MD simulations, the PM3/MM shows that the P(II) conformation is dominant, whereas the AM1/MM predicts that the dominant conformation is alpha(R). Polarizable CHARMM force field gives almost exclusively P(II) conformation and other force fields predict that both alpha-helical and extended (beta and P(II)) conformations are populated with varying extents. Solvation environment around the dipeptide is investigated by examining the radial distribution functions and numbers and lifetimes of hydrogen bonds. Comparing the simulated IR and vibrational circular dichroism spectra with experimental results, we concluded that the dipeptide adopts the P(II) conformation and PM3/MM, AMBER03 with TIP4P water, and AMBER polarizable force fields are acceptable for structure determination of the dipeptide considered in this paper.     

Structure and dynamics of high-spin Ru in aqueous solution: Ab initio QM/MM molecular dynamics simulation

The structural and dynamical properties of high-spin Ru²⁺ in aqueous solution have been theoretically studied using molecular dynamics (MD) simulations. The conventional MD simulation based on pair potentials gives the overestimated average first shell coordination number of 9, whereas the value of 5.9 was observed when the three-body corrected function was included. A combined ab initio quantum mechanical/molecular mechanical (QM/MM) molecular dynamics simulation has been performed to take into account the many-body effects on the hydration shell structure of Ru²⁺. The most important region, the first hydration shell, was treated by ab initio quantum mechanics at UHF level using the SBKJC VDZ ECP basis set for Ru²⁺ and the 6-31G^∗ basis sets for water. An exact coordination number of 6 for the first hydration shell was obtained from the QM/MM simulation. The QM/MM simulation predicts the average Ru²⁺–O distance of 2.42 Å for the first hydration shell, whereas the values of 2.34 and 2.46 Å are resulted from the pair potentials without and with the three-body corrected simulations, respectively. Several other structural properties representing position and orientation of the solvate molecules were evaluated for describing the hydration shell structure of the Ru²⁺ ion in dilute aqueous solution. A mean residence time of 7.1 ps was obtained for water ligands residing in the second hydration shell.     

Development of massive multilevel molecular dynamics simulation program,platypus (PLATform for dYnamic protein unified simulation), for the elucidation of protein functions

A massively parallel program for quantum mechanical‐molecular mechanical (QM/MM) molecular dynamics simulation, called Platypus (PLATform for dYnamic Protein Unified Simulation), was developed to elucidate protein functions. The speedup and the parallelization ratio of Platypus in the QM and QM/MM calculations were assessed for a bacteriochlorophyll dimer in the photosynthetic reaction center (DIMER) on the K computer, a massively parallel computer achieving 10 PetaFLOPs with 705,024 cores. Platypus exhibited the increase in speedup up to 20,000 core processors at the HF/cc‐pVDZ and B3LYP/cc‐pVDZ, and up to 10,000 core processors by the CASCI(16,16)/6‐31G** calculations. We also performed excited QM/MM‐MD simulations on the chromophore of Sirius (SIRIUS) in water. Sirius is a pH‐insensitive and photo‐stable ultramarine fluorescent protein. Platypus accelerated on‐the‐fly excited‐state QM/MM‐MD simulations for SIRIUS in water, using over 4000 core processors. In addition, it also succeeded in 50‐ps (200,000‐step) on‐the‐fly excited‐state QM/MM‐MD simulations for the SIRIUS in water. © 2016 The Authors. Journal of Computational Chemistry Published by Wiley Periodicals, Inc.     

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.     

Quantum and molecular mechanical study of the first proton transfer in the catalytic cycle of cytochrome P450cam and its mutant D251N

In the catalytic cycle of cytochrome P450cam, the hydroperoxo intermediate (Cpd 0) is formed by proton transfer from a reduced oxyheme complex (S5). This process is drastically slowed down when Asp251 is mutated to Asn (D251N). We report quantum mechanical/molecular mechanical (QM/MM) calculations that address this proton delivery in the doublet state through a hydrogen-bond network in the Asp251 channel, both for the wild-type enzyme and the D251N mutant, using four different active-site models. For the wild-type, we find a facile concerted mechanism for proton transfer from protonated Asp251 via Wat901 and Thr252 to the FeOO moiety, with a barrier of about 1 kcal/mol and a high exothermicity of more than 20 kcal/mol. In the D251N mutant with a neutral Asn251 residue, the proton transfer is almost thermoneutral or slightly exothermic in the three models considered. It is still very facile when the Asn251 residue adopts a conformation analogous to Asp251 in the wild-type enzyme, but the barrier increases significantly when the Asn251 side chain flips (as indicated by classical molecular dynamics simulations). This flip disrupts the hydrogen-bond network and hence the proton-transfer pathway, which causes a longer lifetime of S5 in the D251N mutant (consistent with experimental observations). The entry of an additional water molecule into the active site of D251N with flipped Asn251 regenerates the hydrogen-bond network and provides a viable mechanism for proton delivery in the mutant, with a moderate barrier of about 7 kcal/mol.     

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.     

Probing the transition from hydrophilic to hydrophobic solvation with atomic scale resolution

Picosecond and femtosecond X-ray absorption spectroscopy is used to probe the changes of the solvent shell structure upon electron abstraction of aqueous iodide using an ultrashort laser pulse. The transient L(1,3) edge EXAFS at 50 ps time delay points to the formation of an expanded water cavity around the iodine atom, in good agreement with classical and quantum mechanical/molecular mechanics (QM/MM) molecular dynamics (MD) simulations. These also show that while the hydrogen atoms pointed toward iodide, they predominantly point toward the bulk solvent in the case of iodine, suggesting a hydrophobic behavior. This is further confirmed by quantum chemical (QC) calculations of I(-)/I(0)(H(2)O)(n=1-4) clusters. The L(1) edge sub-picosecond spectra point to the existence of a transient species that is not present at 50 ps. The QC calculations and the QM/MM MD simulations identify this transient species as an I(0)(OH(2)) complex inside the cavity. The simulations show that upon electron abstraction most of the water molecules move away from iodine, while one comes closer to form the complex that lives for 3-4 ps. This time is governed by the reorganization of the main solvation shell, basically the time it takes for the water molecules to reform an H-bond network. Only then is the interaction with the solvation shell strong enough to pull the water molecule of the complex toward the bulk solvent. Overall, much of the behavior at early times is determined by the reorientational dynamics of water molecules and the formation of a complete network of hydrogen bonded molecules in the first solvation shell.     

Combined quantum mechanical/molecular mechanical study on the pentacoordinated ferric and ferrous cytochrome P450cam complexes

The pentacoordinated ferric and ferrous cytochrome P450(cam) complexes have been investigated by combined quantum mechanical/molecular mechanical (QM/MM) calculations in the presence of a protein/solvent environment and by QM calculations on the isolated QM regions with use of density functional theory. The B3LYP functional has been found more reliable than the BLYP and BHLYP functionals for estimating the relative state energies. The B3LYP/CHARMM calculations with an all-electron basis set for iron give high-spin ground states for the title complexes, in agreement with experiment. The comparison of the B3LYP/CHARMM results of the entire protein system with the B3LYP calculations on the naked QM regions shows that the amount of stabilization by the protein environment is largest for the intermediate-spin states, followed by the high-spin states of the complexes. The calculation of M?ssbauer parameters in the presence of the enzyme environment confirms the double occupation of the d(xz) orbital in the quintet spin state of the ferrous complex, consistent with the computed QM/MM energies in the enzyme environment, while the d(x)2(-)(y)2 orbital is doubly occupied in the gas-phase quintet state.     

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.     

The elusive oxidant species of cytochrome P450 enzymes: characterization by combined quantum mechanical/molecular mechanical (QM/MM) calculations

The primary oxidant of cytochrome P450 enzymes, Compound I, is hard to detect experimentally; in the case of cytochrome P450(cam), this intermediate does not accumulate in solution during the catalytic cycle even at temperatures as low as 200 K (ref 4). Theory can play an important role in characterizing such elusive species. We present here combined quantum mechanical/molecular mechanical (QM/MM) calculations of Compound I of cytochrome P450(cam) in the full enzyme environment as well as density functional studies of the isolated QM region. The calculations assign the ground state of the species, quantify the effect of polarization and hydrogen bonding on its properties, and show that the protein environment and its specific hydrogen bonding to the cysteinate ligand are crucial for sustaining the Fe-S bond and for preventing the full oxidation of the sulfur.     

Structure-breaking effects of solvated Rb(I) in dilute aqueous solution--an ab initio QM/MM MD approach

Structural properties of the hydrated Rb(I) ion have been investigated by ab initio quantum mechanical/molecular mechanical (QM/MM) molecular dynamics (MD) simulations at the double-zeta HF quantum mechanical level. The first shell coordination number was found to be 7.1, and several other structural parameters such as angular distribution functions, radial distribution functions and tilt- and theta-angle distributions allowed the full characterization of the hydration structure of the Rb(I) ion in dilute aqueous solution. Velocity autocorrelation functions were used to calculate librational and vibrational motions, ion-ligand motions, as well as reorientation times. Different dynamical parameters such as water reorientation, mean ligand residence time, the number of ligand exchange processes, and rate constants were also analyzed. The mean ligand residence time for the first shell was determined as tau = 2.0 ps.     

A general acid-base mechanism for the stabilization of a tetrahedral adduct in a serine-carboxyl peptidase: a computational study

The QM/MM MD and free energy simulations show that serine-carboxyl peptidases (sedolisins) may stabilize the tetrahedral intermediates and tetrahedral adducts primarily through a general acid-base mechanism involving Asp (Asp164 for kumamolisin-As) rather than the oxyanion-hole interactions as in the cases of serine proteases.     

An application of the novel quantum mechanical/molecular mechanical method combined with the theory of energy representation: An ionic dissociation of a water molecule in the supercritical water

A novel quantum chemical approach recently developed has been applied to an ionic dissociation of a water molecule (2H(2)O-->H(3)O(+)+OH(-)) in ambient and supercritical water. The method is based on the quantum mechanical/molecular mechanical (QM/MM) simulations combined with the theory of energy representation (QM/MM-ER), where the energy distribution function of MM solvent molecules around a QM solute serves as a fundamental variable to determine the hydration free energy of the solute according to the rigorous framework of the theory of energy representation. The density dependence of the dissociation free energy in the supercritical water has been investigated for the density range from 0.1 to 0.6 g/cm(3) with the temperature fixed at a constant. It has been found that the product ionic species significantly stabilizes in the high density region as compared with the low density. Consequently, the dissociation free energy decreases monotonically as the density increases. The decomposition of the hydration free energy has revealed that the entropic term (-TDeltaS) strongly depends on the density of the solution and dominates the behavior of the dissociation free energy with respect to the variation of the density. The increase in the entropic term in the low density region can be attributed to the decrease in the translational degrees of freedom brought about by the aggregation of solvent water molecules around the ionic solute.     

Inhibitor binding by metallo-beta-lactamase IMP-1 from Pseudomonas aeruginosa: quantum mechanical/molecular mechanical simulations

The dynamics of the IMP-1 enzyme complexed with three prototypical inhibitors are investigated using a quantum mechanical/molecular mechanical (QM/MM) method based on the self-consistent-charge density-functional tight-binding model. The binding patterns of the inhibitors observed in X-ray diffraction experiments are well reproduced in 600 ps molecular dynamics simulations at room temperature. These inhibitors anchor themselves in the enzyme active site by direct coordination with the two zinc ions, displacing the hydroxide nucleophile that bridges the two zinc ions. In addition, they also interact with several active-site residues and those in two mobile loops. The excellent agreement with experimental structural data validates the QM/MM treatment used in our simulations.     

Performance assessment of semiempirical molecular orbital methods in describing halogen bonding: quantum mechanical and quantum mechanical/molecular mechanical-molecular dynamics study

The performance of semiempirical molecular-orbital methods--MNDO, MNDO-d, AM1, RM1, PM3 and PM6--in describing halogen bonding was evaluated, and the results were compared with molecular mechanical (MM) and quantum mechanical (QM) data. Three types of performance were assessed: (1) geometrical optimizations and binding energy calculations for 27 halogen-containing molecules complexed with various Lewis bases (Two of the tested methods, AM1 and RM1, gave results that agree with the QM data.); (2) charge distribution calculations for halobenzene molecules, determined by calculating the solvation free energies of the molecules relative to benzene in explicit and implicit generalized Born (GB) solvents (None of the methods gave results that agree with the experimental data.); and (3) appropriateness of the semiempirical methods in the hybrid quantum-mechanical/molecular-mechanical (QM/MM) scheme, investigated by studying the molecular inhibition of CK2 protein by eight halobenzimidazole and -benzotriazole derivatives using hybrid QM/MM molecular-dynamics (MD) simulations with the inhibitor described at the QM level by the AM1 method and the rest of the system described at the MM level. The pure MM approach with inclusion of an extra point of positive charge on the halogen atom approach gave better results than the hybrid QM/MM approach involving the AM1 method. Also, in comparison with the pure MM-GBSA (generalized Born surface area) binding energies and experimental data, the calculated QM/MM-GBSA binding energies of the inhibitors were improved by replacing the G(GB,QM/MM) solvation term with the corresponding G(GB,MM) term.