γ-环糊精和波尼松龙的相互作用分子动力学模拟和自由能计算

运用分子动力学（Molecular dynamics,MD）和MM—PBSA（molecular mechanics／Poisson Boltzmann surfaeearea）相结合的方法预测了γ-环糊精（γ-cyclodextrin,γ-CD）和波尼松龙的包结模式．在MD模拟过程中,波尼松龙分别采用A环和D环两种取向从γ-CD大口端进入其空腔．在MD轨迹采样基础上,采用高效MM—PBSA方法计算了两种取向的包结自由能．结果表明,计算包结自由能值和实验包结自由能值非常吻合．进一步分析各个能量项,发现范德华相互作用能为包结的主要驱动力．通过比较两种取向的包结自由能大小,预测D环取向为优势包结模式．     

Structure and thermodynamics of alpha-, beta-, and gamma-cyclodextrin dimers. Molecular dynamics studies of the solvent effect and free binding energies

The alpha-, beta-, and gamma-cyclodextrin (CyDs) dimers were studied by molecular dynamics (MD) simulations in water as an explicit solvent. The relative stability of dimers and the involved molecular interactions were determined. Three possible starting orientations were considered for the dimers: head-to-head, head-to-tail, and tail-to-tail. MD simulations were performed over a period of 5 ns to ensure the stability of the system for both the CyD dimers and monomers. The MM-PBSA methodology was used to obtain the free binding energy of the dimers and to determine the most stable arrangement for each solvated CyD. In a vacuum, MD simulations provided the head-to-head orientation as the most stable orientation for the three CyDs, while in aqueous solution the, the head-to-tail orientation was found to be the most stable for the alpha-CyD dimer and the tail-to-tail orientation the most stable for the beta- and gamma-CyD dimers.     

Molecular Dynamics Simulations and Steered Molecular Dynamics Simulations of Glabridin Bound to Wild Type and V30A Mutant Transthyretin: Ligand-linked Perturbation of Tertiary Conformation

Transthyretin(TTR), as a tetrameric protein, functions as a neuroprotector. The native TTR homotetramer dissociates into dimers and monomers. Dimers and monomers self-assemble into amyloid fibrils, and this process can lead to some diseases. Native TTR homotetramer is a widely accepted model for TTR amyloid formation. In this study, simulations using molecular dynamics(MD) and steered MD(SMD) were performed to explore the mechanisms for glabridin(Glab), a specific inhibitor for TTR binding, for V30A mutant and wild-type(WT) TTR. MD simulation results indicate that, compared with Glab binding to WT and V30A mutant, the WT TTR could lead to the collapse of β-strands from Ser52 to His56 at chain A. This phenomenon facilitated the easy dissociation of chains A and C. Calculations of the binding free energy between the two chains showed that the V30A-Glab TTR complex displayed a lower binding energy than other systems(WT TTR and WT-Glab TTR). Then, SMD simulation was performed to explore the unbinding pathway for Glab through the WT and V30A mutant TTR. The results show that Lys15(chain A) produced a hydrogen bond with Glab at the force peak via the WT TTR tunnel. Meanwhile, in the V30A TTR mutant, the hydrogen bond between Lys15(chain A) and Glab was broken at the force peak. This condition was beneficial for Glab to be taken off from the protein. Our theoretical results will be useful in designing a new specific inhibitor of TTR protein to control the TTR homotetramer dissociation.     

The large scale conformational change of the human DPP III-substrate prefers the "closed" form

Human dipeptidyl peptidase III (DPP III) is a two domain metallo-peptidase from the M49 family. The wide interdomain cleft and broad substrate specificity suggest that this enzyme could experience significant conformational change. Long (>100 ns) molecular dynamics (MD) simulations of DPP III revealed large range conformational changes of the protein, suggesting the pre-existing equilibrium model for a substrate binding. The binding free energy calculations revealed tighter binding of the preferred synthetic substrate Arg-Arg-2-naphtylamide to the "closed" than to the "open" DPP III conformation. Our assumption that Asp372 plays a crucial role in the large scale interdomain closure was proved by the MD simulations of the Asp372Ala variant. During the same simulation time, the variant remained more "open" than the wild type protein. Apparently, Ala was not as efficient as Asp in establishing the interdomain interactions. According to the MM-PBSA calculations, the electrostatic component of the free energy of solvation turned out to be higher for the "closed" protein than for its less compact form. However, the gain in entropy due to water released from the interdomain cleft nicely balanced this negative effect.     

A QM/MM study of the binding of RAPTA ligands to cathepsin B

