A Comparison of Seven Fast but Approximate Methods to Compute the Free Energy of Deprotonation for Amino Acids in Aqueous Solution

Abstract

In recent years, a variety of methods based on statistical mechanics have been successfully applied to calculate free energy differences of chemical reactions from molecular simulation. The accuracy and computational efficiency vary strongly between these methods. Seven approximate but fast methods to calculate free energy differences are compared in terms of accuracy and efficiency with the accurate but expensive thermodynamic integration method as reference, using 28 protonation and deprotonation reactions of aspartic acid in aqueous solution as test cases. At least two simulations are required to obtain an accurate free energy difference between two states of the system. Both, the averaged one-step perturbation method and the linear response method yield the most accurate results, while the latter method shows the fastest convergence.     

Microsecond simulations of the folding/unfolding thermodynamics of the Trp‐cage miniprotein

We study the unbiased folding/unfolding thermodynamics of the Trp‐cage miniprotein using detailed molecular dynamics simulations of an all‐atom model of the protein in explicit solvent using the Amberff99SB force field. Replica‐exchange molecular dynamics simulations are used to sample the protein ensembles over a broad range of temperatures covering the folded and unfolded states at two densities. The obtained ensembles are shown to reach equilibrium in the 1 μs/replica timescale. The total simulation time used in the calculations exceeds 100 μs. Ensemble averages of the fraction folded, pressure, and energy differences between the folded and unfolded states as a function of temperature are used to model the free energy of the folding transition, ΔG(P, T), over the whole region of temperatures and pressures sampled in the simulations. The ΔG(P, T) diagram describes an ellipse over the range of temperatures and pressures sampled, predicting that the system can undergo pressure‐induced unfolding and cold denaturation at low temperatures and high pressures, and unfolding at low pressures and high temperatures. The calculated free energy function exhibits remarkably good agreement with the experimental folding transition temperature (T_f = 321 K), free energy, and specific heat changes. However, changes in enthalpy and entropy are significantly different than the experimental values. We speculate that these differences may be due to the simplicity of the semiempirical force field used in the simulations and that more elaborate force fields may be required to describe appropriately the thermodynamics of proteins. Proteins 2010. © 2010 Wiley‐Liss, Inc.     

次氯酸与不饱和脂肪酸氧化机制和转化产物的理论研究

摘要 目的：揭示次氯酸与不饱和脂肪酸的氧化反应机制及转化产物。方法：运用Gaussian 16软件包,采用密度泛函方法M06-2X(D3),结合6-31+G(d)基组,在SMD液相水模型水平下进行计算。结果：次氯酸与单不饱和脂肪酸油酸的氧化反应是先形成氯鎓离子中间体,氯鎓离子再与水分子反应生成氯醇,第一步氯鎓离子的形成是控速步骤,其反应活化自由能~8 kcal/mol。环氧化合物和短链的醛是两种转化产物,前者由氯醇脱氯化氢而来,而后者由环氧化合物和氯醇通过系列与次氯酸根的反应而得到,生成它们的控速步骤的反应活化自由能分别为23 和24 kcal/mol。选取两个乙基为取代基的乙烯为油酸模型,其与次氯酸反应的活化自由能仅比油酸高1 kcal/mol。计算得到次氯酸与亚油酸、顺-9,反-11 亚油酸、梓树酸和花生四烯酸模型氧化反应生成氯醇的活化自由能分别是~10、13、16和14 kcal/mol。结论：氯鎓离子中间体机制是次氯酸与不饱和脂肪酸氧化反应的主要机制,反应的活化自由能通常低于15 kcal/mol,意味着此氧化反应动力学上容易发生。氧化产物氯醇能转化为环氧化合物和短链的醛,但活化自由能较高,约23和24 kcal/mol。选取距离双键3个碳以内的结构为不饱和脂肪酸模型,它能够很好地反映不饱和脂肪酸的反应活性。     

Microcalorimetric study of glutamine fixation on the glutamine-binding protein of Escherichia coli.

