A Quantum Mechanical and Molecular Mechanical Method Based on CM1A Charges: Applications to Solvent Effects on Organic Equilibria and Reactions

Elucidation of Hydrolysis Mechanisms for Fatty Acid Amide Hydrolase and Its Lys142Ala Variant via QM/MM Simulations

Fatty acid amide hydrolase (FAAH) is a serine hydrolase that degrades anandamide, an endocannabinoid, and oleamide, a sleep-inducing lipid, and has potential applications as a therapeutic target for neurological disorders. Remarkably, FAAH hydrolyzes amides and esters with similar rates; however, the normal preference for esters reemerges when Lys142 is mutated to alanine. To elucidate the hydrolysis mechanisms and the causes behind this variation of selectivity, mixed quantum and molecular mechanics (QM/MM) calculations were carried out to obtain free-energy profiles for alternative mechanisms for the enzymatic hydrolyses. The methodology features free-energy perturbation calculations in Monte Carlo simulations with PDDG/PM3 as the QM method. For wild-type FAAH, the results support a mechanism, which features proton transfer from Ser217 to Lys142, simultaneous proton transfer from Ser241 to Ser217, and attack of Ser241 on the substrate's carbonyl carbon to yield a tetrahedral intermediate, which subsequently undergoes elimination with simultaneous protonation of the leaving group by a Lys142-Ser217 proton shuttle. For the Lys142Ala mutant, a striking multistep sequence is proposed with simultaneous proton transfer from Ser241 to Ser217, attack of Ser241 on the carbonyl carbon of the substrate, and elimination of the leaving group and its protonation by Ser217. Support comes from the free-energy results, which well reproduce the observation that the Lys142Ala mutation in FAAH decreases the rate of hydrolysis for oleamide significantly more than for methyl oleate.

Conformational Equilibrium of 1,2-Dichloroethane in Water: Comparison of PCM and RISM-SCF Methods

The RISM-SCF and polarizable continuum model (PCM) approaches have been applied to study the conformational equilibrium of 1,2-dichloroethane (DCE) in water. Both the electron correlation effect and basis sets play an important role in the relative energies of the gauche and trans conformers in gas and solution phases. Both PCM and RISM-MP2 methods resulted in a consistent trend with the previous experimental and theoretical studies that the population of the gauche conformer increases in going from the gas phase to the aqueous solution. However, the PCM treatment could not describe the solvent effect completely in that the sign of the relative free energy of the gauche and trans forms is opposite to the most recent experimental and theoretical data, while the RISM-MP2 gives the right sign in the free energy difference. We found that the larger excess chemical potential gain (by ca. -4.1 kcal/mol) for the gauche conformer is large enough to result in the gauche preference of DCE in water, though it has to compensate for more solute reorganization energy (approximately 1.6 kcal/mol) and overcome the energy difference (approximately 1.6 kcal/mol) in the gas phase. The radial distribution functions between DCE and the nearest water shows that the electrostatic repulsion between chlorine and oxygen atoms is higher in the trans conformer than in the gauche one, while the attractive interaction between chlorine and hydrogen of water is higher in the gauche conformer.

Cope Elimination: Elucidation of Solvent Effects from QM/MM Simulations

The Cope elimination reactions for threo- and erythro-N,N-dimethyl-3-phenyl-2-butylamine oxide have been investigated using QM/MM calculations in water, THF, and DMSO. The aprotic solvents provide up to million-fold rate accelerations. The effects of solvation on the reactants, transition structures, and rates of reaction are elucidated here using two-dimensional potentials of mean force (PMF) derived from free-energy perturbation calculations in Monte Carlo simulations (MC/FEP). The resultant free energies of activation in solution are in close agreement with experiment. Ab initio calculations at the MP2/6-311+G-(2d,p) level using the PCM continuum solvent model were also carried out; however, only the QM/MM methodology was able to reproduce the large rate increases in proceeding from water to the dipolar aprotic solvents. Solute-solvent interaction energies and radial distribution functions are also analyzed and show that poorer solvation of the reactant in the aprotic solvents is primarily responsible for the observed rate enhancements. It is found that the amine oxide oxygen is the acceptor of three hydrogen bonds from water molecules for the reactant but only one to two weaker ones at the transition state. The overall quantitative success of the computations supports the present QM/MM/MC approach, featuring PDDG/PM3 as the QM method.

Molecular modeling of organic and biomolecular systems using BOSS and MCPRO

