Determinants of reactivity and selectivity in soluble epoxide hydrolase from quantum mechanics/molecular mechanics modeling

Soluble epoxide hydrolase (sEH) is an enzyme involved in drug metabolism that catalyzes the hydrolysis of epoxides to form their corresponding diols. sEH has a broad substrate range and shows high regio- and enantioselectivity for nucleophilic ring opening by Asp333. Epoxide hydrolases therefore have potential synthetic applications. We have used combined quantum mechanics/molecular mechanics (QM/MM) umbrella sampling molecular dynamics (MD) simulations (at the AM1/CHARMM22 level) and high-level ab initio (SCS-MP2) QM/MM calculations to analyze the reactions, and determinants of selectivity, for two substrates: trans-stilbene oxide (t-SO) and trans-diphenylpropene oxide (t-DPPO). The calculated free energy barriers from the QM/MM (AM1/CHARMM22) umbrella sampling MD simulations show a lower barrier for phenyl attack in t-DPPO, compared with that for benzylic attack, in agreement with experiment. Activation barriers in agreement with experimental rate constants are obtained only with the highest level of QM theory (SCS-MP2) used. Our results show that the selectivity of the ring-opening reaction is influenced by several factors, including proximity to the nucleophile, electronic stabilization of the transition state, and hydrogen bonding to two active site tyrosine residues. The protonation state of His523 during nucleophilic attack has also been investigated, and our results show that the protonated form is most consistent with experimental findings. The work presented here illustrates how determinants of selectivity can be identified from QM/MM simulations. These insights may also provide useful information for the design of novel catalysts for use in the synthesis of enantiopure compounds.     

COBRAMM 2.0 — A software interface for tailoring molecular electronic structure calculations and running nanoscale (QM/MM) simulations

We present a new version of the simulation software COBRAMM, a program package interfacing widely known commercial and academic software for molecular modeling. It allows a problem-driven tailoring of computational chemistry simulations with effortless ground and excited-state electronic structure computations. Calculations can be executed within a pure QM or combined quantum mechanical/molecular mechanical (QM/MM) framework, bridging from the atomistic to the nanoscale. The user can perform all necessary steps to simulate ground state and photoreactions in vacuum, complex biopolymer, or solvent environments. Starting from ground-state optimization, reaction path computations, initial conditions sampling, spectroscopy simulation, and photodynamics with deactivation events, COBRAMM is designed to assist in characterization and analysis of complex molecular materials and their properties. Interpretation of recorded spectra range from steady-state to time-resolved measurements. Various tools help the user to set up the system of interest and analyze the results.     

DRF90: a polarizable force field

The direct reaction field (DRF) approach has proven to be a useful tool to investigate the influence of solvents on the quantum/classical behaviour of solute molecules. In this paper, we report the latest extension of this DRF approach, which consists of the gradient of the completely classical energy expressions of this otherwise QM/MM method. They can be used in (completely classical) molecular dynamics (MD) simulations and geometry optimizations, that can be followed by a number of single point QM/MM calculations on configurations obtained in these simulations/optimizations. We report all energy and gradient expressions, and results for a number of interesting (model) systems. They include geometry optimization of the benzene dimer as well as MD simulations of some solvents. The most stable configuration for the benzene dimer is shown to be the parallel-displaced form, which is slightly more stable (0.3 kcal/mol) than the T-shaped dimer.     

The reaction mechanism of paraoxon hydrolysis by phosphotriesterase from combined QM/MM simulations

