Reaction field effects on nitrogen shielding

Modeling of Solute-Solvent Interactions Using an External Electric Field-From Tautomeric Equilibrium in Nonpolar Solvents to the Dissociation of Alkali Metal Halides

An implicit account of the solvent effect can be carried out using traditional static quantum chemistry calculations by applying an external electric field to the studied molecular system. This approach allows one to distinguish between the effects of the macroscopic reaction field of the solvent and specific solute-solvent interactions. In this study, we report on the dependence of the simulation results on the use of the polarizable continuum approximation and on the importance of the solvent effect in nonpolar solvents. The latter was demonstrated using experimental data on tautomeric equilibria between the pyridone and hydroxypyridine forms of 2,6-di-tert-butyl-4-hydroxy-pyridine in cyclohexane and chloroform.

A Benefit of Using the IDSCRF- over UFF-Radii Cavities and Why Joint Correlations of NMR Chemical Shifts Can Be Advantageous: Condensed Pyridines as an IEF-PCM/GIAO/DFT Case Study

Herein, an advantage of the use of IDSCRF- over UFF-radii-based solute cavities in GIAO/DFT calculations is presented for the ¹³C and especially ¹⁵N NMR chemical shifts made for several bicyclic aromatic nitrogen heterocycles in CDCl₃ solution treated within the classical IEF-PCM solvation scheme. Successful application of the IDSCRF-radii in the non 1:1 joint multinuclear ¹H/¹³C and particularly ¹H/¹³C/¹⁵N correlations of the measured δ_H,C(,N) values to those obtained theoretically is also documented for a series of test systems (-268 ≤ δ_N ≤ -72 ppm). The experimentally yet unknown δ_N's were found in this way for the title compounds via a trinuclear eq 1 determined for an optimally chosen value of the multiplication factor of initial raw δ_H data (m_H = 10). Such a simultaneous analysis of the δ_H,C(,N) data is proposed as a novel method to study the solution structure of the other similar conformationally homogeneous (bio)organic compounds. The issue of small spurious imaginary vibrational frequencies computed for a few molecular systems using the Gaussian 09 default UFF-radii is briefly considered as well.

Solvent Dependence of (14)N Nuclear Magnetic Resonance Chemical Shielding Constants as a Test of the Accuracy of the Computed Polarization of Solute Electron Densities by the Solvent

Although continuum solvation models have now been shown to provide good quantitative accuracy for calculating free energies of solvation, questions remain about the accuracy of the perturbed solute electron densities and properties computed from them. Here we examine those questions by applying the SM8, SM8AD, SMD, and IEF-PCM continuum solvation models in combination with the M06-L density functional to compute the (14)N magnetic resonance nuclear shieldings of CH3CN, CH3NO2, CH3NCS, and CH3ONO2 in multiple solvents, and we analyze the dependence of the chemical shifts on solvent dielectric constant. We examine the dependence of the computed chemical shifts on the definition of the molecular cavity (both united-atom models and models based on superposed individual atomic spheres) and three kinds of treatments of the electrostatics, namely the generalized Born approximation with the Coulomb field approximation, the generalized Born model with asymmetric descreening, and models based on approximate numerical solution schemes for the nonhomogeneous Poisson equation. Our most systematic analyses are based on the computation of relative (14)N chemical shifts in a series of solvents, and we compare calculated shielding constants relative to those in CCl4 for various solvation models and density functionals. While differences in the overall results are found to be reasonably small for different solvation models and functionals, the SMx models SM8, and SM8AD, using the same cavity definitions (which for these models means the same atomic radii) as those employed for the calculation of free energies of solvation, exhibit the best agreement with experiment for every functional tested. This suggests that in addition to predicting accurate free energies of solvation, the SM8 and SM8AD generalized Born models also describe the solute polarization in a manner reasonably consistent with experimental (14)N nuclear magnetic resonance spectroscopy. Models based on the nonhomogeneous Poisson equation show slightly reduced accuracy. Scaling the intrinsic Coulomb radii to larger values (as has sometimes been suggested in the past) does not uniformly improve the results for any kind of solvent model; furthermore it uniformly degrades the results for generalized Born models. Use of a basis set that increases the outlying charge diminishes the accuracy of continuum models that solve the nonhomogeneous Poisson equation, which we ascribe to the inability of the numerical schemes for approximately solving the nonhomogeneous Poisson equation to fully account for the effects of electronic charge outside the solute cavity.

Reaction pathways and free energy profiles for spontaneous hydrolysis of urea and tetramethylurea: unexpected substituent effects