The enthalpy variation (ΔH) induced by addition of glutamine to glutamine binding protein isolated from $\text{E. coli}$ has been studied by microcalorimetry. The reaction was very exothermic. The free energy variation (ΔG) was calculated from the dissociation constant (K_D) measured by dialysis techniques. The entropic variation (ΔS) was deduced from ΔG and ΔH values; it was found highly negative, indicating that an important conformational change is occuring. Comparison with others binding proteins and possible significance of such a phenomenon is discussed.     

Enzymatic synthesis of β-lactam antibiotics: A thermodynamic background

The equilibrium constants and the respective standard Gibbs energy changes for hydrolysis of some β-lactam antibiotics have been determined. Native and immobilized penicillin amidase (EC 3.5.1.11) from Escherichia coli has been used as a catalyst. The values of standard Gibbs energy changes corresponding to the pH-independent product of equilibrium concentrations (ΔG⁰_c = ? RT ln K_c) have been calculated. The differences in the structure of the antibiotics nucleus hardly ever affect the value of the pH-independent component of the standard Gibbs energy change (ΔG⁰_c) and value of apparent standard Gibbs energy change at a fixed pH (ΔG^0′_c). At the same time, the value of ΔG⁰_c is more sensitive to the structure of the acyl moiety of the antibiotic; when ampicillin is used instead of benzylpenicillin, ΔG⁰_c increases by ～6.3 kJ mol^?1 (1.5 kcal mol^?1). pH-dependences of the apparent standard Gibbs energy changes for hydrolysis of β-lactam antibiotics have been calculated. The pH-dependences of ΔG^0′_c for hydrolysis of all β-lactam antibiotics have a similar pattern. The thermodynamic pH optimum of the synthesis of these compounds is in the acid pH range (pH < 5.0). The breakage of the β-lactam ring leads to a sharp decrease in the ΔG^0′_c value and a change in the pattern of the pH-dependence. For example, at pH 5.0 ΔG^0′_c decreases from 14.4 kJ mol^?1 for benzylpenicillin to ?1.45 kJ mol^?1 for benzylpenicilloic acid. The reason for these changes is mainly a considerable increase in the pK of the amino group of the nucleus of the antibiotic and, as a consequence, a decrease in the component of standard Gibbs energy change, corresponding to the ionization of the system. The thermodynamic potentials of the enzymatic synthesis of semisynthetic penicillins and cephalosporins on the basis of both free acids and their derivatives (N-acylated amino acids, esters) are discussed. It is shown that with esters of the acids, a high yield of the antibiotic can, in principle, be achieved at higher pH values.     

Linear free energy relationship for the anomeric effect: MP2, DFT and ab initio study of 2-substituted-1,4-dioxanes

The anomeric effect of 2-substituted 1,4-dioxane derivatives was calculated and compared with the values for substituted cyclohexane. The bond lengths, bond angles, torsion angles, and relative energies of axial and equatorial conformers of 2-substituted 1,4-dioxanes were calculated by the second-order Møller–Plesset (MP2), density functional theory (DFT/B3LYP), and Hartree–Fock (HF) methods using 6-31G^∗ basis set. The energy differences between the axial and equatorial conformers, endo and exo-anomeric effects, repulsive non-bond and H-bonding interactions were investigated. A linear free energy relationship (LFER) between calculated (MP2/6-31G^∗) anomeric effect and inductive substituent constants (σ_I) was obtained for 2-substituted-1,4-dioxanes (slope = 6.19 and r² = 0.967). The calculated energy differences indicate lower equatorial orientation for 2-substituted-1,4-dioxanes compared to the 2-substituted-tetrahydropyrans. The contribution of resonance, hyperconjugation, inductive, steric, hydrogen bonding, electrostatic interaction, and level of theory influences the anomeric effect.     

Effective prediction model and determination of binding residues influential for inhibitors targeting HIV-1 integrase-LEDGF/p75 interface by employing solvent accessible surface area energy as key determinant

Abstract