Molecular dynamics simulations employing combined quantum mechanical and molecular mechanical (QM/MM) potentials have been carried out to investigate the reaction mechanism of the hydrolysis of paraoxon by phosphotriesterase (PTE). We used a dual-level QM/MM approach that synthesizes accurate results from high-level electronic structure calculations with computational efficiency of semiempirical QM/MM potentials for free energy simulations. In particular, the intrinsic (gas-phase) energies of the active site in the QM region are determined by using density functional theory (B3LYP) and second-order M?ller-Plesset perturbation theory (MP2) and the molecular dynamics free energy simulations are performed by using the mixed AM1:CHARMM potential. The simulation results suggest a revised mechanism for the phosphotriester hydrolysis mechanism by PTE. The reaction free energy profile is mirrored by structural motions of the binuclear metal center in the active site. The two zinc ions occupy a compact conformation with an average zinc-zinc distance of 3.5 +/- 0.1 A in the Michaelis complex, whereas it is elongated to 5.3 +/- 0.3 A at the transition state and product state. The substrate is loosely bound to the more exposed zinc ion (Znbeta2+) at an average distance of 3.8 A +/- 0.3 A. The P=O bond of the substrate paraoxon is activated by adopting a tight coordination to the Znbeta2+, releasing the coordinate to the bridging hydroxide ion and increasing its nucleophilicity. It was also found that a water molecule enters into the binding pocket of the loosely bound binuclear center, originally occupied by the nucleophilic hydroxide ion. We suggest that the proton of this water molecule is taken up by His254 at low pH or released to the solvent at high pH, resulting in a hydroxide ion that pulls the Znbeta2+ ion closer to form the compact configuration and restores the resting state of the enzyme.     

QM/MM free energy simulations: recent progress and challenges

Due to the higher computational cost relative to pure molecular mechanical (MM) simulations, hybrid quantum mechanical/molecular mechanical (QM/MM) free energy simulations particularly require a careful consideration of balancing computational cost and accuracy. Here, we review several recent developments in free energy methods most relevant to QM/MM simulations and discuss several topics motivated by these developments using simple but informative examples that involve processes in water. For chemical reactions, we highlight the value of invoking enhanced sampling technique (e.g. replica-exchange) in umbrella sampling calculations and the value of including collective environmental variables (e.g. hydration level) in metadynamics simulations; we also illustrate the sensitivity of string calculations, especially free energy along the path, to various parameters in the computation. Alchemical free energy simulations with a specific thermodynamic cycle are used to probe the effect of including the first solvation shell into the QM region when computing solvation free energies. For cases where high-level QM/MM potential functions are needed, we analyse two different approaches: the QM/MM–MFEP method of Yang and co-workers and perturbative correction to low-level QM/MM free energy results. For the examples analysed here, both approaches seem productive although care needs to be exercised when analysing the perturbative corrections.     

Modeling reaction routes from rhodopsin to bathorhodopsin

The quantum mechanical–molecular mechanical (QM/MM) theory was applied to calculate accurate structural parameters, vibrational and optical spectra of bathorhodopsin (BATHO), one of the primary photoproducts of the functional cycle of the visual pigment rhodopsin (RHO), and to characterize reaction routes from RHO to BATHO. The recently resolved crystal structure of BATHO (PDBID: 2G87) served as an initial source of coordinates of heavy atoms. Protein structures in the ground electronic state and vibrational frequencies were determined by using the density functional theory in the PBE0/cc‐pVDZ approximation for the QM part and the AMBER force field parameters in the MM part. Calculated and assigned vibrational spectra of both model protein systems, BATHO and RHO, cover three main regions referring to the hydrogen‐out‐of‐plan (HOOP) motion, the C?C ethylenic stretches, and the C? C single‐bond stretches. The S₀–S₁ electronic excitation energies of the QM part, including the chromophore group in the field of the protein matrix, were estimated by using the advanced quantum chemistry methods. The computed structural parameters as well as the spectral bands match perfectly the experimental findings. A structure of the transition state on the S₀ potential energy surface for the ground electronic state rearrangement from RHO to BATHO was located proving a possible route of the thermal protein activation to the primary photoproduct. Proteins 2010. © 2009 Wiley‐Liss, Inc.     

Hybrid quantum mechanical/molecular mechanical molecular dynamics simulations of HIV-1 integrase/inhibitor complexes