An overview is provided of the capabilities for the current versions of the BOSS and MCPRO programs for molecular modeling of organic and biomolecular systems. Recent applications are noted, particularly for QM/MM studies of organic and enzymatic reactions and for protein-ligand binding.

Solvent effects and mechanism for a nucleophilic aromatic substitution from QM/MM simulations

The nucleophilic aromatic substitution (SNAr) reaction between azide ion and 4-fluoronitrobenzene has been investigated using QM/MM and DFT/PCM calculations in protic and dipolar aprotic solvents. The effects of solvation on the transition structures, the intermediate Meisenheimer complex, and the rate of reaction are elucidated. The large rate increases in proceeding from protic to dipolar aprotic solvents are only reproduced by the QM/MM methodology.

Influence of inter- and intramolecular hydrogen bonding on kemp decarboxylations from QM/MM simulations

The Kemp decarboxylation reaction for benzisoxazole-3-carboxylic acid derivatives has been investigated using QM/MM calculations in protic and dipolar aprotic solvents. Aprotic solvents have been shown to accelerate the rates of reaction by 7-8 orders of magnitude over water; however, the inclusion of an internal hydrogen bond effectively inhibits the reaction with near solvent independence. The effects of solvation and intramolecular hydrogen bonding on the reactants, transition structures, and the rate of reaction are elucidated using two-dimensional potentials of mean force (PMF) derived from free energy perturbation calculations in Monte Carlo simulations (MC/FEP). Free energies of activation in six solvents have been computed to be in close agreement with experiment. Solute-solvent interaction energies show that poorer solvation of the reactant anion in the dipolar aprotic solvents is primarily responsible for the observed rate enhancements over protic media. In addition, a discrepancy for the experimental rate in chloroform has been studied in detail with the conclusion that ion-pairing between the reactant anion and tetramethylguanidinium counterion is responsible for the anomalously slow reaction rate. The overall quantitative success of the computations supports the present QM/MM/MC approach, which features PDDG/PM3 as the QM method.

Potential energy functions for atomic-level simulations of water and organic and biomolecular systems

An overview is provided on the development and status of potential energy functions that are used in atomic-level statistical mechanics and molecular dynamics simulations of water and of organic and biomolecular systems. Some topics that are considered are the form of force fields, their parameterization and performance, simulations of organic liquids, computation of free energies of hydration, universal extension for organic molecules, and choice of atomic charges. The discussion of water models covers some history, performance issues, and special topics such as nuclear quantum effects.

An ab initio and cryospectroscopic study of the hydrogen chloride and boron trifluoride complexes of cyclopropene

Mid-infrared spectra of cyclopropene (c-C3H4) mixtures with HCl or BF3, dissolved in liquefied argon and in liquefied nitrogen have been examined. Evidence was found for the formation of two different isomers of the 1:1 complexes: a pi-type complex in which the interaction occurs with the pi-bond of cyclopropene, and a sigma-type complex in which the electron acceptors interact with one of the CC single bonds. At lower temperatures, indications for the formation of a 1:2 chain-type complex were found. Using spectra recorded between 90 and 124 K in liquid argon, the complexation enthalpies for the pi- and sigma-complexes with HCl were determined to be -8.8(3) and -7.9(3) kJ mol-1, respectively. For the pi-complex with BF3 a value of -7.4(3)kJ mol-1 was found. Structural and spectral information on the 1:1 complexes was obtained from B3LYP/6-311++G(d,p) and MP2(FC)/6-31+G(d) calculations. Using Free Energy Perturbation Monte Carlo simulations to calculate the solvent influences, and statistical thermodynamics to account for zero-point vibrational and thermal contributions, the solution enthalpies were transformed into complexation energies with values of -15.5(6) and -13.8(8) kJ mol-1 for the pi- and sigma-isomers of c-C3H4.HCl, respectively, and of -12.7(8) kJ mol-1 for the pi-type BF3 complex.

Macrophomate Synthase: QM/MM Simulations Address the Diels−Alder versus Michael−Aldol Reaction Mechanism

