Computing Relative Free Energies of Solvation using Single Reference Thermodynamic Integration Augmented with Hamiltonian Replica Exchange

This paper introduces an efficient single-topology variant of Thermodynamic Integration (TI) for computing relative transformation free energies in a series of molecules with respect to a single reference state. The presented TI variant that we refer to as Single-Reference TI (SR-TI) combines well-established molecular simulation methodologies into a practical computational tool. Augmented with Hamiltonian Replica Exchange (HREX), the SR-TI variant can deliver enhanced sampling in select degrees of freedom. The utility of the SR-TI variant is demonstrated in calculations of relative solvation free energies for a series of benzene derivatives with increasing complexity. Noteworthy, the SR-TI variant with the HREX option provides converged results in a challenging case of an amide molecule with a high (13-15 kcal/mol) barrier for internal cis/trans interconversion using simulation times of only 1 to 4 ns.     

Computing Alchemical Free Energy Differences with Hamiltonian Replica Exchange Molecular Dynamics (H-REMD) Simulations

Alchemical free energy calculations play a very important role in the field of molecular modeling. Efforts have been made to improve the accuracy and precision of those calculations. One of the efforts is to employ a Hamiltonian replica exchange molecular dynamics (H-REMD) method to enhance conformational sampling. In this paper, we demonstrated that HREMD method not only improves convergence in alchemical free energy calculations but also can be used to compute free energy differences directly via the Free Energy Perturbation (FEP)algorithm. We show a direct mapping between the H-REMD and the usual FEP equations, which are then used directly to compute free energies. The H-REMD alchemical free energy calculation (Replica exchange Free Energy Perturbation, REFEP) was tested on predicting the pK(a) value of the buried Asp26 in thioredoxin. We compare the results of REFEP with TI and regular FEP simulations. REFEP calculations converged faster than those from TI and regular FEP simulations. The final predicted pK(a) value from the H-REMD simulation was also very accurate, only 0.4 pK(a) unit above the experimental value. Utilizing the REFEP algorithm significantly improves conformational sampling, and this in turn improves the convergence of alchemical free energy simulations.     

Free Energy Perturbation Hamiltonian Replica-Exchange Molecular Dynamics (FEP/H-REMD) for Absolute Ligand Binding Free Energy Calculations

Free Energy Perturbation with Replica Exchange Molecular Dynamics (FEP/REMD) offers a powerful strategy to improve the convergence of free energy computations. In particular, it has been shown previously that a FEP/REMD scheme allowing random moves within an extended replica ensemble of thermodynamic coupling parameters "lambda" can improve the statistical convergence in calculations of absolute binding free energy of ligands to proteins [J. Chem. Theory Comput. 2009, 5, 2583]. In the present study, FEP/REMD is extended and combined with an accelerated MD simulations method based on Hamiltonian replica-exchange MD (H-REMD) to overcome the additional problems arising from the existence of kinetically trapped conformations within the protein receptor. In the combined strategy, each system with a given thermodynamic coupling factor lambda in the extended ensemble is further coupled with a set of replicas evolving on a biased energy surface with boosting potentials used to accelerate the inter-conversion among different rotameric states of the side chains in the neighborhood of the binding site. Exchanges are allowed to occur alternatively along the axes corresponding to the thermodynamic coupling parameter lambda and the boosting potential, in an extended dual array of coupled lambda- and H-REMD simulations. The method is implemented on the basis of new extensions to the REPDSTR module of the biomolecular simulation program CHARMM. As an illustrative example, the absolute binding free energy of p-xylene to the nonpolar cavity of the L99A mutant of T4 lysozyme was calculated. The tests demonstrate that the dual lambda-REMD and H-REMD simulation scheme greatly accelerates the configurational sampling of the rotameric states of the side chains around the binding pocket, thereby improving the convergence of the FEP computations.     

Reference-Quality Free Energy Barriers in Catalysis from Machine Learning Thermodynamic Perturbation Theory

For the first time, we report calculations of the free energies of activation of cracking and isomerization reactions of alkenes that combine several different electronic structure methods with molecular dynamics simulations. We demonstrate that the use of a high level of theory (here Random Phase Approximation—RPA) is necessary to bridge the gap between experimental and computed values. These transformations, catalyzed by zeolites and proceeding via cationic intermediates and transition states, are building blocks of many chemical transformations for valorization of long chain paraffins originating, e.g., from plastic waste, vegetable oils, Fischer–Tropsch waxes or crude oils. Compared with the free energy barriers computed at the PBE+D2 production level of theory via constrained ab initio molecular dynamics, the barriers computed at the RPA level by the application of Machine Learning thermodynamic Perturbation Theory (MLPT) show a significant decrease for isomerization reaction and an increase of a similar magnitude for cracking, yielding an unprecedented agreement with the results obtained by experiments and kinetic modeling.     

