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.     

副本交换分子动力学

副本交换分子动力学(REMD)是一种广泛应用于蛋白质功能性构象变化模拟及相应自由能计算的增强型采样算法。由于REMD理论严格且采样效率高,近年来备受关注,尤其是针对传统REMD方法的发展和优化,显著提高了REMD的采样效率,拓展了其应用范围。但是各种REMD新型方法的最佳适用范围也存在较大区别,使得如何选用合适的REMD方法成为实际应用的难题和挑战。因此,有必要对各种REMD方法及其应用进行阐述,深入比较各方法的优缺点及其实际应用体系。本综述从REMD的原理出发,回顾了近年来各类REMD方法的变形策略,以助于对REMD方法的理解、应用和继续改进。     

Improved Binding Free Energy Predictions from Single-Reference Thermodynamic Integration Augmented with Hamiltonian Replica Exchange

Reliable predictions of relative binding free energies are essential in drug discovery, where chemists modify promising compounds with the aim of increasing binding affinity. Conventional Thermodynamic Integration (TI) approaches can estimate corresponding changes in binding free energies, but suffer from inadequate sampling due to ruggedness of the molecular energy surfaces. Here, we present an improved TI strategy for computing relative binding free energies of congeneric ligands. This strategy employs a specific, unphysical single-reference (SR) state and Hamiltonian replica exchange (HREX) to locally enhance sampling. We then apply this strategy to compute relative binding free energies of twelve ligands in the L99A mutant of T4 Lysozyme. Besides the ligands, our approach enhances hindered rotations of the important V111, as well as V87 and L118 sidechains. Concurrently, we devise practical strategies to monitor and improve HREX-SRTI efficiency. Overall, the HREX-SRTI results agree well (R(2) = 0.76, RMSE = 0.3 kcal/mol) with available experimental data. When optimized for efficiency, the HREX-SRTI precision matches that of experimental measurements.     

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.     

β发卡多肽Trpzip4折叠的副本交换分子动力学模拟

采用副本交换分子动力学对β发卡Trpzip4重折叠进行研究,结果表明,构象空间存在较大能垒时,副本交换分子动力学(REMD)表现出比分子动力学(MD)更优的抽样效率.288 K下,Trpzip4势能、骨架均方根偏差、溶剂可及表面积(SASA)在REMD中呈逐渐降低趋势.采样得到两种特定形态的低势能构象:β发卡和螺旋-卷曲.β发卡中,第3组氢键Asp46容易与Thr49形成氢键;转角Thr49的羟基(OH)倾向于与Asp46的羧基(COO-)或主链上的C=O形成氢键,强化了Asp46与Thr49的相互作用,导致转角形成折回弯度;与Thr49相比,Thr51的羟基倾向与水分子作用,甲基朝内,因此与Asp46的侧链作用微小.螺旋-卷曲中,吲哚环-甲基为主要疏水形式.Trpzip4在折叠成β发卡过程中,转角氢键影响整个β发卡构象的形成.转角氢键具有距离优势和强侧链相互作用,最先形成.     

谷氨酰胺结合蛋白的分子动力学模拟和自由能计算

谷氨酰胺结合蛋白(Glutamine-binding protein, GlnBp)是大肠杆菌透性酶系统中一个细胞外液底物专一性结合蛋白, 对于细胞外液中谷氨酰胺(Gln)的运输和传递至关重要. 本文运用分子动力学(Molecular dynamics, MD)模拟采样, 考察了GlnBp关键残基与底物Gln之间的相互作用和GlnBp两条铰链的功能差别; 并采用MM-PBSA方法计算了GlnBp与底物Gln的结合自由能. 结果表明: Ph13, Phe50, Thr118和Ile69与底物Gln的范德华相互作用和Arg75, Thr70, Asp157, Gly68, Lys115, Ala67, His156与底物Gln的静电相互作用是结合Gln的主要推动力; 复合物的铰链区85～89柔性大, 对构象开合提供了结构基础; 而铰链区181～185柔性小, 其作用更多是在功能上把底物Gln限制在口袋中; 自由能预测值与实验值吻合. 本研究很好地解释了GlnBp结构与功能的关系, 为进一步了解GlnBp的开合及转运Gln的机制提供了重要的结构信息.     

关于副本交换分子动力学模拟复杂化学反应的研究

本文研究用温度副本交换分子动力学(T-REMD)和哈密顿副本交换分子动力学(H-REMD)方法模拟复杂化学反应的问题。使用具有不同活化能和反应能的简单置换反应模型,我们检验了上述两种方法用来预测反应平衡产物的效率和应用范围。T-REMD方法对具有适度活化能(约< 20 kcal·mol^-1)或者反应能量(< 3 kcal·mol^-1)的放热反应是有效的。由于在相空间的不完整采样,对于同时具有高活化能和反应能量的反应其模拟效率有严重障碍,并且对于吸热反应问题更为显着。另一方面,H-REMD对一系列具有不同活化能的反应能的模型表现出色,与T-REMD相比,H-REMD可以使用更少的副本获得优异的结果。     

