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

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

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.     

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.     

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

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

Computational Analysis of Triazole-Based Kojic Acid Analogs as Tyrosinase Inhibitors by Molecular Dynamics and Free Energy Calculations

Molecular docking, molecular dynamics (MD) simulations and the linear interaction energy (LIE) method were used here to predict binding modes and free energy for a set of 1,2,3-triazole-based KA analogs as potent inhibitors of Tyrosinase (TYR), a key metalloenzyme of the melanogenesis process. Initially, molecular docking calculations satisfactorily predicted the binding mode of evaluated KA analogs, where the KA part overlays the crystal conformation of the KA inhibitor into the catalytic site of TYR. The MD simulations were followed by the LIE method, which reproduced the experimental binding free energies for KA analogs with an r² equal to 0.97, suggesting the robustness of our theoretical model. Moreover, the van der Waals contributions performed by some residues such as Phe197, Pro201, Arg209, Met215 and Val218 are responsible for the binding recognition of 1,2,3-triazole-based KA analogs in TYR catalytic site. Finally, our calculations provide suitable validation of the combination of molecular docking, MD, and LIE approaches as a powerful tool in the structure-based drug design of new and potent TYR inhibitors.     

On the Use of Accelerated Molecular Dynamics to Enhance Configurational Sampling in Ab Initio Simulations

We have implemented the accelerated molecular dynamics approach (Hamelberg, D.; Mongan, J.; McCammon, J. A. J. Chem. Phys. 2004, 120 (24), 11919) in the framework of ab initio MD (AIMD). Using three simple examples, we demonstrate that accelerated AIMD (A-AIMD) can be used to accelerate solvent relaxation in AIMD simulations and facilitate the detection of reaction coordinates: (i) We show, for one cyclohexane molecule in the gas phase, that the method can be used to accelerate the rate of the chair-to-chair interconversion by a factor of ～1 × 10(5), while allowing for the reconstruction of the correct canonical distribution of low-energy states; (ii) We then show, for a water box of 64 H(2)O molecules, that A-AIMD can also be used in the condensed phase to accelerate the sampling of water conformations, without affecting the structural properties of the solvent; and (iii) The method is then used to compute the potential of mean force (PMF) for the dissociation of Na-Cl in water, accelerating the convergence by a factor of ～3-4 compared to conventional AIMD simulations.(2) These results suggest that A-AIMD is a useful addition to existing methods for enhanced conformational and phase-space sampling in solution. While the method does not make the use of collective variables superfluous, it also does not require the user to define a set of collective variables that can capture all the low-energy minima on the potential energy surface. This property may prove very useful when dealing with highly complex multidimensional systems that require a quantum mechanical treatment.     

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.     

In Silico Study towards Repositioning of FDA-Approved Drug Candidates for Anticoronaviral Therapy: Molecular Docking,Molecular Dynamics and Binding Free Energy Calculations

The SARS-CoV-2 targets were evaluated for a set of FDA-approved drugs using a combination of drug repositioning and rigorous computational modeling methodologies such as molecular docking and molecular dynamics (MD) simulations followed by binding free energy calculations. Six FDA-approved drugs including, Ouabain, Digitoxin, Digoxin, Proscillaridin, Salinomycin and Niclosamide with promising anti-SARS-CoV-2 activity were screened in silico against four SARS-CoV-2 proteins—papain-like protease (PLpro), RNA-dependent RNA polymerase (RdRp), SARS-CoV-2 main protease (Mpro), and adaptor-associated kinase 1 (AAK1)—in an attempt to define their promising targets. The applied computational techniques suggest that all the tested drugs exhibited excellent binding patterns with higher scores and stable complexes compared to the native protein cocrystallized inhibitors. Ouabain was suggested to act as a dual inhibitor for both PLpro and Mpro enzymes, while Digitoxin bonded perfectly to RdRp. In addition, Salinomycin targeted PLpro. Particularly, Niclosamide was found to target AAK1 with greater affinity compared to the reference drug. Our study provides comprehensive molecular-level insights for identifying or designing novel anti-COVID-19 drugs.     

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