Human immunodeficiency virus (HIV)-1 integrase (IN) is an attractive target for development of acquired immunodeficiency syndrome chemotherapy. In this study, conventional and coupled quantum mechanical and molecular mechanical (QM/MM) molecular dynamics (MD) simulations of HIV-1 IN complexed with 5CITEP (IN-5CITEP) were carried out. In addition to differences in the bound position of 5CITEP, significant differences at the two levels of theory were observed in the metal coordination geometry and the areas involving residues 116-119 and 140-166. In the conventional MD simulation, the coordination of Mg(2+) was found to be a near-perfect octahedral geometry whereas a distorted octahedral complex was observed in QM/MM. All of the above reasons lead to a different pattern of protein-ligand salt link formation that was not observed in the classical MD simulation. Furthermore to provide a theoretical understanding of inhibition mechanisms of 5CITEP and its derivative (DKA), hybrid QM/MM MD simulations of the two complexes (IN-5CITEP and IN-DKA) have been performed. The results reveal that areas involving residues 60-68, 116-119, and 140-149 were substantially different among the two systems. The two systems show similar pattern of metal coordination geometry, i.e., a distorted octahedron. In IN-DKA, both OD1 and OD2 of Asp-64 coordinate the Mg(2+) in a monodentate fashion whereas only OD1 is chelated to the metal as observed in IN-5CITEP. The high potency of DKA as compared to 5CITEP is supported by a strong salt link formed between its carboxylate moiety and the ammonium group of Lys-159. Detailed comparisons between HIV-1 IN complexed with DKA and with 5CITEP provide information about ligand structure effects on protein-ligand interactions in particular with the Lys-159. This is useful for the design of new selective HIV-1 IN inhibitors.     

A computational study of the deacylation mechanism of human butyrylcholinesterase

To investigate the mechanism of the deacylation reaction in the active site of human butyrylcholinesterase (BuChE), we carried out quantum mechanical (QM) calculations on cluster models of the active site built from a crystallographic structure. The models consisted of the substrate butyrate moiety, the catalytic triad of residues (Ser198, Glu325, and His438), the "oxy-anion hole" (Gly116, Gly117, and Ala199), the side chain of Glu197, four water molecules, the side chain of Ser225, and the peptide linkage between Val321 and Asn322. Analyses of the equilibrium geometries, electronic properties, and energies of the QM models gave insights into the catalytic mechanism. In addition, the QM calculations provided the data required to build a molecular mechanics representation of the reactive BuChE region that was employed in molecular dynamics simulations followed by molecular-mechanics-Poisson-Boltzmann (MM-PB) calculations. Subsequently, we combined the QM energies with average MM-PB energies to estimate the free energy of the reactive structures in the enzyme. The rate-determining step corresponds to the formation of a tetrahedral intermediate with a free-energy barrier of approximately 14.0 kcal/mol. The modulation of the BuChE activity, exerted by either neutral molecules (glycerol, GOL) or a second butyrylcholine (CHO) molecule bound to the cation-pi site, does not involve any significant allosteric effect. Interestingly, the presence of GOL or CHO stabilizes a product complex formed between a butyric acid molecule and BuChE. These results are in consonance with the crystallographic structure of BuChE, in which the catalytic Ser198 interacts with a butyric fragment, while the cation-pi site is occupied by one GOL molecule.     

Protein effects in non-heme iron enzyme catalysis: insights from multiscale models

Many non-heme iron enzymes have similar sets of ligands but still catalyze widely different reactions. A key question is, therefore, the role of the protein in controlling reactivity and selectivity. Examples from multiscale simulations, primarily QM/MM, of both mono- and binuclear non-heme iron enzymes are used to analyze the stability of these models and what they reveal about the protein effects. Consistent results from QM/MM modeling are the importance of the hydrogen bond network to control reactivity and electrostatic stabilization of electron transfer from second-sphere residues. The long-range electrostatic effects on reaction barriers are small for many systems. In the systems where large electrostatic effects have been reported, these lead to higher barriers. There is thus no evidence of any significant long-range electrostatic effects contributing to the catalytic efficiency of non-heme iron enzymes. However, the correct evaluation of electrostatic contributions is challenging, and the correlation between calculated residue contributions and the effects of mutation experiments is not very strong. The largest benefits of QM/MM models are thus the improved active-site geometries, rather than the calculation of accurate energies. Reported differences in mechanistic predictions between QM and QM/MM models can be explained by differences in hydrogen bonding patterns in and around the active site. Correctly constructed cluster models can give results with similar accuracy as those from multiscale models, but the latter reduces the risk of drawing the wrong mechanistic conclusions based on incorrect geometries and are preferable for all types of modeling, even when using very large QM parts.     

Recent developments in QM/MM methods towards open-boundary multi-scale simulations