Development of a highly accurate prediction model for protein–ligand inhibition has been a major challenge in drug discovery. Herein, we describe a novel predictive model for the inhibition of HIV-1 integrase (IN)-LEDGF/p75 protein-protein interaction. The model was constructed using energy parameters approximated from molecular dynamics (MD) simulations and molecular mechanics Poisson-Boltzmann surface area (MM-PBSA) calculations. Chemometric analysis using partial least squares (PLS) regression revealed that solvent accessible surface area energy (ΔG_SASA) is the major determinant parameter contributing greatly to the prediction accuracy. PLS prediction model on the ΔG_SASA values collected from 41 complexes yielded a strong correlation between the predicted and the actual inhibitory activities (R² = 0.9666, RMSEC of pIC₅₀ values = 0.0890). Additionally, for the test set of 14 complexes, the model performed satisfactorily with very low pIC₅₀ errors (Q² = 0.5168, RMSEP = 0.3325). A strong correlation between the buried surface areas on the IN protein, when bound with IN-LEDGF/p75 inhibitors, and the respective ΔG_SASA values was also obtained. Furthermore, the current method could identify ‘hot spots’of amino acid residues highly influential to the inhibitory activity prediction. This could present fruitful implications in binding site determination and future inhibitor developments targeting protein-protein interactions.

Communicated by Ramaswamy H. Sarma     

Complete basis set,hybrid-DFT study and NBO interpretation of conformational analysis of 2-methoxytetrahydropyran and its thiopyran and selenopyran analogues in relation to the anomeric effect

2-Methoxytetrahydropyran (1), -thiopyran (2) and -selenopyran (3) have been chosen as model compounds to investigate the origin of the anomeric effect (AE). The impacts of the hyperconjugation, electrostatic and steric interactions on the conformational preferences of compounds 1–3 have been analysed by means of complete basis set-4, hybrid-density functional theory (B3LYP/6-311+G**) based methods and natural bond orbital (NBO) interpretation. Both levels of theory showed that the axial conformations of compounds 1–3 are more stable than their equatorial conformations. The Gibbs free energy difference (G _eq–G _ax) values (i.e. ΔG _eq–ax) between the axial and equatorial conformations increase from compound 1 to compound 2 but decrease from compound 2 to compound 3. Based on the NBO results obtained, the AE associated with the electron delocalisation [i.e. Σ(endo-AE_eq + exo-AE_eq) ? Σ(endo-AE_ax + exo-AE_ax)] increase slightly from compound 1 to compound 2 but decrease from compound 2 to compound 3. Similar trend is also observed for the differences between the calculated total steric exchange energy values [i.e. Δ(TSEE)_eq–ax]. On the other hand, the calculated differences between the dipole moment values of the axial and equatorial conformations [i.e. Δ(μ_eq–μ_ax)] decrease from compound 1 to compound 3. These findings led to the proposal that the AE associated with the electron delocalisation (the hyperconjugation effect) is more significant for the explanation of the conformational preferences of compounds 1–3 than the electrostatic model. The correlations between the AE associated with the electron delocalisation, bond orders, TSEE, ΔG _eq–ax, dipole–dipole interactions, structural parameters and conformational behaviours of compounds 1–3 have been investigated.     

8 Using pseudorotation as a reaction coordinate in free energy simulations of nucleic acids

Backbone sugar groups are central components of nucleic acids. The conformations of the ribose/deoxyribose can be elegantly described using the concept of pseudorotation (Altona and Sundaralingam, 1972), and are dominated by the C2′- and C3′-endo conformers. The free energy barrier of the transition between these two major puckering modes can be probed by NMR relaxation experiments (Johnson and Hoogstraten, 2008), but an atomic picture of the transition path per se is only available for several truncated nucleoside analogues (Brameld & Goddard III, 1999). Here, we implemented a new free energy simulation method for Molecular Dynamics simulations using pseudorotation as the reaction coordinate (Cremer and Pople, 1975). This allowed us to compute the free energy landscape of a complete pseudorotation cycle. The free energy landscape revealed not only the relative stability of C2′- and C3′-endo conformers, but also the main transition path and its free energy barrier. As a validation of our new approach, we calculated free energy surface of the pseudorotation of guanosine monophosphate. The free energy surface revealed that the C2′-endo conformation is ?1?kcal/mol that is more stable and the free energy barrier for the transition is 4.5–5?kcal/mol. These are in excellent agreement with previous NMR measurements (Zhang et al., 2012; Röder et al., 1975). We have further applied this method to other systems that are important in pre-biotic chemistry, including an RNA duplex with unique 2′, 5′-phosphodiester linkages.     