It has been difficult to directly measure the spontaneous hydrolysis rate of urea and, thus, 1,1,3,3-tetramethylurea (Me4U) was used as a model to determine the "experimental" rate constant for urea hydrolysis. The use of Me4U was based on an assumption that the rate of urea hydrolysis should be 2.8 times that of Me4U hydrolysis because the rate of acetamide hydrolysis is 2.8 times that of N,N-dimethyl-acetamide hydrolysis. The present first-principles electronic-structure calculations on the competing non-enzymatic hydrolysis pathways have demonstrated that the dominant pathway is the neutral hydrolysis via the CN addition for both urea (when pH < ~11.6) and Me4U (regardless of pH), unlike the non-enzymatic hydrolysis of amides where alkaline hydrolysis is dominant. Based on the computational data, the substituent shift of the free energy barrier calculated for the neutral hydrolysis is remarkably different from that for the alkaline hydrolysis, and the rate constant for the urea hydrolysis should be ~1.3 × 10(9)-fold lower than that (4.2 × 10(-12) s(-1)) measured for the Me4U hydrolysis. As a result, the rate enhancement and catalytic proficiency of urease should be 1.2 × 10(25) and 3 × 10(27) M(-1), respectively, suggesting that urease surpasses proteases and all other enzymes in its power to enhance the rate of reaction. All of the computational results are consistent with available experimental data for Me4U, suggesting that the computational prediction for urea is reliable.

Modeling in vitro inhibition of butyrylcholinesterase using molecular docking, multi-linear regression and artificial neural network approaches

Butyrylcholinesterase (BChE) has been an important protein used for development of anti-cocaine medication. Through computational design, BChE mutants with ∼2000-fold improved catalytic efficiency against cocaine have been discovered in our lab. To study drug-enzyme interaction it is important to build mathematical model to predict molecular inhibitory activity against BChE. This report presents a neural network (NN) QSAR study, compared with multi-linear regression (MLR) and molecular docking, on a set of 93 small molecules that act as inhibitors of BChE by use of the inhibitory activities (pIC₅₀ values) of the molecules as target values. The statistical results for the linear model built from docking generated energy descriptors were: r(2)=0.67, rmsd=0.87, q(2)=0.65 and loormsd=0.90; the statistical results for the ligand-based MLR model were: r(2)=0.89, rmsd=0.51, q(2)=0.85 and loormsd=0.58; the statistical results for the ligand-based NN model were the best: r(2)=0.95, rmsd=0.33, q(2)=0.90 and loormsd=0.48, demonstrating that the NN is powerful in analysis of a set of complicated data. As BChE is also an established drug target to develop new treatment for Alzheimer's disease (AD). The developed QSAR models provide tools for rationalizing identification of potential BChE inhibitors or selection of compounds for synthesis in the discovery of novel effective inhibitors of BChE in the future.

A density functional and quantum Monte Carlo study of glutamic acid in vacuo and in a dielectric continuum medium

We present density functional theory (DFT) and quantum Monte Carlo (QMC) calculations of the glutamic acid and glutamate ion in vacuo and in various dielectric continuum media within the polarizable continuum model (PCM). In DFT, we employ the integral equation formalism variant of PCM while, in QMC, we use a PCM scheme we have developed to include both surface and volume polarization. We investigate the gas-phase protonation thermochemistry of the glutamic acid using a large set of structural conformations, and find that QMC is in excellent agreement with the best available theoretical and experimental results. For the solvated glutamic acid and glutamate ion, we perform DFT calculations for dielectric constants, ε, between 4 and 78. We find that the glutamate ion in the zwitterionic form is more stable than the non-zwitterionic form over the whole range of dielectric constants, while the glutamic acid is more stable in its non-zwitterionic form at ε = 4. The dielectric constant at which the two glutamic acid species have the same energy depends on the cavity size and lies between 5 and 12.5. We validate these results with QMC for the two limiting values of the dielectric constant, and find qualitative agreement with DFT even though the solvent polarization is less pronounced at the QMC level.

Thermodynamics of chemical reactions with COSMO-RS: the extreme case of charge separation or recombination

Many technically relevant chemical processes in the condensed phase involve as elementary reactive steps the formation of ions from neutral species or, as the opposite, recombination of ions. Such reactions that generate or annihilate charge defy the standard gas phase quantum chemical treatment, and also continuum solvation models are only partially able to account for the right amount of stabilization in solution. In this work, for such types of reaction, a solvation treatment involving the COSMO-RS method is assessed, which leads to improved results, i.e., errors of only around 10 kJ/mol for both protic and aprotic solvents. The examples discussed here comprise protolysis reactions and organo halide heterolysis, for both of which a comparison with reliable experimental data is possible. It is observed that for protolysis, the quality of results does not strongly depend on the quantum chemical method used for energy calculation. In contrast, in the case of heterolytic carbon-chlorine bond cleavage, clearly better results are obtained for higher correlated (coupled cluster) methods or the density functional M06-2X, which is well known for its accuracy if applied to organic chemistry. This hints at least that the right answer is obtained for the right reason and not due to a compensation of errors from gas phase thermodynamics with those from the solvation treatment. Problems encountered with certain critical solvents or upon decomposing Gibbs free energies into heats or entropies of reaction are found to relate mostly to the parameterization of the H-bonding term within COSMO-RS.