Macrophomate synthase (MPS) of the phytopathogenic fungus Macrophoma commelinae catalyzes the transformation of 2-pyrone derivatives to the corresponding benzoate analogues. The transformation proceeds through three separate chemical reactions, including decarboxylation of oxalacetate to produce pyruvate enolate, two C-C bond formations between 2-pyrone and pyruvate enolate that form a bicyclic intermediate, and final decarboxylation with concomitant dehydration. Although some evidence suggests that the second step of the reaction catalyzed by MPS is a Diels-Alder reaction, definite proof that the C-C bond formations occur via a concerted mechanism is still required. An alternative route for formation of the C-C bonds is a stepwise Michael-aldol reaction. In this work, mixed quantum and molecular mechanics (QM/MM) combined with Monte Carlo simulations and free-energy perturbation (FEP) calculations were performed to investigate the relative stabilities of the transition states (TS) for both reaction mechanisms. The key results are that the Diels-Alder TS model is 17.7 and 12.1 kcal/mol less stable than the Michael and aldol TSs in the stepwise route, respectively. Therefore, this work indicates that the Michael-aldol mechanism is the route used by MPS to catalyze the second step of the overall transformation, and that the enzyme is not a natural Diels-Alderase, as claimed by Ose and co-workers (Nature 2003, 422, 185-189; Acta Crystallogr. 2004, D60, 1187-1197). A modified link-atom treatment for the bonds at the QM/MM interface is also presented.

Steric retardation of SN2 reactions in the gas phase and solution

The gas-phase S(N)2 reactions of chloride with ethyl and neopentyl chlorides and their alpha-cyano derivatives have been explored with B3LYP, CBS-QB3, and PDDG/PM3 calculations. Calculations predict that the steric effect of the tert-butyl group raises the activation energy by about 6 kcal/mol relative to methyl in both cases. Solvent effects have been computed with QM/MM Monte Carlo simulations for DMSO, methanol, and water, as well as with a polarizable continuum model, CPCM. Solvents cause a large increase in the activation energies of these reactions but have a very small differential effect on the ethyl and neopentyl substrates and their cyano derivatives. The theoretical results contrast with previous conclusions that were based upon gas-phase rate measurements.

Quantum descriptors for biological macromolecules from linear-scaling electronic structure methods

The characterization of electrostatic and chemical properties at the surface of biological macromolecules is of interest in elucidating the fundamental biological structure-function relationships as well as in problems of rational drug design. This paper presents a set of macromolecular quantum descriptors for the characterization of biological macromolecules in solution that can be obtained with modest computational cost from linear-scaling semi-empirical quantum/solvation methods. The descriptors discussed include: solvent-polarized electrostatic surface potential maps, equilibrated atomic charges, Fukui reactivity indices, approximate local hardness maps, and relative proton potentials. These properties are applied to study the conformational dependence of the electrostatic surface potential of the solvated phosphate-binding protein mutant (T141D), the regioselectivity of the zinc finger domains of HIV-1 nucleocapsid (NC) protein, and the order of pKa values of acidic residues in turkey ovomucoid third domain (OMTKY3) and of the zinc-binding residues in the carboxyl terminal zinc finger of NC. In all cases, insight beyond that obtainable from purely classical models is gained and can be used to rationalize the experimental observations. The macromolecular quantum descriptors presented here greatly extend the arsenal of tools for macromolecular characterization and offer promise in applications to modern structure-based drug design.

Extension of the PDDG/PM3 and PDDG/MNDO semiempirical molecular orbital methods to the halogens

The new semiempirical methods, PDDG/PM3 and PDDG/MNDO, have been parameterized for halogens. For comparison, the original MNDO and PM3 were also reoptimized for the halogens using the same training set; these modified methods are referred to as MNDO' and PM3'. For 442 halogen-containing molecules, the smallest mean absolute error (MAE) in heats of formation is obtained with PDDG/PM3 (5.6 kcal/mol), followed by PM3' (6.1 kcal/mol), PDDG/MNDO (6.6 kcal/mol), PM3 (8.1 kcal/mol), MNDO' (8.5 kcal/mol), AM1 (11.1 kcal/mol), and MNDO (14.0 kcal/mol). For normal-valent halogen-containing molecules, the PDDG methods also provide improved heats of formation over MNDO/d. Hypervalent compounds were not included in the training set and improvements over the standard NDDO methods with sp basis sets were not obtained. For small haloalkanes, the PDDG methods yield more accurate heats of formation than are obtained from density functional theory (DFT) with the B3LYP and B3PW91 functionals using large basis sets. PDDG/PM3 and PM3' also give improved binding energies over the standard NDDO methods for complexes involving halide anions, and they are competitive with B3LYP/6-311++G(d,p) results including thermal corrections. Among the semiempirical methods studied, PDDG/PM3 also generates the best agreement with high-level ab initio G2 and CCSD(T) intrinsic activation energies for S(N)2 reactions involving methyl halides and halide anions. Finally, the MAEs in ionization potentials, dipole moments, and molecular geometries show that the parameter sets for the PDDG and reoptimized NDDO methods reduce the MAEs in heats of formation without compromising the other important QM observables.