Thermodynamics of the first transition in writhe of a small circular DNA by Monte Carlo simulation

Monte Carlo simulations are employed to investigate the thermodynamics of the first transition in writhe of a circular model filament corresponding to a 468 base-pair DNA. Parameters employed in these simulations are the torsional rigidity, C = 2.0 × 10⁻¹⁹ dyne cm², and persistence length, P = 500 Å. Intersubunit interactions are modeled by a screened Coulomb potential. For a straight line of subunits this accurately approximates the nonlinear Poisson-Boltzmann potential of a cylinder with the linear charge density of DNA. Curves of relative free energy vs writhe at fixed linking difference (Δ1) exhibit two minima, one corresponding to slightly writhed circles and one to slightly underwrithed figure-8's, whenever Δ1 lies in the transition region. The free energies of the two minima are equal when Δ1_c = 1.35, which defines the midpoint of the transition. At this midpoint, the free energy barrier between the two minima is found to be ΔG_bar = (0.20) k_BT at 298 K. Curves of mean potential energy vs writhe at fixed linking difference similarly exhibit two minima for Δ1 values in the transition region, and the two minimum mean potential energies are equal when Δ1 = 1.50. At the midpoint writhe, Δ1_c = 1.35, the difference in mean potential energy between the minimum free energy figure-8 and circle states is (1.3) k_BT, and the difference in their entropies is 1.3 k_B. Thus, the entropy of the minimum free energy figure-8 state significantly exceeds that of the circle at the midpoint of the transition. The first transition in writhe is found to occur over a rather broad range of Δ1 values from 0.85 to 1.85. The twist energy parameter (E_T), which governs the overall free energy of supercoiling, undergoes a sigmoidal decrease, while the translational diffusion coefficient undergoes a sigmoidal increase, over this same range. The static structure factor exhibits an increase, which reflects a decrease in radius of gyration associated with the circle to figure-8 transition. © 1996 John Wiley & Sons, Inc.     

Theoretical vibrational spectra of cyclohexanecarboxaldehyde

Vibrational spectra of cyclohexanecarboxaldehyde are calculated with density functional theory using the B3LYP functional together with a 6-311++G** basis set and presented. The results in the case of the infrared spectrum of the mixture of conformers at 300 K are compared with the experimental spectrum and, apart from the intensities of the CH and CO stretches, reasonable agreement is found. These deficiencies can be traced back to the well-known nonlinearities in case of CH stretches and the CO stretch. Potential energy distributions among symmetry coordinates in each normal mode are presented and used to assign specific atomic movements to each of the modes. Potential energy scans for the CHO rotor in both the axial and equatorial conformers are presented and barrier heights are compared with previous Hartree–Fock calculations and experimental data. It is reported that there could be three stable conformers, namely equatorial-gauche (eg), which is the most stable, equatorial-trans (et) and axial-gauche (ag). The optimized energies of all the minima and of the transition states are presented. However, comparison of the calculated spectrum with the experimental one indicates that total energies are slightly in error and that in the mixture of conformers no ag is present and thus the et to eg ratio is also different. Using experimental values for relative energies of conformers, we could obtain spectra in fair agreement with experiment. This indicates that when only total energy differences are calculated, slight errors in them play a role because of the very small relative energies in this case, while properties like geometries and spectra, which depend not on energy differences but on analytically calculated energy derivatives, are not affected.Electronic Supplementary Material available.     

Dependence of Deformability of Geometries and Characteristics of Intramolecular Hydrogen Bonds in Canonical 2′-Deoxyribonucleotides on DNA Conformations

Abstract