We have carried out quantum mechanical (QM) and QM/MM (combined QM and molecular mechanics) calculations, as well as molecular dynamics (MD) simulations to study the binding of a series of six RAPTA (Ru(II)-arene-1,3,5-triaza-7-phosphatricyclo-[3.3.1.1] decane) complexes with different arene substituents to cathepsin B. The recently developed QM/MM-PBSA approach (QM/MM combined with Poisson–Boltzmann solvent-accessible surface area solvation) has been used to estimate binding affinities. The QM calculations reproduce the antitumour activities of the complexes with a correlation coefficient (r ²) of 0.35–0.86 after a conformational search. The QM/MM-PBSA method gave a better correlation (r ² = 0.59) when the protein was fixed to the crystal structure, but more reasonable ligand structures and absolute binding energies were obtained if the protein was allowed to relax, indicating that the ligands are strained when the protein is kept fixed. In addition, the best correlation (r ² = 0.80) was obtained when only the QM energies were used, which suggests that the MM and continuum solvation energies are not accurate enough to predict the binding of a charged metal complex to a charged protein. Taking into account the protein flexibility by means of MD simulations slightly improves the correlation (r ² = 0.91), but the absolute energies are still too large and the results are sensitive to the details in the calculations, illustrating that it is hard to obtain stable predictions when full flexible protein is included in the calculations.     

Use of MM-PBSA in reproducing the binding free energies to HIV-1 RT of TIBO derivatives and predicting the binding mode to HIV-1 RT of efavirenz by docking and MM-PBSA.

In this work, a new ansatz is presented that combines molecular dynamics simulations with MM-PBSA (Molecular Mechanics Poisson-Boltzmann/surface area) to rank the binding affinities of 12 TIBO-like HIV-1 RT inhibitors. Encouraging results have been obtained not only for the relative binding free energies, but also for the absolute ones, which have a root-mean-square deviation of 1.0 kcal/mol (the maximum error is 1.89 kcal/mol). Since the root-mean-square error is rather small, this approach can be reliably applied in ranking the ligands from the databases for this important target. Encouraged by the results, we decided to apply MM-PBSA combined with molecular docking to determine the binding mode of efavirenz SUSTIVA(TM) another promising HIV-1 RT inhibitor for which no ligand-protein crystal structure had been published at the time of this work. To proceed, we define the following ansatz: Five hundred picosecond molecular dynamics simulations were first performed for the five binding modes suggested by DOCK 4.0, and then MM-PBSA was carried out for the collected snapshots. MM-PBSA successfully identified the correct binding mode, which has a binding free energy about 7 kcal/mol more favorable than the second best mode. Moreover, the calculated binding free energy (-13.2 kcal/mol) is in reasonable agreement with experiment (-11.6 kcal/mol). In addition, this procedure was also quite successful in modeling the complex and the structure of the last snapshot was quite close to that of the measured 2,3 A resolution crystal (structure the root-mean-square deviation of the 54 C(alpha) around the binding site and the inhibitor is 1.1 A). We want to point out that this result was achieved without prior knowledge of the structure of the efavirenz/RT complex. Therefore, molecular docking combined with MD simulations followed by MM-PBSA analysis is an attractive approach for modeling protein complexes a priori.     

Structural Investigations and Binding Mechanisms of Oseltamivir Drug Resistance Conferred by the E119V Mutation in Influenza H7N9 Virus

The use of vaccinations and antiviral medications have gained popularity in the therapeutic management of avian influenza H7N9 virus lately. Antiviral medicines are more popular due to being readily available. The presence of the neuraminidase protein in the avian influenza H7N9 virus and its critical role in the cleavage of sialic acid have made it a target drug in the development of influenza virus drugs. Generally, the neuraminidase proteins have common conserved amino acid residues and any mutation that occurs around or within these conserved residues affects the susceptibility and replicability of the influenza H7N9 virus. Herein, we investigated the interatomic and intermolecular dynamic impacts of the experimentally reported E119V mutation on the oseltamivir resistance of the influenza H7N9 virus. We extensively employed molecular dynamic (MD) simulations and subsequent post-MD analyses to investigate the binding mechanisms of oseltamivir-neuraminidase wildtype and E119V mutant complexes. The results revealed that the oseltamivir-wildtype complex was more thermodynamically stable than the oseltamivir-E119V mutant complex. Oseltamivir exhibited a greater binding affinity for wildtype (−15.46 ± 0.23 kcal/mol) relative to the E119V mutant (−11.72 ± 0.21 kcal/mol). The decrease in binding affinity (−3.74 kcal/mol) was consistent with RMSD, RMSF, SASA, PCA, and hydrogen bonding profiles, confirming that the E119V mutation conferred lower conformational stability and weaker protein–ligand interactions. The findings of this oseltamivir-E119V mutation may further assist in the design of compounds to overcome E119V mutation in the treatment of influenza H7N9 virus patients.     