Combined quantum mechanics/molecular mechanics (QM/MM) methods have been widely used in multi-scale modelling and simulations of physical, chemical and biological processes in complex environments. In this review, we provide an overview of the recently developed QM/MM algorithms, with emphasis on our works, towards the ultimate goal of establishing an open boundary between the QM and MM subsystems. The open boundary is characterised by on-the-fly exchanges of partial charges and atoms between the QM and MM subsystems, allowing us to focus on the small QM subsystem of primary interest in dynamics simulations. An open-boundary scheme has the promise to the utilisations of small QM subsystems, high-levels of QM theory and long simulation times, which can potentially lead to new insights.     

Multiscale methods for macromolecular simulations

In this article we review the key modeling tools available for simulating biomolecular systems. We consider recent developments and representative applications of mixed quantum mechanics/molecular mechanics (QM/MM), elastic network models (ENMs), coarse-grained molecular dynamics, and grid-based tools for calculating interactions between essentially rigid protein assemblies. We consider how the different length scales can be coupled, both in a sequential fashion (e.g. a coarse-grained or grid model using parameterization from MD simulations), and via concurrent approaches, where the calculations are performed together and together control the progression of the simulation. We suggest how the concurrent coupling approach familiar in the context of QM/MM calculations can be generalized, and describe how this has been done in the CHARMM macromolecular simulation package.     

Refined models of New Delhi metallo-beta-lactamase-1 with inhibitors: an QM/MM modeling study

New Delhi metallo-beta-lactamase 1 (NDM-1) has been identified as a potential target for the treatment of multi-drug resistance bacterial infections. We used molecular docking, normal MD, SIE, QM/MM MD simulations, QM/MM GBSA binding free energy, and QM/MM GBSA alanine-scanning mutagenesis techniques to investigate interactions of the NDM-1 with 11 inhibitors (Tigecycline, BAL30072, D-captopril, Penicillin G, Ampicillin, Carbenicillin, Cephalexin, Cefaclor, Nitrocefin, Meropenem, and Imipenem). From our normal MD and QM/MM simulations, the correlation coefficients between the predicted binding free energies and experimental values are .88 and .93, respectively. Then simulations, which combined QM/MM/GBSA and alanine-scanning mutagenesis techniques, were performed and our results show that two residues (Lys211 and His250) have the strongest impact on the binding affinities of the 11 NDM-1/inhibitors. Therefore, our approach theoretically suggests that the two residues (Lys211 and His250) are responsible for the selectivity of NDM-1 associated inhibitors.     

How does the protein environment optimize the thermodynamics of thiol sulfenylation? Insights from model systems to QM/MM calculations on human 2-Cys peroxiredoxin

Protein thiol/sulfenic acid oxidation potentials provide a tool to select specific oxidation agents, but are experimentally difficult to obtain. Here, insights into the thiol sulfenylation thermodynamics are obtained from model calculations on small systems and from a quantum mechanics/molecular mechanics (QM/MM) analysis on human 2-Cys peroxiredoxin thioredoxin peroxidase B (Tpx-B). To study thiol sulfenylation in Tpx-B, our recently developed computational method to determine reduction potentials relatively compared to a reference system and based on reaction energies reduction potential from electronic energies is updated. Tpx-B forms a sulfenic acid (R-SO^?) on one of its active site cysteines during reactive oxygen scavenging. The observed effect of the conserved active site residues is consistent with the observed hydrogen bond interactions in the QM/MM optimized Tpx-B structures and with free energy calculations on small model systems. The ligand effect could be linked to the complexation energies of ligand L with CH₃S^? and CH₃SO^?. Compared to QM only calculations on Tpx-B’s active site, the QM/MM calculations give an improved understanding of sulfenylation thermodynamics by showing that other residues from the protein environment other than the active site residues can play an important role.     

Theoretical approaches for dynamical ordering of biomolecular systems

Background

Living systems are characterized by the dynamic assembly and disassembly of biomolecules. The dynamical ordering mechanism of these biomolecules has been investigated both experimentally and theoretically. The main theoretical approaches include quantum mechanical (QM) calculation, all-atom (AA) modeling, and coarse-grained (CG) modeling. The selected approach depends on the size of the target system (which differs among electrons, atoms, molecules, and molecular assemblies). These hierarchal approaches can be combined with molecular dynamics (MD) simulation and/or integral equation theories for liquids, which cover all size hierarchies.