Catalytic mechanism of cytochrome P450 for 5'-hydroxylation of nicotine: fundamental reaction pathways and stereoselectivity

A series of computational methods were used to study how cytochrome P450 2A6 (CYP2A6) interacts with (S)-(-)-nicotine, demonstrating that the dominant molecular species of (S)-(-)-nicotine in CYP2A6 active site exists in the free base state (with two conformations, SR(t) and SR(c)), despite the fact that the protonated state is dominant for the free ligand in solution. The computational results reveal that the dominant pathway of nicotine metabolism in CYP2A6 is through nicotine free base oxidation. Further, first-principles quantum mechanical/molecular mechanical free energy (QM/MM-FE) calculations were carried out to uncover the detailed reaction pathways for the CYP2A6-catalyzed nicotine 5'-hydroxylation reaction. In the determined CYP2A6-(S)-(-)-nicotine binding structures, the oxygen of Compound I (Cpd I) can abstract a hydrogen from either the trans-5'- or the cis-5'-position of (S)-(-)-nicotine. CYP2A6-catalyzed (S)-(-)-nicotine 5'-hydroxylation consists of two reaction steps, that is, the hydrogen transfer from the 5'-position of (S)-(-)-nicotine to the oxygen of Cpd I (the H-transfer step), followed by the recombination of the (S)-(-)-nicotine moiety with the iron-bound hydroxyl group to generate the 5'-hydroxynicotine product (the O-rebound step). The H-transfer step is rate-determining. The 5'-hydroxylation proceeds mainly with the stereoselective loss of the trans-5'-hydrogen, that is, the 5'-hydrogen trans to the pyridine ring. The calculated overall stereoselectivity of ∼97% favoring the trans-5'-hydroxylation is close to the observed stereoselectivity of 89-94%. This is the first time it has been demonstrated that a CYP substrate exists dominantly in one protonation state (cationic species) in solution, but uses its less-favorable protonation state (neutral free base) to perform the enzymatic reaction.

Structural assignment of 6-oxy purine derivatives through computational modeling, synthesis, X-ray diffraction, and spectroscopic analysis

6-Oxy purine derivatives have been considered as potential therapeutic agents in various drug discovery efforts reported in the literature. However, the structural assignment of this important class of compounds has been controversial concerning the specific position of a hydrogen atom in the structure. To theoretically determine the most favorable type of tautomeric form of 6-oxy purine derivatives, we have carried out first-principles electronic structure calculations on the possible tautomeric forms (A, B, and C) and their relative stability of four representative 6-oxy purine derivatives (compounds 1-4). The computational results in both the gas phase and aqueous solution clearly reveal that the most favorable type of tautomeric form of these compounds is A, in which a hydrogen atom bonds with the N1 atom on the purine ring. To examine the computational results, one of the 6-oxy purine derivatives (i.e., compound 4) has been synthesized and its structure has been characterized by X-ray diffraction and spectroscopic analysis. All of the obtained computational and experimental data are consistent with the conclusion that the 6-oxy purine derivative exists in tautomer A. The conclusive structural assignment reported here is expected to be valuable for future computational studies on 6-oxy purine derivative binding with proteins and for computational drug design involving this type of compounds.

Free Energies of Solvation with Surface, Volume, and Local Electrostatic Effects and Atomic Surface Tensions to Represent the First Solvation Shell

Building on the SVPE (surface and volume polarization for electrostatics) model for electrostatic contributions to the free energy of solvation with explicit consideration of both surface and volume polarization effects, on the SMx approach to including first-solvation-shell contributions, and on the linear relationship between the electric field and short-range electrostatic contributions found by Chipman, we have developed a new method for computing absolute aqueous solvation free energies by combining the SVPE method with semiempirical terms that account for effects beyond bulk electrostatics. The new method is called SMVLE, and the elements it contains are denoted by SVPE-CDSL where SVPE denotes accounting for bulk electrostatic interactions between solute and solvent with both surface and volume contributions, CDS denotes the inclusion of solvent cavitation, changes in dispersion energy, and possible changes in local solvent structure by a semiempirical term utilizing geometry-dependent atomic surface tensions as implemented in SMx models, and L represents the local electrostatic effect derived from the outward-directed normal electric field on the cavity surface. The semiempirical CDS and L terms together represent the deviation of short-range contributions to the free energy of solvation from those accounted for by the SVPE term based on the bulk solvent dielectric constant. A solute training set containing a broad range of molecules used previously in the development of SM6 is used here for SMVLE model calibration. The aqueous solvation free energies predicted by the parameterized SMVLE model correlate exceedingly well with experimental values. The square of the correlation coefficient is 0.9949 and the slope is 1.0079. Comparison of the final SMVLE model against the earlier SMx solvation model shows that the parameterized SMVLE model not only yields good accuracy for neutrals but also significantly increases the accuracy for ions, making it the best implicit solvation model to date for aqueous solvation free energies of ions. The semiempirical terms associated with the outward-directed electric field account in a physical way for the improvement in the predictive accuracy for ions. The SMVLE method greatly decreases the need to include explicit water molecules for accurate modeling of solvation free energies of ions.