芬太尼类化合物与阿片μ受体相互作用的分子对接与分子动力学模拟

采用分子对接和分子动力学(MD)模拟方法研究了芬太尼类化合物与阿片μ受体的相互作用机制.先用AutoDock4.0程序将芬太尼类化合物对接到同源模建的阿片μ受体结构中,再用GROMACS程序包在水溶液体系中分别对12个芬太尼激动剂和阿片μ受体蛋白复合物进行了MD模拟研究,优化对接复合物的结构,最后利用MM-PBSA方法,在APBS程序中计算芬太尼类衍生物与阿片μ受体的结合自由能,计算出的受体配合物结合常数(Ki)与其实验值吻合较好,并预测了化合物的活性排序.结果表明,复合物蛋白结构与空载受体蛋白结构有较大差异,特别是胞内区IL2、IL3和跨膜区段TM4骨架构象变化较大,不同的化合物对受体结构影响也有差异,活性较好的化合物会增加蛋白特定区域结构的柔性.芬太尼类化合物可能是通过和受体结合后诱导阿片μ受体构象转变为活性构象,引起一系列的信号传导激活G蛋白,从而引发生理效应.     

Molecular dynamics simulation reveals structural and thermodynamic features of kinase activation by cancer mutations within the epidermal growth factor receptor

The epidermal growth factor receptor (EGFR) is a major target for drugs in treating lung carcinoma as it promotes cell growth and tumor progression. Structural studies have demonstrated that EGFR exists in an equilibrium between catalytically active and inactive forms, and dramatic conformational transitions occur during its activation. It is known that EGFR mutations promote such conformational changes that affect its activation and drug efficacy. The most common point mutation in lung cancer patients is a leucine to arginine substitution at amino acid 834 (L834R). In a recent article, we have studied changes in drug binding affinities due to cancer mutations of EGFR using ensemble molecular dynamics (MD) simulations. Here, we address an enhanced activation mechanism thought to be associated with this mutation. Using extended timescale MD simulations, the structural and energetic properties are studied for both active and inactive conformations of EGFR. The thermodynamic stabilities of these two conformations are characterized by free energy landscapes estimated from molecular mechanics/Poisson-Boltzmann solvent area calculations. Our study reveals that the L834R mutation introduces conformational changes in both states, adjusting the relative stabilities of active and inactive conformations and hence the activation of the EGFR kinase.     

Molecular insight into the interaction between IFABP and PA by using MM-PBSA and alanine scanning methods

The molecular mechanics Poisson-Boltzmann surface area (MM-PBSA) method combined with alanine-scanning mutagenesis is a very important tool for rational drug design. In this study, molecular dynamics (MD) and MM-PBSA were applied to calculate the binding free energy between the rat intestinal fatty acid binding protein (IFABP) and palmitic acid (PA) to gain insight to the interaction details. Equally spaced snapshots along the trajectory were chosen to perform the binding free energy calculation, which yields a result highly consistent with experimental value with a deviation of 0.4 kcal/mol. Computational alanine scanning was performed on the same set of snapshots by mutating the residues in IFABP to alanine and recomputing the DeltaDeltaG(binding). By postprocessing a single trajectory of the wild-type complex, the average unsigned error of our calculated DeltaDeltaG(binding) is below 1.5 kcal/mol for most of the alanine mutations of the noncharged residues (67% in total). To further investigate some particular mutants, three additional dynamical simulations of IFABP Arg126Ala, Arg106Ala, and Arg106Gln mutants were conducted. Recalculated binding free energies are well consistent with the experimental data. Moreover, the ambiguous role of Arg106 caused by the free energy change of the opposite sign when it is mutated to alanine and glutamine respectively is clarified both structurally and energetically. Typically, this can be attributed to the partial electrostatic compensation mainly from Arg56 and the obvious entropy gain in Arg106Ala mutant while not in Arg106Gln mutant. The presented structural model of IFABP-PA complex could be used to guide future studies.     

力场/溶剂水模型对EGFRvⅢ抗原-MR1(scFv)抗体复合物的MM-PBSA计算精度的影响与实验验证