Multiple‐Replica Exchange with Information Retrieval

The parallel tempering simulation method was recently extended to allow for possible exchanges between non‐adjacent replicas. We introduce a multiple‐exchange variant which naturally incorporates the information from all replicas when calculating statistical averages, building on the related virtual‐move method of Coluzza and Frenkel (ChemPhysChem 2005 , 6, 1779). The method is extensively tested on three model systems, namely, a Lennard‐Jones cluster exhibiting a finite size phase transition, the Lennard‐Jones fluid, and the 2D ferromagnetic Ising model. In all cases, the present method performs significantly better and converges faster than conventional parallel tempering Monte Carlo simulations. The standard deviations are also systematically decreased with respect to virtual moves.     

Absolute Protein Binding Free Energy Simulations for Ligands with Multiple Poses,a Thermodynamic Path That Avoids Exhaustive Enumeration of the Poses

We propose a free energy calculation method for receptor–ligand binding, which have multiple binding poses that avoids exhaustive enumeration of the poses. For systems with multiple binding poses, the standard procedure is to enumerate orientations of the binding poses, restrain the ligand to each orientation, and then, calculate the binding free energies for each binding pose. In this study, we modify a part of the thermodynamic cycle in order to sample a broader conformational space of the ligand in the binding site. This modification leads to more accurate free energy calculation without performing separate free energy simulations for each binding pose. We applied our modification to simple model host–guest systems as a test, which have only two binding poses, by using a single decoupling method (SDM) in implicit solvent. The results showed that the binding free energies obtained from our method without knowing the two binding poses were in good agreement with the benchmark results obtained by explicit enumeration of the binding poses. Our method is applicable to other alchemical binding free energy calculation methods such as the double decoupling method (DDM) in explicit solvent. We performed a calculation for a protein–ligand system with explicit solvent using our modified thermodynamic path. The results of the free energy simulation along our modified path were in good agreement with the results of conventional DDM, which requires a separate binding free energy calculation for each of the binding poses of the example of phenol binding to T4 lysozyme in explicit solvent. © 2019 Wiley Periodicals, Inc.     

Computation of Absolute Hydration and Binding Free Energy with Free Energy Perturbation Distributed Replica-Exchange Molecular Dynamics (FEP/REMD)