First-Principles Determination of Molecular Conformations of Indolizidine (-)-235B' in Solution

Indolizidine (-)-235B' is a particularly interesting natural product, as it is the currently known, most potent and subtype-selective open-channel blocker of the alpha4beta2 nicotinic acetylcholine receptor (nAChR). In the current study, extensive first-principles electronic structure calculations have been carried out in order to determine the stable molecular conformations and their relative free energies of the protonated and deprotonated states of (-)-235B' in the gas phase, in chloroform, and in aqueous solution. The (1)H and (13)C NMR chemical shifts calculated using the computationally determined dominant molecular conformation of the deprotonated state are all consistent with available experimental NMR spectra of (-)-235B' in chloroform, which suggests that the computationally determined molecular conformations are reasonable. Our computational results reveal for the first time that two geminal H atoms on carbon-3 (C3) of (-)-235B' have remarkably different chemical shifts (i.e. 3.24 and 2.03 ppm). The computational results help one to better understand and analyze the experimental (1)H NMR spectra of (-)-235B'. The finding of remarkably different chemical shifts of two C3 geminal H atoms in a certain molecular conformation of (-)-235B' may also be valuable in analysis of NMR spectra of other related ring-containing compounds. In addition, the pK(a) of (-)-235B' in aqueous solution is predicted to be ~9.7. All of the computational results provide a solid basis for future studies of the microscopic and phenomenological binding of various receptor proteins with the protonated and deprotonated structures of this unique open-channel blocker of alpha4beta2 nAChRs. This computational study also demonstrates how one can appropriately use computational modeling and spectroscopic analysis to address the structural and spectroscopic problems that cannot be addressed by experiments alone.

First-principles determination of molecular conformations of cyclic adenosine 3',5'-monophosphate in gas phase and aqueous solution

Extensive first-principles electronic structure calculations were performed in this study to explore the possible molecular structures and their concentration distribution of an intracellular second messenger, that is, cyclic adenosine 3',5'-monophosphate (cAMP), and its protonated form (cAMPH) in the gas phase and aqueous solution. The calculations resulted in prediction of four different stable conformers of cAMP and eight different stable conformers of cAMPH and their relative Gibbs free energies in the gas phase and aqueous solution. All of the computational results consistently demonstrate that the predominant conformers of cAMP and cAMPH are always the cAMP-chair-anti and cAMPH-chair2-syn conformers, respectively, in both the gas phase and aqueous solution. It has been demonstrated that the free energy barriers calculated for the intertransformation reactions between different conformers are very low (below approximately 6 kcal/mol) such that the intertransformation reactions between different conformers are very fast so that the concentration distribution of the system can quickly reach the thermodynamic equilibration during the process of binding with a protein. The calculated phenomenological pKa of 3.66 is in good agreement with the experimental pKa of 3.9 reported in literature, suggesting that the computational predictions resulted from this study are reasonable.

Theoretical analysis of solvent effects on nitrogen NMR chemical shifts in oxazoles and oxadiazoles

Using quantum chemistry methods we have evaluated the solvent effects on the (14)N NMR chemical shifts in five oxa- and oxadiazoles dissolved in twelve solvents. These solvents differ in their polarity with the dielectric constants varying from 2 to 80. Moreover, three of them have a hydrogen-bond donor character. All possible hydrogen-bonding in the water solution with the oxygen and nitrogen (hydrogen-acceptor) centers in oxazoles (2) and oxadiazoles (3) have been considered in our studies. It has been shown that both the pure solvent and hydrogen-bonding effects are significant and result in (14)N magnetic shielding increase. In water solutions the pure solvent effect is larger than the hydrogen-bonding effect. In addition, the solvent effect has been analyzed in terms of its direct and indirect contributions. It should be emphasized that our theoretical results for (14)N chemical shifts in oxa- and oxadiazoles remain in a very good agreement with the accurate experimental data.

An efficient implementation for determining volume polarization in self-consistent reaction field theory

An efficient algorithm of the surface and volume polarization for electrostatics (SVPE) method in self-consistent reaction field (SCRF) theory, denoted by SV(1)PE, has been proposed to simulate direct volume polarization potential with a single layer of point charges outside the solute cavity while the indirect effects of volume polarization on surface polarization are still simulated with multiple layers of point charges. The free energies of solvation calculated using the SV(1)PE algorithm (implemented in GAUSSIAN03) reproduce the corresponding values calculated using the standard SVPE implementation within an error of only approximately 0.1% when the solute cavity is defined by the standard 0.001e/a(0) (3) solute charge isodensity contour. The SV(1)PE results are much less sensitive to the used cavity size in comparison with the well-established surface and simulated volume polarization for electrostatics [SS(V)PE] method which simulates volume polarization through an additional surface charge distribution on the cavity surface. The SCRF calculations using the SV(1)PE method are more efficient than those using the original SVPE method.