以表皮生长因子Ⅲ型突变体(EGFRvⅢ)抗原多肽与其抗体(MR1)及其人源化突变体的复合物结构为出发点,采用分子动力学中的6种常用力场及3种常用溶剂水模型,分别对上述抗原-抗体复合物进行100ns的分子动力学模拟与分子力学和连续介质模型计算自由能(MM-PBSA),并在实验上利用等温滴定量热(ITC)仪测定了抗原和抗体相互作用的热力学参数.通过在结构变化、能量变化及野生型与突变体比较等几个方面进行综合分析,给出了最佳的计算模型.对不同力场及水模型计算精度等相关问题进行了探讨.     

Assessing the performance of the g_mmpbsa tools to simulate the inhibition of oseltamivir to influenza virus neuraminidase by molecular mechanics Poisson–Boltzmann surface area methods

The molecular mechanics Poisson–Boltzmann surface area (MM-PBSA) method for GROMACS (g_mmpbsa) is an open-source tool that is capable of reading the trajectories generated by GROMACS and calculating the binding free energy using the MM-PBSA method. However, there are multiple force fields available for users to choose from in the GROMACS software, and there are also different solvent water models to combine with the chosen force fields. These different combinations of parameters may significantly impact the results of g_mmpbsa calculation. Unfortunately, the exact combination of force field and solvent water that can well calculate the free energy of the receptor–ligand binding in GROMACS has not been explored yet. To resolve the above issues, this study mainly explored the molecular dynamics (MD) simulations by GROMACS with the six commonly used force fields and three solvent water models, in combination with g_mmpbsa, to calculate the binding free energies of the influenza virus neuraminidase and its mutants with inhibitor oseltamivir carboxylate and compared the present results with previous published results of Amber software from ours and other researchers. Finally, we provided an optimized calculation model, as well as suggestions that may serve as advice and guidance for future computer-aided designs of drug molecules.     

Fluorescence quenching of dyes by tryptophan: interactions at atomic detail from combination of experiment and computer simulation

Fluorescence spectroscopy and molecular dynamics (MD) simulation are combined to characterize the interaction of two organic fluorescent dyes, rhodamine 6G (R6G) and an oxazine derivative (MR121), with the amino acid tryptophan in aqueous solution. Steady-state and time-resolved fluorescence quenching experiments reveal the formation of essentially nonfluorescent ground-state dye/Trp complexes. The MD simulations are used to elucidate the molecular interaction geometries involved. The MD-derived probability distribution of the distance r between the centers of geometry of the dye and quencher ring systems, P(r), extends to higher distances for R6G than for MR121 due to population in the R6G/Trp system of fluorescent interaction geometries between Trp and the phenyl ring and ester group of the dye. The consequence of this is the experimental finding that under the conditions used in the simulations about 25% of the R6G dye is fluorescent in comparison with 10% of the MR121. Combining the above findings allows determination of the "quenching distance", r, above which no quenching occurs. r is found to be very similar (approximately 5.5 A) for both dye/Trp systems, corresponding to close to van der Waals contact. Both experimental dynamic Stern-Volmer analysis and the MD trajectories demonstrate that the main determinant of the fluorescence intensity is static quenching. The approach presented is likely to be useful in the structural interpretation of data obtained from fluorescent conjugates commonly used for monitoring the binding and dynamics of biomolecular systems.     

Origin of mutational effects at the C3 and G8 positions on hammerhead ribozyme catalysis from molecular dynamics simulations