Correlation between biological activity and binding energy in systems of integrin with cyclic RGD-containing binders: a QM/MM molecular dynamics study

We here report a combined quantum mechanical/molecular mechanical (QM/MM) molecular dynamics (MD) study on the binding interactions between the α(V)β(3) integrin and eight cyclic arginine-glycine-aspartate (RGD) containing peptides. The initial conformation of each peptide within the binding site of the integrin was determined by docking the ligand to the reactive site of the integrin crystal structure with the aid of docking software FRED. The subsequent QM/MM MD simulations of the complex structures show that these eight cyclic RGD-peptides have a generally similar interaction mode with the binding site of the integrin to the cyclo(RGDf-N[M]V) analog found in the crystal structure. Still, there are subtle differences in the interactions of peptide ligands with the integrin, which contribute to the different inhibition activities. The averaged QM/MM protein-ligand interaction energy (IE) is remarkably correlated to the biological activity of the ligand. The IE, as well as a three-variable model which is somewhat interpretable, thus can be used to predict the bioactivity of a new ligand quantitatively, at least within a family of analogs. The present study establishes a helpful protocol for advancing lead compounds to potent inhibitors.     

Photochemical Reaction Dynamics of the Primary Event of Vision Studied by Means of a Hybrid Molecular Simulation

The photoisomerization reaction dynamics of a retinal chromophore in the visual receptor rhodopsin was investigated by means of hybrid quantum mechanical/molecular mechanical (QM/MM) molecular dynamics (MD) simulations. The photoisomerization reaction of retinal constitutes the primary step of vision and is known as one of the fastest reactions in nature. To elucidate the molecular mechanism of the high efficiency of the reaction, we carried out hybrid ab initio QM/MM MD simulations of the complete reaction process from the vertically excited state to the photoproduct via electronic transition in the entire chromophore-protein complex. An ensemble of reaction trajectories reveal that the excited-state dynamics is dynamically homogeneous and synchronous even in the presence of thermal fluctuation of the protein, giving rise to the very fast formation of the photoproduct. The synchronous nature of the reaction dynamics in rhodopsin is found to originate from weak perturbation of the protein surroundings and from dynamic regulation of volume-conserving motions of the chromophore. The simulations also provide a detailed view of time-dependent modulations of hydrogen-out-of-plane vibrations during the reaction process, and identify molecular motions underlying the experimentally observed dynamic spectral modulations.     

Numerical studies of the membrane fluorescent dyes dynamics in ground and excited states

Fluorescence methods are widely used in studies of biological and model membranes. The dynamics of membrane fluorescent markers in their ground and excited electronic states and correlations with their molecular surrounding within the fully hydrated phospholipid bilayer are still not well understood. In the present work, Quantum Mechanical (QM) calculations and Molecular Dynamics (MD) simulations are used to characterize location and interactions of two membrane polarity probes (Prodan; 6-propionyl-2-dimethylaminonaphthalene and its derivative Laurdan; 2-dimethylamino-6-lauroylnaphthalene) with the dioleoylphosphatidylcholine (DOPC) lipid bilayer model. MD simulations with fluorophores in ground and excited states are found to be a useful tool to analyze the fluorescent dye dynamics and their immediate vicinity. The results of QM calculations and MD simulations are in excellent agreement with available experimental data. The calculation shows that the two amphiphilic dyes initially placed in bulk water diffuse within 10 ns towards their final location in the lipid bilayer. Analysis of solvent relaxation process in the aqueous phase occurs on the picoseconds timescale whereas it takes nanoseconds at the lipid/water interface. Four different relaxation time constants, corresponding to different relaxation processes, where observed when the dyes were embedded into the membrane.     

Combined quantum and molecular mechanics calculations on metalloproteins

The combination of quantum mechanics and molecular mechanics (QM/MM) methods is one of the most promising approaches to study the structure, function and properties of proteins. The number of QM/MM applications on metalloproteins is steadily increasing, especially studies with density functional methods on redox-active metal centres. Recent developments include new parameterised methods to treat covalent bonds between the quantum and classical systems, methods to obtain free energy from QM/MM results, and the combination of quantum chemistry and protein crystallography.