Modeling the catalysis of anti-cocaine catalytic antibody: competing reaction pathways and free energy barriers

The competing reaction pathways and the corresponding free energy barriers for cocaine hydrolysis catalyzed by an anti-cocaine catalytic antibody, mAb15A10, were studied by using a novel computational strategy based on the binding free energy calculations on the antibody binding with cocaine and transition states. The calculated binding free energies were used to evaluate the free energy barrier shift from the cocaine hydrolysis in water to the antibody-catalyzed cocaine hydrolysis for each reaction pathway. The free energy barriers for the antibody-catalyzed cocaine hydrolysis were predicted to be the corresponding free energy barriers for the cocaine hydrolysis in water plus the calculated free energy barrier shifts. The calculated free energy barrier shift of -6.87 kcal/mol from the dominant reaction pathway of the cocaine benzoyl ester hydrolysis in water to the dominant reaction pathway of the antibody-catalyzed cocaine hydrolysis is in good agreement with the experimentally derived free energy barrier shift of -5.93 kcal/mol. The calculated mutation-caused shifts of the free energy barrier are also reasonably close to the available experimental activity data. The good agreement suggests that the protocol for calculating the free energy barrier shift from the cocaine hydrolysis in water to the antibody-catalyzed cocaine hydrolysis may be used in future rational design of possible high-activity mutants of the antibody as anti-cocaine therapeutics. The general strategy of the free energy barrier shift calculation may also be valuable in studying a variety of chemical reactions catalyzed by other antibodies or proteins through noncovalent bonding interactions with the substrates.

First-principles calculation of pKa for cocaine, nicotine, neurotransmitters, and anilines in aqueous solution

The absolute pKa values of 24 representative amine compounds, including cocaine, nicotine, 10 neurotransmitters, and 12 anilines, in aqueous solution were calculated by performing first-principles electronic structure calculations that account for the solvent effects using four different solvation models, i.e., the surface and volume polarization for electrostatic interaction (SVPE) model, the standard polarizable continuum model (PCM), the integral equation formalism for the polarizable continuum model (IEFPCM), and the conductor-like screening solvation model (COSMO). Within the examined computational methods, the calculations using the SVPE model lead to the absolute pKa values with the smallest root-mean-square-deviation (rmsd) value (1.18). When the SVPE model was replaced by the PCM, IEFPCM, and COSMO, the rmsd value of the calculated absolute pKa values became 3.21, 2.72, and 3.08, respectively. All types of calculated pKa values linearly correlate with the experimental pKa values very well. With the empirical corrections using the linear correlation relationships, the theoretical pKa values are much closer to the corresponding experimental data and the rmsd values become 0.51-0.83. The smallest rmsd value (0.51) is also associated with the SVPE model. All of the results suggest that the first-principles electronic structure calculations using the SVPE model are a reliable approach to the pKa prediction for the amine compounds.

Vibrational anharmonic calculations in solution: Performance of various DFT approaches

We report anharmonic spectra calculated for formaldehyde in acetonitrile solution using the quartic force field obtained for various DFT/solvent coupled models. A statistical study has been carried out for each mode by using several classes of DFT functionals and comparing them to the reference ab-initio CCSD(T)/cc-pVQZ calculations. Results lead to the recommended use of hybrid functionals associated with the 6-31+G** basis set and the Polarized Continuum model (PCM) to predict the expected shifts relative to the gas phase.

Theoretical studies of the transition-state structures and free energy barriers for base-catalyzed hydrolysis of amides

The transition-state structures and free energy barriers for the rate-determining step (i.e. the formation of a tetrahedral intermediate) of base-catalyzed hydrolysis of a series of amides in aqueous solution have been studied by performing first-principle electronic structure calculations using a hybrid supermolecule-polarizable continuum approach. The calculated results and a revisit of recently reported experimental proton inventory data reveal that the favorable transition-state structure optimized for the tetrahedral intermediate formation of hydroxide ion-catalyzed hydrolysis of formamide may have three solvating water molecules remaining on the attacking hydroxide oxygen and two additional water molecules attached to the carbonyl oxygen of formamide. The calculated results have also demonstrated interesting substituent effects on the optimized transition-state geometries, on the transition-state stabilization, and on the calculated free energy barriers for the base-catalyzed hydrolysis of amides. When some or all of the hydrogen atoms of formamide are replaced by methyl groups, the total number of water molecules hydrogen-bonding with the attacking hydroxide in the transition state decreases from three for formamide to two for N-methylacetamide, N,N-dimethylformamide (DMF), and N,N-dimethylacetamide (DMA). The larger substituents of the amide hinder the solvent water molecules approaching the attacking hydroxide oxygen in the transition state and, therefore, destabilize the transition-state structure and increase the free energy barrier. By using the optimized most favorable transition-state structures, the calculated free energy barriers, i.e., 21.6 (or 21.7), 22.7, 23.1, and 26.0 kcal/mol for formamide, N-methylacetamide, DMF, and DMA, respectively, are in good agreement with the available experimental free energy barriers, i.e., 21.2, 21.5, 22.6, and 24.1 kcal/mol for formamide, N-methylacetamide, DMF, and DMA, respectively.