A series of ten 60 ns molecular dynamics (MD) simulations of the native and mutated full length hammerhead ribozymes in the reactant state and in an activated precursor state (G8:2'OH deprotonated) are reported. Mutant simulations include the C3U, G8A, and G8I single mutants and a C3U/G8A double mutant that exhibits an experimental rescue effect. The results provide critical details into the origin of the observed mutation effects and support a mechanism where the 2'OH of G8 acts as a general acid catalyst that is held in position through Watson-Crick hydrogen bonding between G8 and C3.     

残基突变对P53-DNA结合域肽段构象影响的分子动力学模拟

基于分子动力学模拟方法研究了R249S、R248W 和G245S 突变对P53-DNA 结合域肽段(残基230-258)结构的影响. 采用GROMACS 软件包和GROMOS 43A1分子力场, 分别对野生型wtP53肽段、单点突变型P53-R249S肽段、两点突变型P53-R249S/R248W肽段和三点突变型P53-R249S/R248W/G245S肽段进行了4组独立的分子动力学模拟, 每组体系模拟时间为500 ns. 研究结果表明: R249S单残基突变影响肽段残基形成二级结构的情况, 但不改变肽段三维结构的模式, 同时使该肽段结构相对稳定; R249S和R248W两残基同时突变会加剧R249S突变对肽段的影响, 同时导致三维结构发生较大变化, 构象弯曲呈现双turn 结构, 肽段稳定性进一步增大; G245S突变对肽段的影响与R249S和R248W同时突变对其结构的影响相反, 在两残基突变的基础上, G245S突变会使原突变引起的变化减弱甚至消失, 同时使得该肽段结构稳定性减小. 该研究对认识肿瘤致病分子机制和设计新药物有重要意义.     

β-环糊精和甾类化合物的分子动力学模拟

运用分子动力学(molecular dynamics, MD)和MM-PBSA (molecular mechanics/Poisson Boltzmann surface area)相结合的方法预测了β-环糊精(cyclodextrin, CD)和甾类客体分子包结模式. 通过重原子均方根偏差(root mean square deviation, RMSD)分析可得, 两种包结模式下客体分子都可以和β-CD形成稳定的包结. 在MD轨迹采样基础上, 采用高效MM-PBSA方法计算了两种包结模式下的包结自由能. 计算结果显示, β-CD和三个甾类客体分子包结的主要驱动力为范德华相互作用, 而溶剂化能和熵变则不利于体系的包结. 进一步分析平均构象和包结自由能发现, 对于波尼松龙, D-up (D-ring up orientation)取向为优势包结模式; 而乙炔雌二醇和雌三醇的优势包结模式均为A-up (A-ring up orientation)取向. 通过比较β-CD和三个客体分子的理论包结自由能, 预测包结稳定性的次序为乙炔雌二醇＞雌三醇＞波尼松龙, 和实验结果相一致.     

Insights into the Allosteric Effect of SENP1 Q597A Mutation on the Hydrolytic Reaction of SUMO1 via an Integrated Computational Study

Small ubiquitin-related modifier (SUMO)-specific protease 1 (SENP1) is a cysteine protease that catalyzes the cleavage of the C-terminus of SUMO1 for the processing of SUMO precursors and deSUMOylation of target proteins. SENP1 is considered to be a promising target for the treatment of hepatocellular carcinoma (HCC) and prostate cancer. SENP1 Gln597 is located at the unstructured loop connecting the helices α4 to α5. The Q597A mutation of SENP1 allosterically disrupts the hydrolytic reaction of SUMO1 through an unknown mechanism. Here, extensive multiple replicates of microsecond molecular dynamics (MD) simulations, coupled with principal component analysis, dynamic cross-correlation analysis, community network analysis, and binding free energy calculations, were performed to elucidate the detailed mechanism. Our MD simulations showed that the Q597A mutation induced marked dynamic conformational changes in SENP1, especially in the unstructured loop connecting the helices α4 to α5 which the mutation site occupies. Moreover, the Q597A mutation caused conformational changes to catalytic Cys603 and His533 at the active site, which might impair the catalytic activity of SENP1 in processing SUMO1. Moreover, binding free energy calculations revealed that the Q597A mutation had a minor effect on the binding affinity of SUMO1 to SENP1. Together, these results may broaden our understanding of the allosteric modulation of the SENP1−SUMO1 complex.     

Optimizing a biocatalyst for improved NAD(P)H regeneration: directed evolution of phosphite dehydrogenase

Cofactor regeneration is an important solution to the problem of implementing complex cofactor requiring enzymatic reactions at the industrial scale. NAD(P)H-dependent oxidoreductases are highly valuable biocatalysts, but the high cost of the nicotinamide cofactors necessitates in situ cofactor regeneration for preparative applications. Here we report the use of directed evolution to enhance the industrially important properties of phosphite dehydrogenase for NAD(P)H regeneration. A two-tiered sorting method of selection and screening was used in conjunction with random and rational mutagenesis. Following six rounds of directed evolution, soluble expression in E. coli was increased more than 3-fold, while the turnover rate was increased about 2-fold, effectively lowering the cost of the enzyme by >6-fold. Large-scale production of the final mutant enzyme by fermentation resulted in approximately 6-times higher yield (Units/Liter) than the WT enzyme. The enhancements of PTDH were independent of expression vector and E. coli strain utilized. The advantage of the final mutant over the WT enzyme was demonstrated using the industrially relevant bioconversion of trimethylpyruvate to L-tert-leucine. The mutations discovered are discussed in the context of a three dimensional structural model and the resulting changes in kinetics and soluble expression. The engineered phosphite dehydrogenase has great potential for NAD(P)H regeneration in industrial biocatalysis.