The molecular structure and deformability (with respect to average geometry) of methyl ethers of canonical 2′-deoxyribonucleotides thymidine-5′-phosphate (mTMP), 2-deoxycytidine-5′- phosphate (mCMP), 2-deoxyadenosine-5′-phosphate (mAMP) and 2′-deoxyguanosine-5′- phosphate (mGMP) in different types of DNA have been calculated using B3LYP/cc-pvdz method. Comparison of energy at equilibrium conformations of nucleotides and conformations with torsion angles of backbone fixed to average values for different types of DNA reveals that incorporation of nucleotides to A-DNA macromolecules requires the minimum amount of deformation energy. Therefore, this type of DNA should be the least strained from viewpoint of intramolecular deformations of monomers. Modeling of environmental effects within the PCM approach reveals that the immersion of nucleotides in polar medium results in significant decrease of energy differences between anti conformers of all DNTs and syn conformers of mGMP This also leads to reduction by almost a half nucleotides' deformation energy facilitating formation of DNA macromolecule. Change of DNTs conformation causes switch between different types of intramolecular H bonds. Every type of DNA possesses unique set of intramolar hydrogen bonds in nucleotides.     

Computation of the bond dissociation enthalpies and free energies of hydroxylic antioxidants using the ab initio Hartree–Fock method

Abstract

Introduction

A new method for calculating theoretical bond dissociation enthalpy (BDE) and bond dissociation free energy (BDFE) of hydroxylic antioxidants is forwarded. BDE and BDFE may be understood as activation energies accompanying the formation of transition states, which may undergo downhill homolytic dissociation. The new method does not involve the complete fission of O–H bonds.

Method

Theoretical gas phase BDE values were calculated with the ab initio unrestricted Hartree–Fock (UHF) method, as changes in enthalpy between ground singlet states (GS) and triplet dissociative states (DS). Similarly, gas phase BDFEs were estimated from the corresponding changes in Gibbs free energy. The results were then compared with reliable experimental reports.

Results

The proposed theoretical approach of BDE and BDFE determination was tested using 10 simple phenols, 5 flavonoids, and l-ascorbic acid derivatives. The agreement between our calculated gas phase results and the adopted experimental values were generally within 0.5 kcal mol^?1, with a very few exceptions.

Discussion

Generally, steric interactions as well as intramolecular hydrogen bonding involving the dissociating OH group should be minimized in the GS. The DS are both electronically and vibrationally exited transition states. They have one unpaired electron on the carbon atom, which bears the homolytically dissociating OH group and are second order saddle points with a fixed <C–O–H bond angel of 180°.

Conclusion

It was concluded that ab initio UHF was well suited for the estimation of gas phase BDE and BDFE. The method presented has a good potential for application across a range of hydroxylic antioxidants. Currently, work is underway to extend its application in other class of antioxidants.     

Dynamics and Stability of GCAA Tetraloops with 2-Aminopurine and Purine Substitutions

Abstract

The contributions of various interactions in the GGCGCAAGCC hairpin containing a GCAA tetraloop were studied by computer simulations using the substitutions of functional groups. The guanosine (G) in the first tetraloop position or in the C-G closing base pair was replaced by 2-aminopurine (AP), and the individual tetraloop's adenosines (A) were replaced by purine (PUR). These substitutions eliminated particular hydrogen bonds thought to stabilize the GCAA tetraloop. For each substitution, molecular dynamics (MD) simulations were carried out in an aqueous solution with sodium counterions, using the CHARMM27 force field. The MD simulations showed that the substitutions in the first (G→AP) and the third (A→PUR) position of the GCAA tetraloop did not significantly influence the conformation of the hairpin. A long-lived bridging water molecule observed in the GCAA loop was present in both modified loops. The substitutions made in the last loop position (A→PUR) or in the C-G base pair closing the tetraloop (G→AP) to some extent influenced the loop structure and dynamics. These loops did not display the long- lived bridging water molecules. When the second A in the GCAA loop was replaced by PUR, the first A in the loop was observed in the anti or in the syn orientation about the gly- cosyl bond. The G to AP substitution in C-G base pair led to a change of their arrangement from the Watson-Crick to wobble. The MD simulations of the hairpin with C-AP wobble closing base pair showed increased conformational dynamics of the hairpin. The changes of hairpin formation free energy associated with the substitutions of individual bases were calculated by the free energy perturbation method. Our theoretical estimates suggest a larger destabilization for the G to AP substitutions in GCAA loop than for the substitutions of individual A's by PUR, which is in accordance with experimental tendency. The calculations predicted a similar free energy change for G to AP substitutions in the GCAA tetraloop and in the C-G closing base pair.     