Nuclear magnetic shielding constants of liquid water: Insights from hybrid quantum mechanics/molecular mechanics models

We present a gauge-origin independent method for the calculation of nuclear magnetic shielding tensors of molecules in a structured and polarizable environment. The method is based on a combination of density functional theory (DFT) or Hartree-Fock wave functions with molecular mechanics. The method is unique in the sense that it includes three important properties that need to be fulfilled in accurate calculations of nuclear magnetic shielding constants: (i) the model includes electron correlation effects, (ii) the model uses gauge-including atomic orbitals to give gauge-origin independent results, and (iii) the effect of the environment is treated self-consistently using a discrete reaction-field methodology. The authors present sample calculations of the isotropic nuclear magnetic shielding constants of liquid water based on a large number of solute-solvent configurations derived from molecular dynamics simulations employing potentials which treat solvent polarization either explicitly or implicitly. For both the (17)O and (1)H isotropic shielding constants the best predicted results compare fairly well with the experimental data, i.e., they reproduce the experimental solvent shifts to within 4 ppm for the (17)O shielding and 1 ppm for the (1)H shielding.

New formulation and implementation for volume polarization in dielectric continuum theory

In the use of dielectric continuum theory to model bulk solvation effects on the electronic structure and properties of a solute, volume polarization contributions due to quantum mechanical penetration of the solute charge density outside the cavity nominally enclosing it are known to be significant. This work provides a new formulation and implementation of methods for solution of the requisite Poisson equation. In previous formulations the determination of the surface polarization contribution required evaluation of the difficult to calculate electric field generated by the volume polarization. It is shown that this problematic quantity can be eliminated in favor of other more easily evaluated quantities. That formal advance also opens the way for a more efficient apparatus to be implemented for calculation of the direct contribution of volume polarization to the solvation energy. The new formulation and its practical implementation are described, and illustrative numerical results are given for several neutral and ionic solutes to study the convergence and precision in practice.

A polarizable continuum approach for the study of heterogeneous dielectric environments

We present a computational method, exploiting some features of the polarizable continuum model (PCM) to describe heterogeneous media; it belongs to the family of electrostatic embedding mixed methods, such as the more common quantum-mechanical (QM)/molecular mechanics approaches, with the electrostatic long range effects accounted for by a polarized continuum instead of atomic point charges. Provided effective dielectric constants are determined for the various parts of the system, the method is much faster than its atomistic counterpart, and allows for high-level QM calculations on the fragment of interest, using all the highly efficient computational tools developed for homogeneous PCM. Two case studies (the calculation of the pKa of solvent exposed acidic residues in a model protein, and the calculation of the electron spin resonance spectrum of a typical spin probe partially embedded in a membrane) are analyzed in some detail, to illustrate the application of the method to complex systems.

Cation electric field is related to hydration energy

Electronic structure calculations on ions that use dielectric continuum theory to mimic solvent around the bare ionic solutes are often prone to make large errors in the hydration energies. It is found for cations that much of the error can be accounted for by a simple linear correlation with the maximum value found anywhere on the dielectric cavity surface of the solute potential or, even better, the outgoing normal electric field, thus mirroring analogous results previously obtained for anions. This correlation allows for significantly improved estimates of cation hydration energies while still retaining cavities of physically reasonable size in determination of the bulk dielectric contributions.

Comparison of frozen-density embedding and discrete reaction field solvent models for molecular properties

We investigate the performance of two discrete solvent models in connection with density functional theory (DFT) for the calculation of molecular properties. In our comparison we include the discrete reaction field (DRF) model, a combined quantum mechanics and molecular mechanics (QM/MM) model using a polarizable force field, and the frozen-density embedding (FDE) scheme. We employ these solvent models for ground state properties (dipole and quadrupole moments) and response properties (electronic excitation energies and frequency-dependent polarizabilities) of a water molecule in the liquid phase. It is found that both solvent models agree for ground state properties, while there are significant differences in the description of response properties. The origin of these differences is analyzed in detail and it is found that they are mainly caused by a different description of the ground state molecular orbitals of the solute. In addition, for the calculation of the polarizabilities, the inclusion of the response of the solvent to the polarization of the solute becomes important. This effect is included in the DRF model, but is missing in the FDE scheme. A way of including it in FDE calculations of the polarizabilities using finite field calculations is demonstrated.

Modeling multiple species of nicotine and deschloroepibatidine interacting with alpha4beta2 nicotinic acetylcholine receptor: from microscopic binding to phenomenological binding affinity