谷氨酰胺结合蛋白(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的机制提供了重要的结构信息.     

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

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

重要性采样方法与自由能计算

分子动力学模拟与自由能计算已经在化学、生物学与材料学等领域得到广泛的应用。然而,由于在传统分子动力学模拟的时间尺度内,体系很难跨越较高的自由能能垒,在相空间内的采样大大受限,采样困难使自由能计算难以收敛。增强采样是解决这一问题的有效途径,重要性采样方法就是其中一类。本文综述了四种广泛应用的重要性采样方法--伞状采样方法、metadynamics方法、自适应偏置力方法和温度加速分子动力学方法的原理和进展,其中重点概述了自适应偏置力方法的最新发展--扩展自适应偏置力方法和扩展广义自适应偏置力方法,并对这四种重要性采样方法的优缺点进行了比较。最后,讨论和展望了重要性采样与自由能计算方法面临的挑战和前景,并提出了对自适应偏置力方法可能的改进,如与加速分子动力学（aMD）或弦方法结合以提高在高维度空间中的采样效率。     

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

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

Molecular Dynamics Simulation of Aqueous Solutions Using Interaction Energy Components: Application to the Solvation Gibbs Energy

Molecular dynamics simulations of aqueous solutions of the solutes acetamide (AcNH₂), acetic acid (AcOH), and acetaldehyde (AcH) were made using Lennard–Jones 12-6-1 potentials to describe the solute–solvent interactions. The Morokuma decomposition scheme and the ESIE solute atomic charges were used to reproduce the exchange, polarization, and electrostatic components of the solute–water interaction energy. A nonlinear perturbation was incorporated into the “slow-growth” technique in order to improve the results for the solvation Gibbs energy that were found to be in agreement with the available experimental and theoretical values.     

Inhibition Mechanism of Hydroxyproline-like Small Inhibitors to Disorder HIF-VHL Interaction by Molecular Dynamic Simulations and Binding Free Energy Calculations

Protein-protein interactions are vital for a wide range of biological processes. The interactions between the hypoxia-inducible factor and von Hippel Lindau (VHL) are attractive drug targets for ischemic heart disease. In order to disrupt this interaction, the strategy to target VHL binding site using a hydroxyproline-like (pro-like) small molecule has been reported. In this study, we focused on the inhibition mechanism between the pro-like inhibitors and the VHL protein, which were investigated via molecular dynamics simulations and binding free energy calculations. It was found that pro-like inhibitors showed a strong binding affinity toward VHL. Binding free energy calculations and free energy decompositions suggested that the modification of various regions of pro-like inhibitors may provide useful information for future drug design.     

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.     

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形成疏水作用区. 计算结果和实验结果吻合.     

Conformational Transitions and Convergence of Absolute Binding Free Energy Calculations

The Binding Energy Distribution Analysis Method (BEDAM) is employed to compute the standard binding free energies of a series of ligands to a FK506 binding protein (FKBP12) with implicit solvation. Binding free energy estimates are in reasonably good agreement with experimental affinities. The conformations of the complexes identified by the simulations are in good agreement with crystallographic data, which was not used to restrain ligand orientations. The BEDAM method is based on λ -hopping Hamiltonian parallel Replica Exchange (HREM) molecular dynamics conformational sampling, the OPLS-AA/AGBNP2 effective potential, and multi-state free energy estimators (MBAR). Achieving converged and accurate results depends on all of these elements of the calculation. Convergence of the binding free energy is tied to the level of convergence of binding energy distributions at critical intermediate states where bound and unbound states are at equilibrium, and where the rate of binding/unbinding conformational transitions is maximal. This finding mirrors similar observations in the context of order/disorder transitions as for example in protein folding. Insights concerning the physical mechanism of ligand binding and unbinding are obtained. Convergence for the largest FK506 ligand is achieved only after imposing strict conformational restraints, which however require accurate prior structural knowledge of the structure of the complex. The analytical AGBNP2 model is found to underestimate the magnitude of the hydrophobic driving force towards binding in these systems characterized by loosely packed protein-ligand binding interfaces. Rescoring of the binding energies using a numerical surface area model corrects this deficiency. This study illustrates the complex interplay between energy models, exploration of conformational space, and free energy estimators needed to obtain robust estimates from binding free energy calculations.     

Enhanced Lipid Diffusion and Mixing in Accelerated Molecular Dynamics

Accelerated molecular dynamics (aMD) is an enhanced sampling technique that expedites conformational space sampling by reducing the barriers separating various low-energy states of a system. Here, we present the first application of the aMD method on lipid membranes. Altogether, ～1.5 μs simulations were performed on three systems: a pure POPC bilayer, a pure DMPC bilayer, and a mixed POPC:DMPC bilayer. Overall, the aMD simulations are found to produce significant speedup in trans-gauche isomerization and lipid lateral diffusion versus those in conventional MD (cMD) simulations. Further comparison of a 70-ns aMD run and a 300-ns cMD run of the mixed POPC:DMPC bilayer shows that the two simulations yield similar lipid mixing behaviors, with aMD generating a 2-3-fold speedup compared to cMD. Our results demonstrate that the aMD method is an efficient approach for the study of bilayer structural and dynamic properties. On the basis of simulations of the three bilayer systems, we also discuss the impact of aMD parameters on various lipid properties, which can be used as a guideline for future aMD simulations of membrane systems.     

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.