The Thermodynamics of Folding of a β Hairpin Peptide Probed Through Replica Exchange Molecular Dynamics Simulations

Peptides that possess a well defined native state are ideal model systems to study the folding of proteins. They possess many of the complexities of larger proteins, yet their small size renders their study computationally tractable. Recent advances in sampling techniques, including replica exchange molecular dynamics, now permit a full characterization of the thermodynamics of folding of small peptides. These simulations not only yield insight into the folding of larger proteins, but equally importantly, they allow, through comparison with experiment, an objective test of the accuracy of force fields, water models and of different numerical schemes for dealing with electrostatic interactions. In this account, we present a molecular dynamics simulation of a small β-hairpin peptide using the replica exchange algorithm and illustrate how this enhanced sampling scheme enables a thorough characterization of the native and unfolded states, and sheds new light into its folding mechanism.     

烷基芳基磺酸盐的分子动力学模拟与自由能微扰计算

为研究不同结构的表面活性剂分子在溶液中胶束化能力的差异, 采用分子动力学方法模拟三种烷基芳基磺酸盐在真空和水溶液环境下的结构与相互作用. 利用自由能微扰(FEP)方法计算了水合自由能, 发现与用传统热力学表面张力法测定自制的烷基芳基磺酸盐结果一致. 研究表明: 烷基芳基磺酸盐在水溶液中的胶束化过程是自发进行的, 随着分子结构中芳环向长烷基链中间位置移动, 胶束化能力和胶束稳定性均下降; 疏水基周围水分子的“冰山结构”会影响胶束的稳定性, 而水分子中氢键的生存周期是反映冰山结构变化的重要指标; 同时, 亲水基与水分子间形成氢键的数目会增强或减弱分子脱离胶束体的趋势, 从而影响胶束结构的稳定性.     

葡萄糖苷酶抑制剂作用机理的分子动力学模拟和自由能计算

含有锍离子的葡萄糖苷酶抑制剂如kotalanol (SK)和它除去磺酸基团后的衍生物(DSK), 是潜在的毒副作用较小的治疗II 型糖尿病的候选药物. α-葡萄糖苷酶抑制活性实验显示, DSK活性比SK略高, 而将二者环上的S原子替换成NH后(分别称为DSN和SN), DSN的活性要比SN高1500倍左右. 本文用分子动力学模拟, 结合自由能计算和自由能分解的方法对上述四个抑制剂的作用机理进行了研究. 研究结果表明活性的巨大差异是由NH基团取代效应和磺酸基团立体效应共同作用的结果, 由于N―C键长比S―C键长短, NH基团取代导致烷基链的翻转, 同时, 磺酸基团限制了链的翻转, 因此改变了抑制剂的结合模式. 计算结果与实验基本一致.本文的研究结果有助于进一步理解含锍离子的葡萄糖苷酶抑制剂的结合机理, 并为设计更有潜力的葡萄糖苷酶抑制剂提供了有价值的信息.     

FKBP12蛋白与其抑制剂结合的分子动力学模拟和结合自由能计算

FKBP12 (FK506-binding protein-12)是一种具有神经保护和促神经再生作用的蛋白. 采用分子动力学模拟取样, 运用MM-GBSA方法计算了FKBP12和3个抑制剂(GPI-1046, 308和107)的绝对结合自由能, GPI-1046的结合能最小, 308小于107的结合能. 通过能量分解的方法考察了FKBP12蛋白的主要残基与抑制剂之间的相互作用和识别, 计算结果表明: 3个抑制剂具有相似的结合模式, Ile56和Tyr82主要表现为氢键作用, Tyr26, Phe46, Val55, Ile56, Trp59, Tyr82, Tyr87和Phe99形成疏水作用区. 计算结果和实验结果吻合.     

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.     

Dissipative Particle Dynamics Simulations of Polymer Brushes: Comparison with Molecular Dynamics Simulations

Summary: The structure of polymer brushes is investigated by dissipative particle dynamics (DPD) simulations that include explicit solvent particles. With an appropriate choice of the DPD interaction parameters , we obtain good agreement with previous molecular dynamics (MD) results where the good solvent behavior has been modeled by an effective Lennard–Jones potential. The present results confirm that DPD simulation techniques can be applied for large length scale simulations of polymer brushes. A relation between the different length scales and is established.

Polymer brush at a solid–liquid interface.     

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.     

新型二氟甲基磷酸类酪氨酸蛋白磷酸酯酶1B抑制剂的分子动力学模拟和结合自由能计算

通过分子对接建立了一系列含二氟甲基磷酸基团(DFMP)或二氟甲基硫酸基团(DFMS)的抑制剂与酪氨酸蛋白磷酸酯酶1B(PTP1B)的相互作用模式, 并通过1 ns的分子动力学模拟和molecular mechanics/generalized Born surface area (MM/GBSA)方法计算了其结合自由能. 计算获得的结合自由能排序和抑制剂与靶酶间结合能力排序一致; 通过基于主方程的自由能计算方法, 获得了抑制剂与靶酶残基间相互作用的信息, 这些信息显示DFMP/DFMS基团的负电荷中心与PTP1B的221位精氨酸正电荷中心之间的静电相互作用强弱决定了此类抑制剂的活性, 进一步的分析还显示位于DFMP/DFMS基团中的氟原子或其他具有适当原子半径的氢键供体原子会增进此类抑制剂与PTP1B活性位点的结合能力.     