A variety of molecular modeling, molecular docking, and first-principles electronic structure calculations were performed to study how the alpha4beta2 nicotinic acetylcholine receptor (nAChR) binds with different species of two typical agonists, (S)-(-)-nicotine and (R)-(-)-deschloroepibatidine, each of which is distinguished by different free bases and protonation states. On the basis of these results, predictions were made regarding the corresponding microscopic binding free energies. Hydrogen-bonding and cation-pi interactions between the receptor and the respective ligands were found to be the dominant factors differentiating the binding strengths of different microscopic binding species. The calculated results and analyses demonstrate that, for each agonist, all the species are interchangeable and can quickly achieve a thermodynamic equilibrium in solution and at the nAChR binding site. This allows quantitation of the equilibrium concentration distributions of the free ligand species and the corresponding microscopic ligand-receptor binding species, their pH dependence, and their contributions to the phenomenological binding affinity. The predicted equilibrium concentration distributions, pK(a) values, absolute phenomenological binding affinities, and their pH dependence are all in good agreement with available experimental data, suggesting that the computational strategy from the microscopic binding species and affinities to the phenomenological binding affinity is reliable for studying alpha4beta2 nAChR-ligand binding. This should provide valuable information for future rational design of drugs targeting nAChRs. The general strategy of the "from-microscopic-to-phenomenological" approach for studying interactions of alpha4beta2 nAChRs with (S)-(-)-nicotine and (R)-(-)-deschloroepibatidine may also be useful in studying other types of ligand-protein interactions involving multiple molecular species of a ligand and in associated rational drug design.

First-principle studies of intermolecular and intramolecular catalysis of protonated cocaine

We have performed a series of first-principles electronic structure calculations to examine the reaction pathways and the corresponding free energy barriers for the ester hydrolysis of protonated cocaine in its chair and boat conformations. The calculated free energy barriers for the benzoyl ester hydrolysis of protonated chair cocaine are close to the corresponding barriers calculated for the benzoyl ester hydrolysis of neutral cocaine. However, the free energy barrier calculated for the methyl ester hydrolysis of protonated cocaine in its chair conformation is significantly lower than for the methyl ester hydrolysis of neutral cocaine and for the dominant pathway of the benzoyl ester hydrolysis of protonated cocaine. The significant decrease of the free energy barrier, approximately 4 kcal/mol, is attributed to the intramolecular acid catalysis of the methyl ester hydrolysis of protonated cocaine, because the transition state structure is stabilized by the strong hydrogen bond between the carbonyl oxygen of the methyl ester moiety and the protonated tropane N. The relative magnitudes of the free energy barriers calculated for different pathways of the ester hydrolysis of protonated chair cocaine are consistent with the experimental kinetic data for cocaine hydrolysis under physiologic conditions. Similar intramolecular acid catalysis also occurs for the benzoyl ester hydrolysis of (protonated) boat cocaine in the physiologic condition, although the contribution of the intramolecular hydrogen bonding to transition state stabilization is negligible. Nonetheless, the predictability of the intramolecular hydrogen bonding could be useful in generating antibody-based catalysts that recruit cocaine to the boat conformation and an analog that elicited antibodies to approximate the protonated tropane N and the benzoyl O more closely than the natural boat conformer might increase the contribution from hydrogen bonding. Such a stable analog of the transition state for intramolecular catalysis of cocaine benzoyl-ester hydrolysis was synthesized and used to successfully elicit a number of anticocaine catalytic antibodies.

Reaction Pathways and Free Energy Barriers for Alkaline Hydrolysis of Insecticide 2-Trimethylammonioethyl Methylphosphonofluoridate and Related Organophosphorus Compounds: Electrostatic and Steric Effects

Reaction pathways and free energy barriers for alkaline hydrolysis of the highly neurotoxic insecticide 2-trimethylammonioethyl methylphosphonofluoridate and related organophosphorus compounds were studied by performing first-principles electronic structure calculations on representative methylphosphonofluoridates, (RO)CH3P(O)F, in which R = CH2CH2N+(CH3)3, CH3, CH2CH2C(CH3)3, CH2CH2CH(CH3)2, CH(CH3)CH2N+(CH3)3, and CH(CH3)CH2N(CH3)2. The dominant reaction pathway was found to be associated with a transition state in which the attacking nucleophile OH- and the leaving group F- are positioned on opposite sides of the plane formed by the three remaining atoms attached to the phosphorus in order to minimize the electrostatic repulsion between these two groups. The free energy barriers calculated for the rate-determining step of the dominant pathway are 12.5 kcal/mol when R = CH2CH2N+(CH3)3, 15.5 kcal/mol when R = CH3, 17.9 kcal/mol when R = CH2CH2C(CH3)3, 16.5 kcal/mol when R = CH2CH2CH(CH3)2, 13.4 kcal/mol when R = CH(CH3)CH2N+(CH3)3, and 18.7 kcal/mol when R = CH(CH(3))CH(2)N(CH(3))(2). The calculated free energy barriers are in good agreement with available experimentally derived activation free energies, i.e. 14.7 kcal/mol when R = CH(3), 13.4 kcal/mol when R = CH2CH2N+(CH3)3, and 13.9 kcal/mol when R = CH(CH3)CH2N+(CH3)3. A detailed analysis of the calculated energetic results and available experimental data suggests that the net charge of the molecule (M) being hydrolyzed is a prominent factor affecting the free energy barrier (DeltaG) for the alkaline hydrolysis of phosphodiesters, phosphonofluoridates, and related organophosphorus compounds. The electrostatic interactions between the attacking nucleophile OH- and the molecule M being hydrolyzed favor such an order of the free energy barrier: DeltaG(M(+)+OH-) < DeltaG(M0+OH-) < DeltaG(M(-)+OH-), where M+, M0, and M- represent the cationic, neutral, and anionic molecules, respectively. The change of the substituent R in (RO)CH(3)P(O)F from CH3 to CH2CH2N+(CH3)3 is associated with both the electrostatic and steric effects on the free energy barrier, but the electrostatic effect dominates the substituent shift of the free energy barrier. This helps to better understand why the alkaline hydrolysis of (RO)CH3P(O)F with R = CH2CH2N+(CH3)3 and CH(CH3)CH2N+(CH3)3 is significantly faster than that with R = CH3. The effect of electrostatic interaction also helps to understand why the rate constants for the alkaline hydrolysis of phosphodiesters, such as intramolecular second messenger adenosine 3',5'-phosphate (cAMP), are generally smaller than those for the alkaline hydrolysis of the phosphonofluoridates and related phosphotriesters.