Distributed Replica (REPDSTR) is a powerful parallelization technique enabling simulations of a group of replicas in a parallel/parallel fashion, where each replica is distributed to different nodes of a large cluster [Theor. Chem. Acc. 109: 140 (2003)]. Here, we use the framework provided by REPDSTR to combine a staged free energy perturbation protocol with different values of the thermodynamic coupling parameters with replica-exchange molecular dynamics (FEP/REMD. The structure of REPDSTR, which allows multiple parallel input/output (I/O), facilitates the treatment of replica-exchange to couple the N window simulations. As a result, each of the N synchronous window simulations benefit from the sampling carried out by the N-1 others. As illustrative examples of the FEP/REMD strategy, calculations of the absolute hydration and binding free energy of small molecules were performed using the biomolecular simulation program CHARMM adapted for the IBM Blue Gene/P platform. The computations show that a FEP/REMD strategy significantly improves the sampling and accelerate the convergence of absolute free energy computations.     

Laplacian of the Hamiltonian Kinetic Energy Density as an Indicator of Binding and Weak Interactions

The kinetic energy is the center of a controversy between two opposite points of view about its role in the formation of a chemical bond. One school states that a lowering of the kinetic energy associated with electron delocalization is the key stabilization mechanism of covalent bonding. In contrast, the opposite school holds that a chemical bond is formed by a decrease in the potential energy due to a concentration of electron density within the binding region. In this work, a topographic analysis of the Hamiltonian Kinetic Energy Density (KED) and its laplacian is presented to gain more insight into the role of the kinetic energy within chemical interactions. This study is focused on atoms, diatomic and organic molecules, along with their dimers. In addition, it is shown that the laplacian of the Hamiltonian KED exhibits a shell structure in atoms and that their outermost shell merge when a molecule is formed. A covalent bond is characterized by a concentration of kinetic energy, potential energy and electron densities along the internuclear axis, whereas a charge-shift bond is characterized by a fusion of external concentration shells and a depletion in the bonding region. In the case of weak intermolecular interactions, the external shell of the molecules merge into each other resulting in an intermolecular surface comparable to that obtained by the Non-covalent interaction (NCI) analysis.     

Revision of AMBER Torsional Parameters for RNA Improves Free Energy Predictions for Tetramer Duplexes with GC and iGiC Base Pairs

All-atom force fields are important for predicting thermodynamic, structural, and dynamic properties of RNA. In this paper, results are reported for thermodynamic integration calculations of free energy differences of duplex formation when CG pairs in the RNA duplexes r(CCGG)(2), r(GGCC)(2), r(GCGC)(2), and r(CGCG)(2) are replaced by isocytidine-isoguanosine (iCiG) pairs. Agreement with experiment was improved when ε/ζ, α/γ, β, and χ torsional parameters in the AMBER99 force field were revised on the basis of quantum mechanical calculations. The revised force field, AMBER99TOR, brings free energy difference predictions to within 1.3, 1.4, 2.3, and 2.6 kcal/mol at 300 K, respectively, compared to experimental results for the thermodynamic cycles of CCGG → iCiCiGiG, GGCC → iGiGiCiC, GCGC → iGiCiGiC, and CGCG → iCiGiCiG. In contrast, unmodified AMBER99 predictions for GGCC → iGiGiCiC and GCGC → iGiCiGiC differ from experiment by 11.7 and 12.6 kcal/mol, respectively. In order to test the dynamic stability of the above duplexes with AMBER99TOR, four individual 50 ns molecular dynamics (MD) simulations in explicit solvent were run. All except r(CCGG)(2) retained A-form conformation for ≥82% of the time. This is consistent with NMR spectra of r(iGiGiCiC)(2), which reveal an A-form conformation. In MD simulations, r(CCGG)(2) retained A-form conformation 52% of the time, suggesting that its terminal base pairs may fray. The results indicate that revised backbone parameters improve predictions of RNA properties and that comparisons to measured sequence dependent thermodynamics provide useful benchmarks for testing force fields and computational methods.     

蛋白质相互作用: 界面分析, 结合自由能计算与相互作用设计

蛋白质相互作用在生命活动中起着重要作用. 研究蛋白质间相互作用的本质有助于了解生命活动中这些基本单元的作用. 本文主要综述了近期蛋白质相互作用研究的进展, 包括蛋白质相互作用界面的基本性质, 蛋白质结合自由能的计算方法, 不同相互作用在蛋白质结合/解离中的角色和差异, 以及上述知识在蛋白质相互作用设计中的应用. 蛋白质相互作用界面的特性, 例如界面大小、保守性以及结构的动态性质, 使得具有生物功能的蛋白质相互作用界面区别于非特异性的晶体堆积界面. 生物功能界面的一个重要结构特征是界面上存在着关键残基以及相对独立的相互作用模块. 利用多种方法, 如MM-PBSA、统计平均势以及不同的相互作用自由能模型, 可以在不同的精度上计算蛋白质相互作用自由能. 利用蛋白质相互作用界面的特点, 从不同的角度进行蛋白质相互作用对的设计与改造, 近年来已经有了不少成功的例子, 但还存在着很大的挑战. 我们认为在今后的蛋白质相互作用设计中, 考虑各种因素对蛋白质结合与解离的动力学过程的影响将有助于提高人类控制蛋白质相互作用的能力.     

不同形貌Ag_3PO_4的表面自由能及原位光催化过程的热力学和动力学参数

<正>目前,在利用表面能研究纳米催化剂形貌效应时,多以催化剂表面原子排布或不同晶面的夹角与密勒指数间的关系来定性描述表面能高低[1,2],或利用密度泛函理论对表面能进行理论计算[3,4],但其数值可能会因计算方法和模型不符而与真实值存在偏差,而表面能的实验测定至今仍无精确、可靠、通用的方法[5].因此发展一种准确获取表面能的普适性方法并与材料的催化性能相关联非常必要.微量热技术能对体系变化过程进行实时监测,通过体系变化过程的热力学和动力学参数获取纳米材料的     

Alchemical Grid Dock (AlGDock): Binding Free Energy Calculations between Flexible Ligands and Rigid Receptors

Alchemical Grid Dock (AlGDock) is open-source software designed to compute the binding potential of mean force—the binding free energy between a flexible ligand and a rigid receptor—for a small organic ligand and a biological macromolecule. Multiple BPMFs can be used to rigorously compute binding affinities between flexible partners. AlGDock uses replica exchange between thermodynamic states at different temperatures and receptor–ligand interaction strengths. Receptor–ligand interaction energies are represented by interpolating precomputed grids. Thermodynamic states are adaptively initialized and adjusted on-the-fly to maintain adequate replica exchange rates. In demonstrative calculations, when the bound ligand is treated as fully solvated, AlGDock estimates BPMFs with a precision within 4 kT in 65% and within 8 kT for 91% of systems. It correctly identifies the native binding pose in 83% of simulations. Performance is sometimes limited by subtle differences in the important configuration space of sampled and targeted thermodynamic states. © 2019 Wiley Periodicals, Inc.     

A Thermodynamic Study on the Binding of Human Serum Albumin with Lanthanum Ion

Thermodynamics of the interaction between lanthanum(III) ion, La³⁺, and human serum albumin (HSA), was investigated at pH 7.0 and 300 K in Tris‐HCl buffer by isothermal titration calorimetry. A solvation model was used to reproduce the enthalpies of HSA interaction with La³⁺. The solvation parameters recovered from our new model were attributed to the structural change of HSA and its biological activity. The interaction of HSA with La³⁺ showed a set of two binding sites with negative cooperativity.     

Calculation of the Absolute Free Energy of Binding and Related Entropies with the HSMD-TI Method: The FKBP12-L8 Complex

The hypothetical scanning molecular dynamics (HSMD) method is used here for calculating the absolute free energy of binding, ΔA(0) of the complex of the protein FKBP12 with the ligand SB2 (also denoted L8) - a system that has been studied previously for comparing the performance of different methods. Our preliminary study suggests that considering long-range electrostatics is imperative even for a hydrophobic ligand such as L8. Therefore the system is modeled by the AMBER force field using Particle Mesh Ewald (PME). HSMD consists of three stages applied to both the ligand-solvent and ligand-protein systems. (1) A small set of system configurations (frames) is extracted from an MD trajectory. (2) The entropy of the ligand in each frame is calculated by a reconstruction procedure. (3) The contribution of water and protein to ΔA(0) is calculated for each frame by gradually increasing the ligand-environment interactions from zero to their full value using thermodynamic integration (TI). Unlike the conventional methods, the structure of the ligand is kept fixed during TI, and HSMD is thus free from the end-point problem encountered with the double annihilation method (DAM); therefore, the need for applying restraints is avoided. Furthermore, unlike the conventional methods, the entropy of the ligand and water is obtained directly as a byproduct of the simulation. In this paper, in addition to the difference in the internal entropies of the ligand in the two environments, we calculate for the first time the external entropy of the ligand, which provides a measure for the size of the active site. We obtain ΔA(0) = -10.7 ±1.0 as compared to the experimental values -10.9 and -10.6 kcal/mol. However, a protein/water system treated by periodic boundary conditions grows significantly with increasing protein size and the computation of ΔA(0) would become expensive by all methods. Therefore, we also apply HSMD to FKBP12-L8 described by the GSBP/SSBP model of Roux's group (implemented in the software CHARMM) where only part of the protein and water around the active site are considered and long-range electrostatic effects are taken into account. For comparison this model was also treated by the double decoupling method (DDM). The two methods have led to comparable results for ΔA(0) which are somewhat lower than the experimental value. The ligand was found to be more confined in the active site described by GSBP/SSBP than by PME where its entropy in solvent is larger than in the active site by 1.7 and by 5.5 kcal/mol, respectively.     

Structural and Thermodynamic Analysis of Daunomycin Binding with Desoxyhexanucleotides with Different Base Sequences by NMR Spectroscopy

Complexation of the anthracycline antibiotic daunomycin (DAU) with self-complementary desoxyhexanucleotides with different base sequences in the chain, 5-d(CGTACG) and 5-d(CGCGCG), is studied in aqueous solution. Homonuclear 2D-¹H NMR spectroscopy (TOCSY and NOESY) and heteronuclear ¹H–³¹P NMR spectroscopy (HMBC) are used for complete assignment of the nonexchangeable proton resonances and the phosphorus signals and for qualitative determination of preferable DAU binding sites, respectively. Daunomycin is preferably intercalated into the first d(CG)-site of each hexanucleotide, and the aminosaccharide ring of DAU lies in the minor groove of the hexamer duplex, partly overlapping the third base pair. Quantitative analysis of DAU complexation with desoxyhexanucleotides was done by analyzing the concentration and temperature dependences of the DAU proton chemical shifts; the equilibrium reaction constants and the thermodynamic parameters of the formation of the 1:1, 1:2, 2:1, and 2:2 DAU complexes with hexamers, as well as the limiting values of the DAU proton chemical shifts were determined in aqueous solution. The antibiotic preferably binds with the triplet sections of the nucleotide sequence containing two neighboring CG-pairs of nitrogen bases flanked by the AT-pair in 5-d(CGTACG) compared to the triplets consisting of a sequence of three CG-pairs in the 5-d(CGCGCG) hexamer. The specific binding of daunomycin depends on the character of the hydrophobic interactions of the aminosaccharide ring of the antibiotic in the minor groove of the DNA double helix. Calculations have been carried out to determine the most probable spatial structures of the 1:2 DAU–desoxyhexanucleotide complexes; the results agree well with X-ray diffraction data.     

多肽抑制剂抑制淀粉质多肽42构象转换的分子动力学模拟和结合自由能计算

应用分子动力学模拟和结合自由能计算方法研究了多肽抑制剂KLVFF、VVIA和LPFFD抑制淀粉质多肽42 (Aβ42)构象转换的分子机理. 结果表明, 三种多肽抑制剂均能够有效抑制Aβ42的二级结构由α-螺旋向β-折叠的构象转换. 另外, 多肽抑制剂降低了Aβ42分子内的疏水相互作用, 减少了多肽分子内远距离的接触, 有效抑制了Aβ42的疏水塌缩, 从而起到稳定其初始构象的作用. 这些抑制剂与Aβ42之间的疏水和静电相互作用(包括氢键)均有利于它们抑制Aβ42的构象转换. 此外, 抑制剂中的带电氨基酸残基可以增强其和Aβ42之间的静电相互作用(包括氢键), 并降低抑制剂之间的聚集, 从而大大增强对Aβ42构象转换的抑制能力. 但脯氨酸的引入会破坏多肽的线性结构, 从而大大降低其与Aβ42 之间的作用力. 上述分子模拟的结果揭示了多肽抑制剂KLVFF、VVIA和LPFFD抑制Aβ42构象转换的分子机理, 对于进一步合理设计Aβ的高效短肽抑制剂具有非常重要的理论指导意义.     

微量热法研究碱土金属离子与18-冠-6结合的热力学常数

通过微量热法研究了18-冠-6与碱土金属离子结合过程. 根据反应后的平衡关系, 求出反应平衡常数、焓、熵和吉布斯自由能变化. 金属离子由原来溶剂化的游离态变为与18-冠-6以静电力结合, 体系的能量降低, 放出了热量;同时由原来的自由运动变为被限制在冠醚孔穴内, 自由度减小, 熵值减少. 冠醚内的氧原子能够紧紧地围绕着半径为138 pm的Ba2+, 而离子半径较小的其它离子, 围绕得不够紧密, 与氧原子之间的作用力变小. 从Be2+、Mg2+、Ca2+、Sr2+到Ba2+, 与冠醚结合的放热增加, 平衡常数也逐渐增大.     

A New Approach for Thermodynamic Studies on Binding of Some Metal Ions with Human Growth Hormone

A thermodynamic study on the interaction between Mg²⁺, Ni²⁺ and Co²⁺ ions (M²⁺) and human growth hormone (hGH) was made at 27 °C in aqueous NaCl solutions using isothermal titration calorimetry. Gholamreza Rezaei Behbehani’s solvation model (GRB) was used to model the enthalpies of M²⁺+hGH interactions over the studied range of metal ion concentrations. The solvation parameters derived from the solvation model were attributed to a structural change of hGH due to its interactions with the metal ion.     

加压毛细管电色谱手性配体交换法对游离亮氨酸异构体分离的研究

以反相C18为固定相,通过对流动相中的手性选择体、流动相组成与pH值及操作电压等条件的考察,成功建立了毛细管电色谱（CEC）手性配体交换法分离游离亮氨酸异构体的方法。优化的色谱条件：C18反相毛细管柱;以2mmol／L天冬氨酰苯丙氨酸甲酯为手性选择体,Cu^2＋为配体离子,流动相为20mmol／L乙酸缓冲液（pH5．5）,5％（体积分数）的甲醇为有机改性剂;流速为0．02mL／min;操作压力为6．9&#215;10^6 Pa;检测波长为UV 280nm;工作电压＋5kV;室温。两异构体可在15min内基线分离,检出限为10μmol／L（S／N=3）,并在0．05～1．5mmol／L范围内呈良好的线性关系,日让、日内精密度均小于3．5％。