Using Selectively Applied Accelerated Molecular Dynamics to Enhance Free Energy Calculations

Accelerated molecular dynamics (aMD) has been shown to enhance conformational space sampling relative to classical molecular dynamics; however, the exponential reweighting of aMD trajectories, which is necessary for the calculation of free energies relating to the classical system, is oftentimes problematic, especially for systems larger than small poly peptides. Here, we propose a method of accelerating only the degrees of freedom most pertinent to sampling, thereby reducing the total acceleration added to the system and improving the convergence of calculated ensemble averages, which we term selective aMD. Its application is highlighted in two biomolecular cases. First, the model system alanine dipeptide is simulated with classical MD, all-dihedral aMD, and selective aMD, and these results are compared to the infinite sampling limit as calculated with metadynamics. We show that both forms of aMD enhance the convergence of the underlying free energy landscape by 5-fold relative to classical MD; however, selective aMD can produce improved statistics over all-dihedral aMD due to the improved reweighting. Then we focus on the pharmaceutically relevant case of computing the free energy of the decoupling of oseltamivir in the active site of neuraminidase. Results show that selective aMD greatly reduces the cost of this alchemical free energy transformation, whereas all-dihedral aMD produces unreliable free energy estimates.     

Protecting High Energy Barriers: A New Equation to Regulate Boost Energy in Accelerated Molecular Dynamics Simulations

Molecular dynamics (MD) is one of the most common tools in computational chemistry. Recently, our group has employed accelerated molecular dynamics (aMD) to improve the conformational sampling over conventional molecular dynamics techniques. In the original aMD implementation, sampling is greatly improved by raising energy wells below a predefined energy level. Recently, our group presented an alternative aMD implementation where simulations are accelerated by lowering energy barriers of the potential energy surface. When coupled with thermodynamic integration simulations, this implementation showed very promising results. However, when applied to large systems, such as proteins, the simulation tends to be biased to high energy regions of the potential landscape. The reason for this behavior lies in the boost equation used since the highest energy barriers are dramatically more affected than the lower ones. To address this issue, in this work, we present a new boost equation that prevents oversampling of unfavorable high energy conformational states. The new boost potential provides not only better recovery of statistics throughout the simulation but also enhanced sampling of statistically relevant regions in explicit solvent MD simulations.     

Molecular Dynamics Simulations of Proteins with Chemically Modified Disulfide Bonds

Proteins that are used as therapeutic drugs act in the extracellular microenvironment. They usually have a small number of intramolecular disulfide bonds to help maintain their tertiary structure in the vascular circulation. In general, most cysteine residues are part of a disulfide bond with free sulfhydrals being uncommon. We have studied whether the site-specific chemical reduction of disulfides and the incorporation of a 3-carbon methylene bridge between the cysteines in interferon-α 2a would change the structure of this protein. Bridging of both of the disulfide bonds of interferon-α 2a was studied using two different molecular simulation protocols: (1) molecular dynamics, and (2) stochastic dynamics. We have shown that the disulfide bonds in interferon-α 2a can be reduced and chemically modified without significantly altering the tertiary structure of the protein. This offers the novel possibility of chemically modifying therapeutically important proteins without affecting their biological properties.     

基于分子动力学模拟和连续介质模型的自由能计算方法*

近些年,基于分子动力学模拟和连续介质模型的自由能计算方法受到了越来越多的关注,其中MM/PBSA就是最具代表性的方法.在MM/PBSA中,体系的焓变采用分子力学(MM)的方法计算得到;溶剂效应中极性部分对自由能的贡献通过解Poisson-Boltzmann(PB)方程的方法计算得到;溶液效应中非极性部分对自由能的贡献则通过分子表面积(SA)计算得到.本文结合我们科研组的工作,就近几年MM/PBSA方法的最新进展做了较为详细的阐述,同时对MM/PBSA的发展前景进行了展望.     

Free Energy Calculations using a Swarm‐Enhanced Sampling Molecular Dynamics Approach

Free energy simulations are an established computational tool in modelling chemical change in the condensed phase. However, sampling of kinetically distinct substates remains a challenge to these approaches. As a route to addressing this, we link the methods of thermodynamic integration (TI) and swarm‐enhanced sampling molecular dynamics (sesMD), where simulation replicas interact cooperatively to aid transitions over energy barriers. We illustrate the approach by using alchemical alkane transformations in solution, comparing them with the multiple independent trajectory TI (IT‐TI) method. Free energy changes for transitions computed by using IT‐TI grew increasingly inaccurate as the intramolecular barrier was heightened. By contrast, swarm‐enhanced sampling TI (sesTI) calculations showed clear improvements in sampling efficiency, leading to more accurate computed free energy differences, even in the case of the highest barrier height. The sesTI approach, therefore, has potential in addressing chemical change in systems where conformations exist in slow exchange.