De(side chain) model of epothilone: bioconformer interconversions DFT study

Using ab initio methods, we have studied conformations of the de(sidechain)de(dioxy)difluoroepothilone model to quantify the effect of stability change between the exo and endo conformers of the epoxy ring. The DFT minimization of the macrolactone ring reveals four low energy conformers, although MP2 predicted five stable structures. The model tested with DFT hybride functional (B3LYP/6–31+G(d,p)) exhibits the global minimum for one of the exo forms (C), experimentally observed in the solid state, but unexpectedly with the MP2 electron correlation method for the virtual endo form (W). Using the QST3 technique, several pathways were found for the conversion of the low energy conformers to the other low energy exo representatives, as well as within the endo analog subset. The potential energy relationships obtained for several exo forms suggest a high conformational mobility between three, experimentally observed, conformers. The high rotational barrier, however, excludes direct equilibrium with experimental EC-derived endo form S. The highest calculated transition state for the conversion of the most stable exo M interligand to the endo S form is approximately a 28 kcal/mol above the energy of the former. The two-step interconversion of the exo H conformer to the endo S requires at least 28 kcal/mol. Surprisingly, we found that the transition state energy of the H form to the virtual endo W has the acceptable value of about 9 kcal/mol and the next energy barrier for free interconversion of endo W to endo S is 13 kcal/mol. Figure DFT Model of Epothilone Interconversions Electronic supplementary material  The online version of this article (doi:) contains supplementary material, which is available to authorized users.     

Binding of 1-substituted carbazolyl-3,4-dihydro-β-carbolines with DNA: Molecular dynamics simulation and MM-GBSA analysis

Molecular Mechanics-Generalized Born-Solvent Accessibility free energy calculations were used to analyse DNA binding affinity of 1-substituted carbazolyl-3,4-dihydro-β-carboline molecules. In this study, DNA structure with sequence of d(CGATCG)2 was used for simulations. 15 ns molecular dynamics simulations of the studied complexes were performed. The calculated free energy was compared with experimental antitumor activity (IC₅₀). The predicted free energies decreased with the increase of IC₅₀ values. It was shown that molecules 1–6 bind to DNA via intercalation mode, while molecules 7–9 bind through groove binding mode. Also, it was found that the vdW energy term (ΔE_vdW) and the non-polar desolvation energy (ΔG_SA) are the favorable terms for binding energy, whereas net electrostatic energies (ΔE_ele + ΔG_GB) and conformational entropy energy (TΔS) are unfavorable ones.     

Conformational diversity of the quercetin molecule: a quantum-chemical view

Abstract

This paper focuses on the comprehensive conformational analysis of the quercetin molecule with a broad range of the therapeutic and biological actions. All possible conformers of these molecule, corresponding to the local minima on the potential energy hypersurface, have been obtained by the sequential rotation of its five hydroxyl groups and also by the rotation of its (A?+?C) and B rings relatively each other. Altogether, it was established 48 stable conformers, among which 24 conformers possess planar structure and 24 conformers – nonplanar structure. Their structural, symmetrical, energetical and polar characteristics have been investigated in details. Quantum-mechanical calculations indicate that conformers of the quercetin molecule are polar structures with a dipole moment, which varies within the range from 0.35 to 9.87 Debay for different conformers. Relative Gibbs free energies of these conformers are located within the range from 0.0 to 25.3?kcal·mol^?1 in vacuum under normal conditions. Impact of the continuum with ε?=?4 leads to the decreasing of the Gibbs free energies (–0.19–18.15?kcal·mol^?1) and increasing of the dipole moment (0.57–12.48?D). It was shown that conformers of the quercetin molecule differ from each other by the intramolecular specific contacts (two or three), stabilizing all possible conformers of the molecule – H-bonds (both classical ОН…О and so-called unusual С′Н…О and ОН…С′) and attractive van-der-Waals contacts О…О. Obtained conformational analysis for the quercetin molecule enables to provide deeper understanding of the ‘structure-function’ relationship and also to suggest its mechanisms of the therapeutic and biological actions.

Communicated by Ramaswamy H. Sarma