The determinants of carboxyl pKa values in turkey ovomucoid third domain

A computational methodology for protein pK(a) predictions, based on ab initio quantum mechanical treatment of part of the protein and linear Poisson-Boltzmann equation treatment of the bulk solvent, is presented. The method is used to predict and interpret the pK(a) values of the five carboxyl residues (Asp7, Glu10, Glu19, Asp27, and Glu43) in the serine protease inhibitor turkey ovomucoid third domain. All the predicted pK(a) values are within 0.5 pH units of experiment, with a root-mean-square deviation of 0.31 pH units. We show that the decreased pK(a) values observed for some of the residues are primarily due to hydrogen bonds to the carboxyl oxygens. Hydrogen bonds involving amide protons are shown to be particularly important, and the effect of hydrogen bonding is shown to be nonadditive. Hydrophobic effects are also shown to be important in raising the pK(a). Interactions with charged residues are shown to have relatively little effect on the carboxyl pK(a) values in this protein, in general agreement with experiment.

Fundamental reaction mechanism for cocaine hydrolysis in human butyrylcholinesterase

Butyrylcholinesterase (BChE)-cocaine binding and the fundamental pathway for BChE-catalyzed hydrolysis of cocaine have been studied by molecular modeling, molecular dynamics (MD) simulations, and ab initio calculations. Modeling and simulations indicate that the structures of the prereactive BChE/substrate complexes for (-)-cocaine and (+)-cocaine are all similar to that of the corresponding prereactive BChE/butyrylcholine (BCh) complex. The overall binding of BChE with (-)-cocaine and (+)-cocaine is also similar to that proposed with butyrylthiocholine and succinyldithiocholine, i.e., (-)- or (+)-cocaine first slides down the substrate-binding gorge to bind to Trp-82 and stands vertically in the gorge between Asp-70 and Trp-82 (nonprereactive complex) and then rotates to a position in the catalytic site within a favorable distance for nucleophilic attack and hydrolysis by Ser-198 (prereactive complex). In the prereactive complex, cocaine lies horizontally at the bottom of the gorge. The fundamental catalytic hydrolysis pathway, consisting of acylation and deacylation stages similar to those for ester hydrolysis by other serine hydrolases, was proposed on the basis of the simulated prereactive complex and confirmed theoretically by ab initio reaction coordinate calculations. Both the acylation and deacylation follow a double-proton-transfer mechanism. The calculated energetic results show that within the chemical reaction process the highest energy barrier and Gibbs free energy barrier are all associated with the first step of deacylation. The calculated ratio of the rate constant (k(cat)) for the catalytic hydrolysis to that (k(0)) for the spontaneous hydrolysis is approximately 9.0 x 10(7). The estimated k(cat)/k(0) value of approximately 9.0 x 10(7) is in excellent agreement with the experimentally derived k(cat)/k(0) value of approximately 7.2 x 10(7) for (+)-cocaine, whereas it is approximately 2000 times larger than the experimentally derived k(cat)/k(0) value of approximately 4.4 x 10(4) for (-)-cocaine. All of the results suggest that the rate-determining step of the BChE-catalyzed hydrolysis of (+)-cocaine is the first step of deacylation, whereas for (-)-cocaine the change from the nonprereactive complex to the prereactive complex is rate-determining and has a Gibbs free energy barrier higher than that for the first step of deacylation by approximately 4 kcal/mol. A further analysis of the structural changes from the nonprereactive complex to the prereactive complex reveals specific amino acid residues hindering the structural changes, providing initial clues for the rational design of BChE mutants with improved catalytic activity for (-)-cocaine.