Accuracy of free energies of hydration using CM1 and CM3 atomic charges

Absolute free energies of hydration (DeltaGhyd) have been computed for 25 diverse organic molecules using partial atomic charges derived from AM1 and PM3 wave functions via the CM1 and CM3 procedures of Cramer, Truhlar, and coworkers. Comparisons are made with results using charges fit to the electrostatic potential surface (EPS) from ab initio 6-31G* wave functions and from the OPLS-AA force field. OPLS Lennard-Jones parameters for the organic molecules were used together with the TIP4P water model in Monte Carlo simulations with free energy perturbation theory. Absolute free energies of hydration were computed for OPLS united-atom and all-atom methane by annihilating the solutes in water and in the gas phase, and absolute DeltaGhyd values for all other molecules were computed via transformation to one of these references. Optimal charge scaling factors were determined by minimizing the unsigned average error between experimental and calculated hydration free energies. The PM3-based charge models do not lead to lower average errors than obtained with the EPS charges for the subset of 13 molecules in the original study. However, improvement is obtained by scaling the CM1A partial charges by 1.14 and the CM3A charges by 1.15, which leads to average errors of 1.0 and 1.1 kcal/mol for the full set of 25 molecules. The scaled CM1A charges also yield the best results for the hydration of amides including the E/Z free-energy difference for N-methylacetamide in water.

Contributions of conformational compression and preferential transition state stabilization to the rate enhancement by chorismate mutase

The rate enhancement provided by the chorismate mutase (CM) enzyme for the Claisen rearrangement of chorismate to prephenate has been investigated by application of the concept of near attack conformations (NACs). Using a combined QM/MM Monte Carlo/free-energy perturbation (MC/FEP) method, 82% and 100% of chorismate conformers were found to be NAC structures in water and in the CM active site, respectively. Consequently, the conversion of non-NACs to NACs does not contribute to the free energy of activation from preorganization of the substrate into NACs. The FEP calculations yielded differences in free energies of activation that well reproduce the experimental data. Additional calculations indicate that the rate enhancement by CM over the aqueous phase results primarily from conformational compression of NACs by the enzyme and that this process is enthalpically controlled. This suggests that preferential stabilization of the transition state in the enzyme environment relative to water plays a secondary role in the catalysis by CM.

Investigation of solvent effects for the Claisen rearrangement of chorismate to prephenate: mechanistic interpretation via near attack conformations

Solvent effects on the rate of the Claisen rearrangement of chorismate to prephenate have been examined in water and methanol. The preequilibrium free-energy differences between diaxial and diequatorial conformers of chorismate, which had previously been implicated as the sole basis for the observed 100-fold rate increase in water over methanol, have been reframed using the near attack conformation (NAC) concept of Bruice and co-workers. Using a combined QM/MM Monte Carlo/free-energy perturbation (MC/FEP) method, 82%, 57%, and 1% of chorismate conformers were found to be NAC structures (NACs) in water, methanol, and the gas phase, respectively. As a consequence, the conversion of non-NACs to NACs provides no free-energy contributions to the overall relative reaction rates in water versus methanol. Free-energy perturbation calculations yielded differences in free energies of activation for the two polar protic solvents and the gas phase. The rate enhancement in water over the gas phase arises from preferential hydration of the transition state (TS) relative to the reactants via increased hydrogen bonding and long-range electrostatic interactions, which accompany bringing the two negatively charged carboxylates into closer proximity. More specifically, there is an increase of 1.3 and 0.6 hydrogen bonds to the carboxylate groups and the ether oxygen, respectively, in going from the reactant to the TS in water. In methanol, the corresponding changes in hydrogen bonding with first shell solvent molecules are small; the rate enhancement arises primarily from the enhanced long-range interactions with solvent molecules. Thus, the reaction occurs faster in water than in methanol due to greater stabilization of the TS in water by specific interactions with first shell solvent molecules.

Comparison of formation of reactive conformers for the SN2 displacements by CH3CO2- in water and by Asp124-CO2- in a haloalkane dehalogenase

The S(N)2 displacement of Cl(-) from 1,2-dichloroethane by acetate (CH(3)CO(2)(-)) in water and by the carboxylate of the active site aspartate in the haloalkane dehalogenase of Xanthobacter autothropicus have been compared by using molecular dynamics simulations. In aqueous solution, six families of contact-pair structures (I-VI) were identified, and their relative concentrations and dissociation rate constants were determined. The near attack conformers (NACs) required for the S(N)2 displacement reaction are members of the IV (CH(3)COO(-)...CH(2)(Cl)CH(2)Cl) family and are formed in the sequence II-->III-->IV-->NAC. The NAC subclass is defined by the COO(-)...CCl contact distance of < or = 3.41 A and the COO(-)...CCl angle of 157-180 degrees. The mole percentage of NACs is 0.16%, based on the 1 M standard state. This result may be compared with 13.4 mole percentage of NACs in the Michaelis complex in the enzyme. It follows that NAC formation in the enzyme is favored by 2.6 kcal/mol. Because reaction coordinates from S to TS, both in water and in the enzyme, pass via NAC (i.e., S --> NAC --> TS), the reduction in the S --> NAC barrier by 2.6 kcal/mol accounts for approximately 25% of the reduction of total barrier in the S --> TS (10.7 kcal/mol). The remaining 75% of the advantage of the enzymatic reaction revolves around the efficiency of NAC --> TS step. This process, based on previous studies, is discussed briefly.

Fluoride-selective binding in a new deep cavity calix[4]pyrrole: experiment and theory

A new "super-extended cavity" tetraacetylcalix[4]pyrrole derivative was synthesized and characterized, and X-ray crystal structures of complexes bound to fluoride and acetonitrile were obtained. The binding behavior of this receptor was investigated by NMR titration, and the complex was found to exclusively bind fluoride ions in DMSO-d(6). This unusual binding behavior was investigated by Monte Carlo free energy perturbation simulations and Poisson calculations, and the ion specificity was seen to result from the favorable electrostatic interactions that the fluoride gains by sitting lower in the phenolic cavity of the receptor. The effect of water present in the DMSO on the calculated free energies of binding was also investigated. Owing to the use of a saturated ion solution, the effect of contaminating water is small in this case; however, it has the potential to be very significant at lower ion concentrations. Finally, the adaptive umbrella WHAM protocol was investigated and optimized for use in binding free energy calculations, and its efficiency was compared to that of the free energy perturbation calculations; adaptive umbrella WHAM was found to be approximately two times more efficient. In addition, structural evidence demonstrates that the protocol explores a wider conformational range than free energy perturbation and should therefore be the method of choice. This paper represents the first complete application of this methodology to "alchemical" changes.

Prediction of drug solubility from structure

The aqueous solubility of a drug is an important factor affecting its bioavailability. Numerous computational methods have been developed for the prediction of aqueous solubility from a compound's structure. A review is provided of the methodology and quality of results for the most useful procedures including the model implemented in the QikProp program. Viable methods now exist for predictions with less than 1 log unit uncertainty, which is adequate for prescreening synthetic candidates or design of combinatorial libraries. Further progress with predictive methods would require an experimental database of highly accurate solubilities for a large, diverse collection of drug-like molecules.

Influence of intermolecular hydrogen bonds on the tautomerism of pyridine derivatives

The effect of the dimerization, by hydrogen-bond (HB) complexation, on the tautomerism of 2-hydroxypyridine and a series 2-aminopyridines has been carried using ab initio methods. The results obtained for 2-hydroxypyridine fit satisfactorily with the experimental data and show that the 2-pyridone/2-pyridone homodimer complex is the most stable. For 2-aminopyridines, the effect of the substituent on the amino group has been investigated. For the monomers studied, the most favorable tautomer is the 2H; however, with electronegative substituents, the 1H/1H homodimers are more stable than the corresponding 2H/2H ones. The atom in molecule methodology has been used to characterize the HBs formed. Exponential relationships have been found between the electron density and its laplacian at the HB critical point vs the HB distance.

Rationale for the observed COX-2/COX-1 selectivity of celecoxib from Monte Carlo simulations

Computational studies have yielded an analysis of the contributions to the free energy difference between the binding of celecoxib to COX-1 and to COX-2. The energetic and structural results point to the Ile to Val mutation at residue 523 as the key contributor to COX-2 selectivity; unfavorable steric contact between a sulfonamide oxygen and the delta methyl group of Ile523 destabilizes the complex with COX-1. The His to Arg change at residue 513 is less significant.