gaff力场分子模型

反应力场(reaxff)的分子动力学模拟后处理1.引言1.1 概述反应力场（ReaxFF）是一种用于分子动力学模拟的强大工具，它结合了量子力学的精确性和经典力场的计算效率。

ReaxFF模型能够模拟化学反应的过程，涵盖了分子的构型变化、键断裂和形成以及原子间相互作用的动力学行为。

因此，ReaxFF在材料科学、催化剂设计、生物化学以及能源领域等众多领域中得到了广泛的应用。

本文旨在探索反应力场（ReaxFF）在分子动力学模拟后处理中的重要性和实用性。

通过对分子体系进行动力学模拟并从模拟结果中提取信息，后处理方法能够帮助我们深入理解化学反应的机理和动力学过程。

这对于研究新材料的合成、理论催化研究、药物设计等具有重要意义。

在本文的正文部分，我们将首先介绍反应力场（ReaxFF）模型的基本原理和特点。

通过详细讨论其原子间相互作用的表达式以及参数的调整方法，我们将全面了解ReaxFF模型在化学反应模拟中的优势和局限性。

接着，我们将探讨分子动力学模拟的基本原理，包括算法和计算流程。

了解这些基础知识将有助于我们更好地理解后续章节中的后处理方法。

在结论部分，我们将重点介绍反应力场（ReaxFF）的分子动力学模拟在不同领域的应用情况。

通过实例分析，我们将展示ReaxFF模型对材料性质、催化机制和生物化学反应的模拟结果，并探讨其在相关研究中的意义和贡献。

同时，我们将讨论后处理方法在研究中的重要性，并介绍常用的后处理技术，如分子动力学轨迹分析、聚类算法等。

这些方法可以帮助我们从庞大的模拟数据中提取高质量的信息，从而更好地理解化学反应的本质。

综上所述，本文将以反应力场（ReaxFF）的分子动力学模拟后处理为核心内容进行探讨。

通过深入分析ReaxFF模型的原理和应用，以及后处理方法的意义和方法，我们将对化学反应的机理和动力学行为有更深入的认识。

希望本文能为广大研究者在相关领域的研究提供一些参考和借鉴。

1.2文章结构文章结构是指文章的整体组织框架，它有助于读者理解文章的逻辑结构和内容安排。

力场与拓扑之二：如何选择力场已有9886 次阅读 2015-12-9 13:01 |系统分类:科研笔记2015-12-08 22:27:19待补充, 参考Sobereva力场瞎总结, 磷脂膜模拟的力场瞎总结.性能模拟有机小分子热力学性质用Charmm generalized >= OPLS-AA >= GAFF，但实际上GAFF已经很好了。

它们计算各种有机分子的密度、蒸发焓都很准确，但是介电常数、等温压缩系数计算得都一般。

GAFF对于带有硝基的分子不好。

OPLS和GAFF对于苯甲醛、甲酸，以及有两个及以上Br或Cl相距较近的情况都不好。

对蛋白质构象的模拟: ff99SB+ildn+nmr >CHARMM27 >OPLS »f99。

•Berger：专门用于磷脂的力场•PFF=Polarizable Force Field•VAMM=Virtual atom molecular mechanics材料力场•cvff(consistent valence forcefield)：参数用于有机分子、蛋白质模拟，函数形式略复杂。

cvff_aug是对其扩展，可以用于研究硅酸盐、铝硅酸盐、磷酸盐、泥土•CFF(consistent family of forcefield)：包括CFF91和CFF95。

适用面很广，涵盖有机无机小分子、聚合物、多糖和生物大分子，还支持金属。

函数形式挺复杂。

参数由从头算获得，非键参数从CVFF弄来，不适合凝聚相模拟。

•pcff：基于CFF91，适用范围做了扩展，主要用于聚合物和有机材料，也能用于无机材料，还有糖、核酸、脂的参数。

•COMPASS=Condensed-phase Optimized Molecular Potentials for Atomistic Simulation Studies：在pcff基础上改进的新版本，同样由从头算获得参数，在凝聚相模拟方面大有改善。

环糊精和雌二醇的分子动力学模拟和MM PBSA 自由能计算姚雪霞(南京农业大学工学院,江苏南京210031)摘要 运用分子动力学和MM PBS A 相结合的方法预测 环糊精和雌二醇包结模式。

通过MD 轨迹分析,不论是A up 取向还是D up 取向,雌二醇都可以和 CD 形成稳定的包结;通过氢键占有率的分析,发现 CD 的构象发生一定程度的变化。

MM PBS A 计算结果表明,A up 取向包结自由能更低,为优势包结模式。

进一步分析各个能量项可得,范德华相互作用能为包结的主要驱动力。

关键词 MM PBS A;分子动力学; -环糊精;雌二醇中图分类号 O 641 文献标识码 A 文章编号 0517-6611(2008)29-12557-03M olecu l ar Dy na m ics Si m ul ati on and MM PBS A Free Energy Cal cul ati on of cyc l odextr i n a nd Estrad i ol YAO Xue xi a (Co ll ege o fEng i neeri ng ,N anji ng Ag ri cult ura lUn i versit y ,N an ji ng ,Ji angsu 210031)Abstract The possi ble bi ndi ng m ode bet ween cyc l odex tri n ( C D )and estradi o lwas pred i cted by using mo lecul ar dynam i cs (M D )and MM PBSA (m olecul ar mechanics/Po i ss on Boltz mann surf ace area)method .Bot h A up orient a ti on and D up orient a tion of estradi o l coul d for m stabl e complexesw ith CD based onMD trajectory anal y sis .The conf or m ati on o f CD had so m e extent change by i nvesti gati ng t he occupa ti on of hydrogen bonds of CD.T he co mputed results by MM PBS A techniq ue showed t hat t he bi nd f ree energy ofA up or i entati onw as l ower than that o fD up or i entati on ,theref ore A up o ri entati on belonged t o preferenti a lb i ndi ng mode .By further analyzing each energy it em ,it coul d be obtai ned t hat van derW aals'i nteracti on energy w as t he ma i n dri v i ng f orce f or compl ex .K ey words MM PBS A;M o l ecular dyna m i cs ; cycl odex tri n ;Estrad i ol基金项目 南京农业大学青年科技创新基金(K J 07033)。

gaff力场分子模型（原创实用版）目录1.GAFF 力场分子模型简介2.GAFF 力场的发展历程3.GAFF 力场的特点和优势4.GAFF 力场在我国的应用现状与前景正文一、GAFF 力场分子模型简介GAFF（General Amino Acid Force Field）力场分子模型是一种用于描述生物大分子（如蛋白质和核酸）中分子间相互作用力的物理模型。

在生物物理学、生物信息学和药物设计等领域具有广泛应用。

该模型通过大量实验数据和理论计算参数化，以模拟生物大分子在各种条件下的结构和动力学行为。

二、GAFF 力场的发展历程GAFF 力场分子模型的发展始于 20 世纪 90 年代，由美国加州大学洛杉矶分校的研究团队提出。

该团队基于早期的力场模型，如 MM3、MM4 和MM5 等，结合大量实验数据，对力场参数进行了优化和改进，最终形成了GAFF 力场。

此后，GAFF 力场经过多次更新和完善，现已成为生物大分子研究领域的重要工具。

三、GAFF 力场的特点和优势1.高精度：GAFF 力场分子模型通过大量的实验数据和理论计算参数化，具有较高的模拟精度，能够较好地反映生物大分子在实际环境中的结构和动力学特征。

2.广泛适用性：GAFF 力场适用于多种生物大分子，包括蛋白质、核酸、多糖等，为研究不同类型生物大分子的相互作用提供了便利。

3.灵活性：GAFF 力场提供了多种算法和参数设置，可以根据研究者的需求和计算资源进行选择，适应不同应用场景的需求。

4.开源性：GAFF 力场分子模型的源代码和相关工具对外公开，便于研究者进行二次开发和定制，促进了该领域的技术创新和应用拓展。

四、GAFF 力场在我国的应用现状与前景近年来，我国在生物物理学、生物信息学和药物设计等领域取得了显著进展，GAFF 力场分子模型在这些成果中发挥了重要作用。

目前，我国已有多个研究团队在 GAFF 力场的基础上，开展了针对特定生物大分子或疾病的研究，取得了一系列原创性成果。

CHEMICAL INDUSTRY AND ENGINEERING PROGRESS 2018年第37卷第2期·702·化 工 进展Ponatinib 抑制野生型和T315I 突变型BCR-ABL1激酶的分子模拟陈清，王坚毅（广西大学化学化工学院，广西 南宁 530004）摘要：药物普纳替尼（Ponatinib ）对野生型和T315I 突变型BCR-ABL1激酶均有较强抑制作用，然而其抑制机理仍然未知。

为了更好地了解Ponatinib 与BCR-ABL1激酶的作用机理，本文首先运用分子动力学模拟方法研究Ponatinib 对BCR-ABL1激酶构象变化的影响，然后采用MM-PBSA 方法计算复合物的结合自由能。

结果表明：Ponatinib 诱导野生型BCR-ABL1激酶P-loop 和铰链区相互靠近，导致野生型激酶结合口袋关闭，有利于激酶与药物结合；而Ponatinib 诱导T315I 突变型BCR-ABL1激酶P-loop 和铰链区相互远离，导致T315I 突变型激酶结合口袋打开，对激酶与药物结合不利；Ponatinib 与野生型以及T315I 突变型BCR-ABL1激酶的结合自由能分别为–58.57 kcal/mol 和–43.54kcal/mol ，Ponatinib 对野生型BCR-ABL1激酶的抑制能力明显强于T315I 突变体，与文献报道的实验抑制活性结果一致。

研究结果对认识靶点蛋白和抑制剂分子机制以及设计新药有重要意义。

关键词：分子模拟；酶；T315I 突变；普纳替尼；动力学理论；活性中图分类号：Q754 文献标志码：A 文章编号：1000–6613（2018）02–0702–06 DOI ：10.16085/j.issn.1000-6613.2017-0891Inhibition of Ponatinib on wild-type and T315I mutant BCR-ABL1 kinasebased on molecular simulationCHEN Qing ，WANG Jianyi（School of Chemistry and Chemical Engineering ，Guangxi University ，Nanning 530004，Guangxi ，China ）Abstract ：The drug Ponatinib had a strong inhibitory effect on wild-type and T315I mutant BCR-ABL1 kinase ，but its inhibition mechanism was still unknown. To have a better understanding of the interaction mechanism between Ponatinib and BCR-ABL1 kinase ，the effect of Ponatinib on the conformational changes of BCR-ABL1 kinase was studied by molecular dynamics simulation at first ，and then the binding free energy of complex was calculated by the MM-PBSA method. The results indicated that the P-loop and hinge region of wild-type BCR-ABL1 kinase approached each other induced by Ponatinib ，leading to the closure of the binding pocket of wild-type BCR-ABL1 kinase ，which would be favorable to the binding of drug. In contrast ，the P-loop and hinge region of T315I mutant BCR-ABL1 kinase move far away from each other induced by Ponatinib ，resulting in the opening of the binding pocket of T315I mutant BCR-ABL1 kinase ，which was no advantage for kinase binding with drug. The binding free energy of Ponatinib bound to wild-type and T315I mutant BCR-ABL1 kinase were –58.57kcal/mol and –43.54kcal/mol ，respectively. The results showed that the inhibitory ability of Ponatinib on wild-type BCR-ABL1 kinase was significantly stronger than that of T315I mutant ，which was consistent with the experimental inhibition activity reported in the literature. 。

㊀㊀2023年4月第38卷第2期㊀㊀㊀㊀㊀㊀㊀㊀JOURNAL OF LIGHT INDUSTRY㊀Vol.38No.2Apr.2023㊀收稿日期:2022-05-23基金项目:国家自然科学基金青年基金项目(31901676,32001352);广州市基础与应用基础研究项目(202102080221);贵州省科技计划项目[黔科合支撑(2020)1Y171号];阳西县级科技项目(21011)作者简介:傅亮(1968 ),男,湖南省益阳市人,暨南大学副教授,博士,主要研究方向为食品科学㊂E-mail :tfuliang@ 通信作者:陈永生(1985),男,河南省漯河市人,暨南大学副研究员,博士,主要研究方向为食品化学㊁食品营养学㊂E-mail :chysh11@傅亮,吕金羚,张锦,等.分子模拟技术在食品组分互作体系及安全领域的应用研究进展[J].轻工学报,2023,38(2):1-13.FU L,LYU J L,ZHANG J,et al.Research progress on the application of molecular simulation technology in food component interaction system and safety[J].Journal of Light Industry,2023,38(2):1-13.DOI:10.12187/2023.02.001分子模拟技术在食品组分互作体系及安全领域的应用研究进展傅亮1,吕金羚1,张锦1,庄国栋1,3,朱勇2,陈永生11.暨南大学理工学院,广东广州510632;2.贵州大学酿酒与食品工程学院,贵州贵阳550025;3.广东药科大学中药学院,广东广州511400摘要:在食品加工和储存过程中,组分分子之间会发生复杂的相互作用,传统的实验方法难以从分子层面对影响食品营养价值和安全性的互作机制进行清晰直观的阐释㊂分子模拟技术基于分子模型研究分子行为以预测或解释宏观实验现象,具有客观㊁高效㊁成本低廉等优点,是连接微观与宏观尺度的重要桥梁㊂对该技术在食品组分互作体系㊁食品安全领域中的应用研究进展进行了综述,认为:分子模拟技术能够从原子/分子层面揭示食品营养物质分子间的互作机制㊁食品污染物的毒理机制及抑菌机理,表征分子构象的动态变化㊁电子转移㊁共价键断裂和重建等可视化数据,为改善食品的加工工艺和储存条件㊁提高食品的安全性等生产实践提供指导性建议㊂然而,近年来该技术在应用过程中存在大多数互作体系的研究相对独立㊁高反应性的分子在模拟过程中被认为是纯分离物或惰性物质等问题,今后应深入研究并逐步完善环境诱导下食品组分分子的构象变化㊁功能特性与食品品质三者之间的联系,以期为分子模拟技术在食品不同领域的实际应用提供参考㊂关键词:分子模拟技术;食品组分;互作机制;食品安全中图分类号:TS201.2㊀㊀文献标识码:A㊀㊀文章编号:2096-1553(2023)02-0001-130　前言食品能够给人体提供丰富的营养物质,在食品加工和储存过程中,营养物质之间及其与环境之间会发生复杂的变化,影响食品的感官品质㊁营养价值和安全性[1]㊂目前,主要通过一些宏观信息对分子间的相互作用进行研究,但使用的某些精密仪器或多或少存在对样品纯度要求高㊁实验耗时久㊁分离能力差㊁灵敏度低等局限性,因此需要结合多种研究方法从不同角度进行研究分析㊂有研究[2]表明,食品中营养物质的结构与功能特性关系密切,因此研究人员亟需从微观层面获取更多的信息以提高对研究体系的认识水平㊂分子模拟(Molecular Simulation,MS)技术是一㊃1㊃㊀2023年4月第38卷第2期㊀种利用计算机强大的数据处理和图像显示能力,基于原子水平的分子模型,模拟分析分子的结构与行为㊁体系的各种理化性质并指导进一步实验的技术[3]㊂与传统实验技术相比,分子模拟技术具有可预测㊁高效率㊁高安全性㊁低成本等优点,可以弥补传统宏观实验方法的不足,受到业界极大的关注[4],广泛应用于微观层面的食品中营养物质分子结构动态变化㊁食品组分分子间互作机制㊁食品污染毒性机理等研究中[5-6]㊂本文拟综述分子模拟技术在食品组分互作体系中的应用,以及解决组分互作造成的食品安全问题的研究进展,并讨论该技术在应用过程中存在的问题和未来的发展方向,以期为该技术在食品领域的广泛应用提供参考㊂1㊀分子模拟技术的发展自第一次模拟实现刚性球之间弹性碰撞以来,分子模拟技术得到了很大的发展㊂20世纪60年代,研究者开始利用分子模拟技术逐渐深入了解简单和复杂流体的结构和动力学㊂相关的技术程序不断推陈出新,有大分子模拟程序CHARMM,也有在此基础上开发的AMBER程序和GROMOS程序㊂20世纪90年代后,分子模拟技术快速发展[7],目前,该技术已经可以模拟小分子及蛋白质㊁多糖等多种复杂体系的相互作用,模拟程序GROMACS㊁AutoDock㊁SYBYL等也陆续问世㊂特别是1998年和2013年诺贝尔化学奖的授予,分别肯定了科学家在量子化学计算方法与 发展多尺度模型研究复杂化学体系 领域的卓越贡献[8],这表明分子模拟技术已经与试管实验同样重要㊂分子模拟技术逐渐成为食品领域研究的重要工具,在ScienceDirect上分别以 molecularsimulation(分子模拟) 及 molecular-simulation+food(分子模拟+食品) 为关键词进行检索,其出版物数量在2000 2021年间呈逐年飙升的趋势㊂由此可见,该技术在食品领域的前景十分广阔,其应用研究代表了目前食品研究的部分热点发展方向㊂2㊀应用于食品领域的分子模拟技术主要方法㊀㊀食品领域常用的分子模拟技术方法主要包括量子力学法㊁蒙特卡罗法㊁分子对接和分子动力学法[9]㊂2.1㊀量子力学法量子力学法是基于薛定谔方程近似解研究物质世界微观粒子的结构㊁性质及其运动规律的模拟计算方法,主要分为从头算法㊁半经验方法及密度泛函理论(Density Functional Theory,DFT)[10]㊂其中,从头算法基于量子力学定律和光速㊁普朗克常数等物理常数的数值,计算过程中不考虑电子相关能㊁高估电子排斥㊁耗时久;半经验方法基于实验数据的参数优化求解过程,能对分子系统进行合理的定性描述和准确的定量预测,计算效率高㊁成本低;DFT基于仅有3个变量的电子密度,计算速度快㊁精度高,不仅考虑了电子相关能,还能探究电子间的交联和电子对的瞬时互作㊂化学反应即共价键的生成和断裂,本质上是原子表层电子的分布变化,理论上只能用量子力学法阐释㊂相比之下,分子力学法的本质是假设原子为刚性球体,通过解牛顿力学方程来描述体系变化,不考虑电子转移,不适合描述化学反应中共价键的变化㊂目前,量子力学法广泛应用于食品活性成分的作用机理和构-效关系的研究中[11],其中蛋白方面的拓扑文件常用AMBER等力场,小分子拓扑的构建则常用Gaussian㊁ORCA等程序运算生成㊂此外,在食品酶的催化机制㊁食品营养因子㊁风味成分的化学反应等涉及共价键变化的研究中,也常使用基于DFT的量子力学法或将量子力学法与分子动力学联用进行分析[11-14]㊂2.2㊀蒙特卡罗法蒙特卡罗法是一类把随机变量作为研究对象,求解数理及工程技术问题近似解的数值模拟法,常用于食品安全风险评估[15]㊁测量不规则形状食品的体积[16]㊁分析食品中细菌群体行为的变异性[17]等㊂蒙特卡罗法适用于模拟常见的关键事件,且模拟结果与实际情况随着抽样模拟次数的增加逐渐逼近,计算精度也随之提高㊂试验风险分析能综合反映多个因素的共同影响,不受因素间关系的影响,避免主观因素造成的风险分析误差,从而达到减化繁琐操作或避免危险工作的目的㊂但如果关键事件很少发生,使用蒙特卡罗法评估事件发生的概率时可能存在大量不切实际的迭代,导致模拟过程耗时长㊁效㊃2㊃㊀傅亮,等:分子模拟技术在食品组分互作体系及安全领域的应用研究进展率低㊂2.3㊀分子对接分子对接是研究分子(受体和配体)之间非键相互作用,预测结合模式(作用位点和结合构象)和亲和力的模拟方法[5-18],可为揭示分子间的作用机制及设计新的配体提供微观层面的理论辅助和决策基础㊂分子对接按系统简化的程度和方式可分为刚性对接㊁柔性对接和半柔性对接㊂刚性对接过程中配体结构与受体结合位点的构象稳定不变,计算速度快,但模拟精度不足;柔性对接虽允许研究体系的构象自由改变,但变量随体系原子数的增多呈几何级数增加,计算量和时间成本消耗较大,对模拟平台要求较高;半柔性对接在对接精度和计算量上均采用折中策略,控制配体构象变化,使受体构象不变或只有部分残基变化,广泛应用于研究大分子-小分子的识别和互作机制㊂分子对接在食品领域常用于研究蛋白质-配体的结合模式和亲和力及食品污染物的毒理机制等,目前常用的对接软件有AutoDock㊁DOCK㊁Gold 等[19]㊂2.4㊀分子动力学法分子动力学法[4]是分子模拟技术中最接近真实条件的模拟方法,基于经典牛顿力学中的分子力场模拟计算分子体系动态行为㊁热力学及动力学性质[20],常与分子对接法联用以描述分子构象的动态变化过程,进而获得分子间的结合模式㊁结合亲和力㊁瞬态稳定性等互作机制㊂相比之下,分子力学法只能基于经验或半经验力场描述体系的静态性质㊂㊀㊀㊀分子动力学法在模拟计算时不考虑分子力场的参数和电子运动,只考虑化学键的旋转伸缩和键角的改变,并选取在合适力场下势能最低的稳定构型进一步分析㊂因此,分子力场的选择对于分子动力学是否能获得最接近真实状态的模拟计算结果至关重要㊂分子动力学法模拟计算速度快㊁精度高,不仅可以处理具有数万粒子的系统,还能描述不同温度下分子体系热力学状态的变化㊂但该方法模拟的时间和空间尺度有限,为了获取更多的动态信息,一般通过分子力学-蒙特卡罗法联用来优化体系结构以缩短模拟时长㊂目前,分子动力学法常用的分子力场是AMBER㊁CHARMM㊁OPLS 等[1,4,21](见表1),常用的模拟程序有GROMACS㊁AMBER㊁CHARMM 等[21-22](见表2)㊂3㊀分子模拟技术在食品组分互作体系中的应用3.1㊀脂质及其衍生物3.1.1㊀理化反应的分子机理及其产物分析㊀油脂中含有大量的甘油三酯,易受加工工艺和储存环境的影响而发生水解㊁热解㊁氧化等复杂变化,产生烷烃等系列终产物,影响油脂的营养价值和功能特性[4]㊂分子模拟技术能监测复杂脂质体系中分子结构的动态变化,阐明变化机理和产物类别㊂G.Bella 等[14]基于DFT 计算研究以甘油三酯为主的植物油混合物体系中主要成分的化学位移,为核磁共振等传统光谱的区域拥挤现象提供了合理的可视化㊀㊀表1㊀分子动力学模拟中常用的力场Table 1㊀The frequently-used force field in molecular dynamics simulation力场获取方法应用范围兼容程序网站许可情况AMBER Antechamber,acpype生物大分子㊁一些有机小分子GROMACS,CHARMM https:ʊ /AmberModels.php免费GROMOS ATB,PRODRG2蛋白质㊁核酸㊁烷烃的凝聚相GROMACShttp:ʊvienna-ptm.univie.ac.at /?page_id =100免费OPLS LigParGen,MKTOP,TPPmktop 有机小分子㊁蛋白质的凝聚相GROMACS,CHARMM,NAMD,LAMMPShttp:ʊ /oplsaam.html免费CHARMM CGenFF,SwissParam蛋白质㊁脂类(膜)㊁核酸㊁小分子药物GROMACS,CHARMMhttp:ʊ /charmm_ff.shtml#gromacs免费ReaxFF ADF,Materials studio,LAMMPS高温㊁燃烧㊁爆炸等极端过程中的化学变化过程ADF,Materials studio,LAMMPS https:ʊ /doc /ReaxFF /Included_Forcefields.html收费GAFF Antechamber,acpype 有机小分子AMBER,GROMACS,CHARMMhttps:ʊ /AmberModels.php 免费MARTINImartinize.py磷脂㊁胆固醇㊁蛋白质㊁碳水化合物㊁聚合物GROMACShttp:ʊcgmartini.nl /index.php 免费㊃3㊃㊀2023年4月第38卷第2期㊀表2㊀分子动力学模拟常用的程序Table 2㊀The frequently-used programs in molecular dynamics simulation程序应用范围优点兼容力场网站许可情况GROMACS 磷脂双层膜㊁蛋白质㊁药物分子等复杂的生物系统;聚合物㊁无机物等非生物体系能处理数百万原子的复杂系统,效率高GROMOS,OPLS,AMBER,CHARMM http:ʊ /免费AMBER 蛋白质㊁碳水化合物㊁核酸等生物大分子支持建模和蒙特卡罗法AMBER,GAFFhttp:ʊ /工具和力场免费;模拟程序软件包收费CHARMM 生物小分子和大分子基于经验能量计算出体系的相互作用能和构想能㊁局域最小化㊁动力学行为等CHARMM,AMBER,GAFF,OPLShttp:ʊ /力场免费;模拟程序软件包收费LAMMPS 软㊁固态材料能模拟金属系统OPLS,ReaxFFhttps:ʊ /免费NAMD 蛋白质㊁核酸㊁细胞膜等体系全原子模拟,脚本语言简单㊁并行能力良好,能介观或跨尺度模拟,模拟速度快CHARMM,x-plor,AMBER,Dreiding http:ʊ /Research /namd /收费Discovery Studio 生命大分子的建模和模拟,药物设计功能模块化且高度集成的综合模拟平台CHARMm,CFF,MMFFhttp:ʊwww.3dsbiovia.com /products /collaborative-science /biovia-discovery-studio /收费Materials Studio材料科学领域模拟聚合㊁催化㊁表面重构等多种化学反应过程将材料模拟引入桌面,能同时考虑周期性和非周期性溶剂化效应COMPASS,PCFF,CVFF,ReaxFFhttps:ʊwww.3dsbiovia.com /products /collaborative-science /biovia-materials-studio /收费预测㊂但由于力场的兼容性有限,仅用一种力场无法详细描述脂质体系涉及的所有理化反应机理㊂Z.Q.Zhang 等[23]采用基于ReaxFF 力场的分子动力学模拟研究植物油中不饱和甘油三酯的主要热解机制,发现初始分解过程中C O 键的断裂和脱羧发生在C C 键β位的C C 键裂解之前,最终热解产物主要包括烷烃㊁烯烃㊁CO 2㊁H 2等㊂虽然基于ReaxFF 力场的分子动力学模拟能在原子水平准确地描述化学反应过程中的电子转移㊁化学键断裂和重建情况,与实验结果吻合度高,但其只能在燃烧或高温等极端条件下模拟仅包含C㊁H㊁O 等元素的特定体系化学变化,缺少对酯化㊁分解等脂质体系经常发生的化学反应的描述㊂由此可见,一些力场在脂质体系中的应用还有很大的拓展空间㊂3.1.2㊀乳液体系的动态监测及其稳定的分子机制分析㊀乳液是水分散在油中(W /O)或油分散在水中(O /W)的双向分散体㊂牛奶㊁奶油㊁冰淇淋等许多食品中均存在乳液或在加工过程中经历乳化,乳液体系在聚集方面的稳定性决定了食品的品质㊂分子模拟技术能动态监测乳液的聚集情况㊁预测体系不同组分对乳液稳定性的影响㊂M.Koroleva 等[24]基于Langevin 动力学开发了研究W /O 乳液中液滴聚集的仿真程序,可以模拟监测不同液滴数量的聚集体随时间的变化情况,为模拟乳液液滴的运动提供了更真实的描述㊂乳液也常作为生物活性小分子的递送载体以提高小分子的生物利用率㊂分子模拟技术能在分子水平上揭示功能性乳液体系中分子间的稳定机制,可作为传统实验结果的补充和验证㊂以卵磷脂为油相制备的姜黄素负载乳液稳定性高且溶解度好,分子动力学模拟表明这是姜黄素通过氢键和水桥与脱脂淀粉形成稳定络合物体系造成的结果[25]㊂3.1.3㊀脂质衍生物和其他食品组分的互作机制分析㊀鸡蛋黄中包含中性脂肪㊁卵磷脂和胆固醇,常用于面包㊁馒头等烘焙食品中㊂S.Y.Sang 等[26]研究发现,蛋黄卵磷脂能促进淀粉糊化,降低凝胶强度,差示扫描量热法得出这是双螺旋交联间距增大造成的结果㊂但由于差示扫描量热法计算淀粉糊化度的经验公式不固定,且仅能通过淀粉糊化全过程中糊化温度和吸热的变化等宏观数据判断,因此又通过分子动力学模拟进一步发现卵磷脂与直链淀粉在交联过程中的单螺旋结构,深化了对卵磷脂促进糊化㊃4㊃㊀傅亮,等:分子模拟技术在食品组分互作体系及安全领域的应用研究进展的分子机制的认识㊂生物表面活性剂比合成表面活性剂更安全㊁环保㊂分子模拟技术能通过构建分子模型深入研究脂质衍生物类表面活性剂的构-效关系,揭示表面活性剂诱导蛋白质构象变化行为的分子机制㊂C.Russell等[27]研究发现,鼠李糖脂比甘油单酯或卵磷脂对卵清蛋白的结合亲和力更大,具有从油-水界面取代蛋白质的能力,是3种表面活性剂中唯一一种能够以高物质的量比为乳液提供热稳定性的表面活性剂㊂综上所述,分子模拟技术已广泛应用于阐明脂质及其衍生物体系中存在的理化反应机理和分子互作机制㊂随着对分子模拟技术的深入了解,研究人员能通过构建粗粒度模型延长复杂脂质系统的模拟尺度㊂但脂质及其衍生物易发生氧化㊁热解等反应并生成大量结构未知的中间产物和反应产物,这些产物也会与其他潜在食品组分互作进而改变体系的结构和功能特性㊂因此,在常见的食品加工和储存条件下,脂质及其衍生物与蛋白质㊁碳水化合物等其他食品组分的微观结构变化和功能活性之间的联系仍有待进一步研究,热加工-脂质-蛋白质㊁脂质-蛋白质-碳水化合物等复杂三元体系模型的构建也有助于深入探究脂质体系变化程度对食品品质的影响㊂3.2㊀蛋白质食品蛋白质的营养功能特性与其结构关系密切㊂但蛋白质结构复杂,相对分子质量大的蛋白质复合物不易结晶㊁难以分离,因此X射线晶体衍射技术或冷冻电镜技术不适用于解析其晶体结构㊂核磁共振法虽能在近生理状态下获得蛋白质的动态结构信息,但对样品浓度要求高,无法分析相对分子质量大㊁不溶性的蛋白结构,且图谱解析耗时久[28]㊂分子模拟技术能弥补这些传统方法的不足,从分子层面深入探讨周围环境对蛋白质结构和功能的影响机制㊁蛋白质-配体的互作机制,并提供蛋白质的构象动态信息㊂3.2.1㊀酶的作用机理分析㊀分子模拟技术广泛应用于动态分析酶促反应的全过程,设计新的酶,分析酶的底物㊁产物㊁抑制剂和激活剂,研究酶的固定化和失活[29-34]㊂壳聚糖酶OU01在催化水解过程中经历了开放-封闭-开放的构象转变[29],主要应用于工业目的新酶的设计㊂分子模拟技术可与传统光谱实验相结合,分析不同特异性酶水解降低β-乳球蛋白(β-LG)抗原性机制的差异[30],找到能显著降低β-LG抗原性的地衣芽孢杆菌蛋白酶(结合能-5.66kcal/mol),指导乳制品产业优化加工方法,并促进低致敏性乳制品的发展㊂然而,分子对接提供的溶剂条件有限,模拟过程中一般不考虑蛋白质的构象变化,也无法完全模拟如溶剂化作用等的真实情况㊂F.Xie等[33]在研究木质素的α-淀粉酶激活作用时发现,分子对接模拟所获得的结合能(-5.84kcal/mol)与298K条件下热力学分析获得的吉布斯自由能变ΔG(-14.11kcal/mol)有偏差㊂3.2.2㊀蛋白质的界面性质分析㊀蛋白质是两亲性分子,可吸附在界面上,通过降低界面张力增加食品的稳定性㊂有研究[6]表明,蛋白质吸附行为可导致其构象和功能性质发生变化㊂因此,深入研究蛋白质的界面性质对于控制加工和储存过程中富含蛋白质食品的品质是十分必要的㊂目前研究界面性质的方法主要是张力法和流变测量,但这些方法只能通过界面凝胶或晶体的形成㊁X射线检测吸附层的厚度等宏观信息来间接推断蛋白质构象的改变㊂分子模拟技术通过建立原子模型能够在分子层面模拟蛋白质在界面上的构象动态行为,为研究蛋白质的结构-功能-界面行为三者之间的联系提供详细的直观信息㊂H.Schestkowa等[35]通过分子动力学模拟发现,油-水界面蛋白质分子暴露的疏水域对吸附速率的影响大于静电排斥,这也印证了悬滴实验的结果㊂D.L.Cheung[2]通过分子动力学模拟发现,肌红蛋白衍生肽1-55在空气-水界面上会变换3种构象(1种紧凑㊁2种延伸),肽56-131在界面上形成致密的天然构象(见图1),肽1-55是很好的乳化剂和发泡剂,而肽56-131与之完全相反㊂由此可见,蛋白质的吸附构象与功能性质之间存在密切联系,分子模拟技术能直观再现蛋白质界面吸附过程中的构象变化,与实验结果相结合则可进一步阐明蛋白质的界面吸附行为与其功能特性的联系,为改善富含蛋白质食品的营养功能特性及加工和储存过程中的稳定性提供有效信息㊂这些研究也表明,蛋白质晶体和界面环境复杂多变,而目前分子模㊃5㊃㊀2023年4月第38卷第2期㊀㊀㊀㊀图1㊀肌红蛋白多肽在空气-水界面吸附过程中的构象变化[2]Fig.1㊀Conformational changes of myoglobin peptides during adsorption at the air-water interface [2]拟技术的模拟尺度有限,蛋白质㊁界面等分子模型的精确性还有很大的提升空间㊂3.2.3㊀食源性活性肽的感官品质分析㊀食源性活性肽具有降压㊁降糖等特性,但在功能性食品中可能会因味道的限制而影响其作为营养素的使用㊂分子模拟技术能建立食源性活性肽的结构㊁功能特性和感官品质之间的内在联系㊂P.Zhou 等[36]通过定量构效关系(QSAR)的非线性建模揭示了短肽的血管紧张素Ⅰ转换酶(ACE)抑制与苦味之间呈显著正相关,但三肽和四肽的相关性小㊂在此基础上,利用量子力学-分子动力学计算分析发现,四肽能与ACE 有效结合,氨基酸残基继续增多也不会实质性地增强肽对ACE 的亲和力,表明三肽和四肽可以作为具有抗高血压活性和低苦味功能性食品的潜在配料㊂传统的Q 值法仅考虑肽的平均疏水性和苦味之间的关系,相较而言,QSAR 能全面分析氨基酸残基的空间性质,包括疏水性㊁空间位置㊁氨基酸组成㊁净电荷㊁侧链分子大小㊁α-螺旋的倾向性等对苦味的影响,所建立的预测模型对大多数肽具有高拟合能力和可靠的稳定性,对建立食源性活性肽的构-效关系和感官品质之间的联系有显著作用㊂3.2.4㊀蛋白质与其他食品组分的互作机制分析㊀蛋白质不能单独作用,其与配体分子的相互作用和特异性识别是蛋白质发挥其生物功能的重要途径之一,也一直是食品科学领域研究的前沿和热点[37]㊂防冻蛋白质能减少解冻损失,最大程度地保护冷冻食品的品质㊂传统的NMR 和CD 光谱技术很难表征复杂且灵活的晶体结构,导致防冻蛋白质的防冻机制相关研究不足,严重阻碍了其在食品工业领域的应用㊂J.H.Wu 等[38]通过分子动力学模拟构建了分子结构模型,发现新型丝胶蛋白肽可以通过氢键和疏水相互作用吸附在冰表面,预防冰晶生长㊂防冻蛋白肽防冻性能的关键是碳水化合物的羟基及其与周围水分子形成水桥的能力,这些水合能力通过阻止冰晶再生引起的机械性细胞损伤来保持冷冻食品的品质[39]㊂除此之外,分子模拟技术也常用于表征蛋白-蛋白互作过程中的构象变化机制㊂K.K.Darmawan 等[40]用分子动力学模拟发现,尽管蛋白质之间在胃液pH 值条件下会相互排斥,但β-乳球蛋白与载脂蛋白-乳铁蛋白(apo-Lf)之间的主要相互作用力依然存在,可能形成能够截留额外生物活性营养物质的纳米级絮凝剂样复合物㊂此外,持久存在的相互作用力也使apo-Lf 抗菌肽区域结构刚性增强,表明其在胃消化过程中对有害细菌具有潜在作用能力㊂由此可见,分子模拟技术能为深入理解蛋白质的初始变性过程和潜在功能特性提供直观的分子构象信息㊂综上所述,分子模拟技术已广泛应用于在分子水平上研究酶与底物结合的活性位点㊁蛋白质的界面吸附行为和防冻机制,有助于研究人员深入研究蛋白质与其他食品组分之间的互作机制,理解蛋白质结构与功能之间的关系,优化食品的加工和储存环境,指导开发更有益于人体健康的功能性食品㊂但在真实环境下,大部分蛋白质的结构仅仅维持在稳定边缘,并不具备完美的物理学结构㊂氨基酸序列的组合㊁无序区域及主侧链空间构象的变化使蛋白质的结构非常复杂,所构建模型中存在的细小偏差都会影响模拟过程中分子间的识别契合及相互作㊃6㊃。

GAFF立场在Gromacs中的实现GAFF立场在GMX软件中的实现步骤1 小分子的GAFF参数可以用Antechamber来构建，电荷目前有2个比较合理的处理小分子，AM1-bcc和RESP电荷。

AM1-bcc可以用antechamber直接得到，也可以用 chimera来实现。

具体可以结合自己会的来做即可。

Resp需借助第三方程序来处理得到（高斯+amebr）。

拿到小分子GAFF参数后，需要用脚本来进行转化成GMX识别的参数文件，因为GMX和amber 格式的参数文件很不一样，这样就导致程序识别不了GAFF的参数文件。

转化可以用amb2gmx.pl和acpype.py脚本.这样就可以得到小分子的gro和top文件。

gro文件GMX可以识别，但top文件GMX还是不能识别，需要自己手动改动。

需要对GMX的立场和参数文件有一定的熟悉度。

最后得 gro与itp文件。

实现命令如下：a.产生小分子的高斯输入文件：（假设小分子名字为lig.mol2)antechamber -i lig.mol2 -fi mol2 -o -fogcrt -pf y -gm "%mem=4096MB" -gn "%nproc=8" -nc 1 -gk "#HF/6-31G* SCF=tight Test Pop=MK iop(6/33=2) iop(6/42=6)"b.用高斯计算小分子的静电势：nohup g09 > lig.log &c.把高斯的输出结果做resp电荷拟合：antechamebr -ilig.out -fi gout -o lig.mol2 -fo mol2 -c resp -rn LIGd.小分子参数的检查与坐标改动（因为高斯计算结果不会改动小分子的构想，但会改变小分子的移动，即相对运动。

消除相对移动只需要把坐标改成原始坐标即可，然后再产生GAFF缺失的参数文件）parmchk -i lig.mol2 -f mol2 -o lig.frcmodlig.frcmod即是小分子的缺失的参数文件e.生成小分子的amber参数文件（lig.prmtop与lig.inpcrd）先写个leap.in脚本文件，内容如下：source leaprc.ff99SBsource leaprc.gaffloadamberparamslig.frcmodlig=loadmol2 lig.mol2check ligsaveamberparmlig.prmtoplig.inpcrdsavepdblig lig.pdbquit然后运行tleap命令即可：tleap -f leap.in这样就拿到了小分子的参数文件（lig.prmtop与lig.inpcrd）f.调用脚本程序（amb2gmx.pl和acpype.py），尽量用acpype.py,比较稳定，需要安装python2.7版本环境，才能运行。

第62卷 第1期吉林大学学报(理学版)V o l .62 N o .1 2024年1月J o u r n a l o f J i l i nU n i v e r s i t y (S c i e n c eE d i t i o n )J a n 2024d o i :10.13413/j .c n k i .jd x b l x b .2023170新型靶向αv β6多肽的设计㊁合成与亲和力评价李跃鹏,王远强(重庆理工大学药学与生物工程学院,重庆400054)摘要:通过计算机辅助药物设计策略及固相合成法设计并筛选具有全新结构的αv β6多肽配体,并用酶联免疫吸附法(E L I S A )测定多肽配体与αv β6的结合亲和力,建立αv β6多肽配体的筛选方案.首先用S y b y l -X 1.3对αv β6多肽配体及天然配体进行分子对接;其次用A m b e r 16对分子对接结果进行分子动力学模拟,确定多肽配体与αv β6蛋白之间的结合模式,并以R G D L X X L (X 为任意氨基酸)为多肽配体核心结构,通过逐步延伸氨基酸构建虚拟肽库,筛选获得长度为7~10个氨基酸的多肽配体,进一步通过相似氨基酸替换的方法设计筛选不同于R G D L X X L 核心的新的多肽配体;最后用固相合成法合成新设计的多肽配体,利用间接E L I S A 法测定多肽配体-αv β6的结合亲和力.已有多肽配体和天然配体的分子对接以及分子动力学模拟结果表明,αv β6与配体的结合主要通过A s p 218和多肽配体之间形成氢键,以及M g 2+和多肽配体形成金属螯合作用完成.结合虚拟组合筛选与相似氨基酸替换,发现G R T D L G T L L F R ,G R R T D L A T I H G ,R T D V G R V R G R G 和R G D V G R V G R 等多肽均满足该结合模式,R T D V G R V R G R G 与αv β6的亲和力为10.76μm o l /L .可见R T D V G R V R G R G 与αv β6亲和力良好,是一条全新的αv β6多肽配体.关键词:整合素αv β6;多肽配体;药物设计;亲和力;酶联免疫吸附法中图分类号:R 914.2 文献标志码:A 文章编号:1671-5489(2024)01-0181-08D e s i g n ,S y n t h e s i s a n dA f f i n i t yE v a l u a t i o no f N o v e l T a r g e t e d αv β6P o l y p e p t i d e s L IY u e p e n g ,WA N G Y u a n q i a n g(S c h o o l o f P h a r m a c y a n dB i o e n g i n e e r i n g ,C h o n g q i n g U n i v e r s i t y o f T e c h n o l o g y ,C h o n g q i n g 400054,C h i n a )收稿日期:2023-05-02.第一作者简介:李跃鹏(1997 ),男,汉族,硕士研究生,从事药物设计和生物信息学的研究,E -m a i l :l i y u e p e n g l m c @f o x m a i l .c o m.通信作者简介:王远强(1979 ),男,汉族,博士,教授,从事靶向药物筛选与评价和多肽抗菌药物开发的研究,E -m a i l :w a n g y q n n @c q u t .e d u .c n .基金项目:重庆市科学技术局自然科学基金(批准号:C S T B 2022N S C Q -M S X 1327)和重庆市人力资源和社会保障局留学人员回国创业创新支持计划项目(批准号:c x 2020012).A b s t r a c t :W e d e s i g n e d a n ds c r e e n e d αv β6p o l y p e p t i d e l i g a n d sw i t hn e ws t r u c t u r eb y u s i n g c o m p u t e r -a i d e dd r u g d e s i g ns t r a t e g y a n ds o l i d -p h a s es y n t h e s i s m e t h o d ,a n dd e t e r m i n e dt h eb i n d i n g a f f i n i t y b e t w e e n t h e p o l y p e p t i d e l i g a n d s a n d αv β6b y u s i n g e n z y m e -l i n k e d i mm u n o s o r b e n t a s s a y (E L I S A )t o e s t a b l i s ha s c r e e n i n g s c h e m e f o r αv β6p o l y p e p t i d e l i g a n d s .F i r s t l y ,S y b y l -X1.3w a su s e d t o p e r f o r m m o l e c u l a r d o c k i n g o f t h e αv β6p o l y p e p t i d e l i g a n d s a n dn a t u r a l l i g a n d s .S e c o n d l y ,A m b e r 16w a s u s e d t o p e r f o r m m o l e c u l a rd y n a m i c ss i m u l a t i o nt od e t e r m i n et h eb i n d i n g m o d eb e t w e e nt h e p o l y p e p t i d e l i g a n d s a n d αv β6p r o t e i n ,a n dR G D L X X L (X w a sa n y a m i n oa c i d )w a su s e da s t h ec o r es t r u c t u r eo f t h e p o l y p e p t i d el i g a n d ,t h ev i r t u a l p e p t i d el i b r a r y w a sc o n s t r u c t e db ygr a d u a le x t e n s i o no fa m i n o281吉林大学学报(理学版)第62卷a c i d s,p o l y p e p t i d e l i g a n d sw i t hal e n g t ho f7 10a m i n oa c i d sw e r es c r e e n e d,a n dn e w p o l y p e p t i d el i g a n d sd i f f e r e n tf r o m t h ec o r eo f R G D L X X L w e r ed e s i g n e da n ds c r e e n e d b y s i m i l a ra m i n oa c i d s u b s t i t u t i o nm e t h o d.F i n a l l y,t h e n e w l y d e s i g n e d p o l y p e p t i d e l i g a n d sw e r e s y n t h e s i z e db y s o l i d-p h a s e s y n t h e s i sm e t h o d,a n dt h eb i n d i n g a f f i n i t y o f p o l y p e p t i d el i g a n d-αvβ6w a sd e t e r m i n e db y i n d i r e c t E L I S A m e t h o d.T h e m o l e c u l a rd o c k i n g a n d m o l e c u l a rd y n a m i c ss i m u l a t i o nr e s u l t so f p o l y p e p t i d e l i g a n d s a n dn a t u r a l l i g a n d s s h o wt h a t t h e b i n d i n g o fαvβ6t o t h e l i g a n d i sm a i n l y a c h i e v e d t h r o u g h t h e h y d r o g e nb o n df o r m a t i o n b e t w e e n A s p218a n dt h e p o l y p e p t i d el i g a n d,a n dt h e m e t a lc h e l a t i o n b e t w e e n M g2+a n dt h e p o l y p e p t i d e l i g a n d.C o m b i n e dw i t hv i r t u a l c o m b i n a t i o ns c r e e n i n g a n ds i m i l a r a m i n oa c i ds u b s t i t u t i o n,w ef i n dt h a t p o l y p e p t i d e ss u c h a s G R T D L G T L L F R,G R R T D L A T I H G, R T D V G R V R G R G a n d R G D V G R V G R a l l m e e tt h i s b i n d i n g p a t t e r n,a n d t h e a f f i n i t y b e t w e e n R T D V G R V R G R Ga n dαvβ6i s10.76μm o l/L.T h e r e f o r e,R T D V G R V R G R Ga n dαvβ6h a v ea g o o d a f f i n i t y a n d a r e an e wαvβ6p o l y p e p t i d e l i g a n d.K e y w o r d s:i n t e g r i nαvβ6;p o l y p e p t i d e;d r u g d e s i g n;a f f i n i t y;e n z y m e-l i n k e d i m m u n o s o r b e n t a s s a y(E L I S A)整合素是细胞黏附分子家族中的重要成员之一,是一组跨膜糖蛋白受体,广泛分布于细胞表面.整合素的主要功能是介导细胞与细胞㊁细胞与细胞外基质之间的相互黏附,并介导细胞与细胞外基质之间的双向信号传导,对细胞的黏附㊁增殖㊁转移㊁凋亡起重要的调控作用.整合素是由α亚基和β亚基通过非共价键组成的跨膜异二聚体糖蛋白,目前已发现18种α亚基和8种β亚基,可形成24种不同的整合素受体[1].其中,αvβ6是唯一由β6亚基参与形成的异源二聚体整合素.整合素αvβ6在健康的成人上皮中不表达或表达很弱,但在创伤愈合㊁纤维化和炎症等生理或病理过程中,甚至在一些恶性肿瘤中,它的表达水平显著上调,该特征使整合素αvβ6在各种疾病诊断与治疗中成为一个非常有吸引力的靶点[2-3].目前文献报道的靶向αvβ6配体包括小分子配体和多肽类配体,其中多肽类配体主要包括定向工程化噬菌体展示获得的多肽R01和S02[4],基于向日葵胰蛋白酶抑制剂噬菌体展示筛选得到的多肽S F L A P3[5],源于口蹄疫病毒外壳的多肽A20F M D v2[6],噬菌体展示文库筛选得到的多肽D L X X L[7],以及天然配体P r o-T G F-β1和T G F-β3[8].其中[18F]氟苯甲酰基标记的肽[18F]F B A-A20F M D v2已在人类临床试验中用作P E T放射性示踪剂,用于特发性肺纤维化的诊断和治疗评估.研究表明,尽管多肽配体没有直接的抗癌活性,但可以作为分子成像载体用于癌症等疾病诊断,以及通过与药物或纳米粒的偶联进行整合素靶向肿瘤治疗.靶向整合素放射药物已是比较成熟的研究方向,此类药物结合了两种关键元素:一种靶向化合物/配体和一种放射性同位素[9].靶向放射性药物可以与肿瘤组织中的特异性靶点结合并聚集,通过放射性射线产生电离辐射作用,由于正常组织细胞与肿瘤组织细胞对射线的敏感性不同,因此可以实现肿瘤的靶向治疗或根据放射性射线对肿瘤进行精准诊断.如177L u-P S MA-617(N o v a r t i s)即将被批准用于前列腺癌的诊疗[10],由北京大学王凡课题组研发的中国首个核医学肿瘤显像诊断1类新药99m T c-3P R G D2已完成三期临床[11].因此开发具有靶向αvβ6作用的多肽配体具有重要意义.针对目前αvβ6结合多肽基本含有R G D L X X L(X为任意氨基酸)核心结构,且该结构处于专利保护范围内,本文拟筛选设计非专利范围内全新结构的αvβ6结合多肽配体.先通过分子模拟方法分析αvβ6与多肽配体的结合模式,再通过计算机辅助药物设计方法发现一种全新结构的αvβ6多肽配体,该配体具有良好的结合亲和力,可为后续的靶向αvβ6多肽和多肽-核素复合物研究提供理论与实际研究指导,是一个具有发展前景的先导化合物.1材料与方法1.1试剂和仪器F m o c-氨基酸(上海麦克林生化科技股份有限公司),2-氯代三苯甲基氯-树脂(阿拉丁试剂(上海)有限公司),二氯甲烷(D C M )㊁二甲基甲酰胺(D M F )㊁二异丙基乙胺(D I E A )和乙腈等(国药集团化学试剂有限公司),人整合素αv β6蛋白(M e d C h e m E x p r e s s L L C )㊁链霉亲和素-辣根过氧化物酶(S t r e p t a v i d i n -H P R )和3,3ᶄ,5,5ᶄ-四甲基联苯胺(T M B )显色液(上海碧云天生物技术有限公司),A 20F M D v 2(N A V P N L R G D L Q V L A Q K V A R T )㊁A 20F M D v 2-B i o t i n (B i o t i n -N A V P N L R G D L Q V L A Q -K V A R T )和A 20F M D v 2-G R D (N A V P N L G R D L Q V L A Q K V A R T )(核欣(苏州)医药科技有限公司),酶联免疫吸附试剂盒(江苏酶免实业有限公司).液相制备色谱仪(L C -20A P 型,日本岛津公司);三重四级杆液-质联用仪(A G I L E N T6470型,美国安捷伦科技有限公司);全波长酶标仪(M u l t i s k a nS k y H i g h 型,美国赛默飞世尔科技公司).1.2 分子模拟方法1.2.1 分子对接采用S y b y l -X 1.3进行分子对接,靶蛋白αv β6源自蛋白质结构数据库(P r o t e i n D a t aB a n k ,h t t p s ://w w w.r c s b .o r g ,I D :5F F O ),配体来源于文献㊁有机小分子生物活性数据库(P u b C h e m ,h t t p s ://p u b c h e m.n c b i .n l m.n i h .g o v )和药物化学数据库(C h e m b l ,h t t p s ://w w w .e b i .a c .u k /c h e m b l /)等.使用S y b y l 软件的S u r f l e x -D o c k 模块,将蛋白晶体结构中天然配体所在位置定义为活性结合口袋.利用S u r f l e x -D o c kG e o m X (S F X C )模块进行高精度分子对接,分子对接参数:配体分子初始构像数目为10,分子片段最大构像数目为20,分子最大可旋转键数据为100,其他参数默认.综合结合构像和对接打分构建αv β6-多肽复合物用于分子动力学模拟.1.2.2 分子动力学模拟和结合自由能计算本文所有分子动力学模拟均在L i n u x 工作站用A m b e r 16[12]计算.αv β6使用F F 14S B 力场,多肽配体使用G A F F 力场和AM 1-B C C 电荷,将多肽-蛋白体系浸没于0.15m o l /L 氯化钠溶液(T I P 3P )的立方体盒子中,复合物距立方体盒子边缘为1.00n m ,体系包含98个N a +,78个C l -,36206个水分子.模拟过程如下:首先进行5步能量优化,避免可能的分子碰撞;其次进行两步加热,使体系温度从0逐渐升温到303.15K ,并进行溶液密度调整和体系平衡;最后在303.15K 下对系统(N P T )进行100n s 的动力学模拟,时间步长为2f s,压力恒定为1个大气压.结合自由能以及能量分解计算由A m b e r 的MM P B S A.p y 程序完成[13],选取最后20n s 进行计算并取其平均结构进行结合模式分析,分析αv β6与配体的结合模式.当结合自由能的值为负值时,体系是稳定的,且该值越小配体-受体结合亲和力越高.1.2.3 组合虚拟筛选本文虚拟筛选由S Y B Y L -X1.3完成.首先以R G D L X X L 为基础,在5,6号位通过20种天然氨基酸的随机组合建立一个含400条多肽的多肽文库,通过R o s e t t a 软件[14]对多肽进行批量结构优化构建虚拟肽库,根据1.2.1节中定义的活性口袋进行虚拟筛选,综合结合构像以及对接打分筛选出最佳结合七肽.然后将筛选出的七肽作为核心逐步延长多肽链,同理筛选出最佳结合八肽㊁九肽和十肽用于进一步研究.结合模式分析结果表明,现有多肽配体与αv β6的结合主要为九肽和十肽核心结构,因此将筛选出的九肽㊁十肽和骨架R G D L X X L 进行氨基酸替换,如亮氨酸可替换为缬氨酸,二者等电点相近,分别为5.98和5.96,侧链结构相似,分别为异丙基和异丁基,二者均为带疏水性侧链的脂肪族氨基酸,氨基酸替换时在改变多肽结构的同时尽量不改变其理化性质.在保持多肽关键结合氨基酸不变的情况下,进一步构建虚拟肽库,通过与上一步相同的虚拟筛选㊁精准对接㊁分子动力学模拟获得最佳结合多肽配体,利用分子动力学模拟和结合自由能能量分解计算验证结构改造的合理性.1.3 多肽固相合成方法用D C M 激活树脂后,通过D I E A 将第一个氨基酸耦联到树脂上,然后用哌啶/D M F 保护第一个氨基酸,通过茚三酮实验确认氨基酸耦合后用D M F 和甲醇洗涤树脂.使用预活化的第二个氨基酸㊁肽耦合试剂H B T U 和D I E A 溶液进行耦联,确认耦联后,重复洗涤,重复去保护和耦合循环,直到获得所需的肽链长度.最后将树脂结合的多肽转移到圆底瓶中,通过n (T F A )ʒn (H 2O )ʒn (E D T )ʒn (T I S )=94.5ʒ2ʒ2.5ʒ1的组合溶液同时去除树脂和保护基团.再利用制备型高效液相色谱结合梯381 第1期 李跃鹏,等:新型靶向αv β6多肽的设计㊁合成与亲和力评价481吉林大学学报(理学版)第62卷度洗脱对多肽进行分离纯化,其中乙腈/水为流动相,检测波长为214n m.最后用液质联用仪(L C-M S)确认目标产物.1.4E L I S A法测定结合亲和力人源αvβ6蛋白经包被㊁封闭后,在孔板中加入待测配体溶液.其中空白组仅加入100μL T B S缓冲溶液;空白对照中加入A20F M D v2-B i o t i n和T B S缓冲溶液各50μL;阳性对照加入A20F M D v2-B i o t i n和A20F M D v2(梯度)各50μL;阴性对照加入A20F M D v2-B i o t i n和A20F M D v2-G R D各50μL;实验组加入A20F M D v2-B i o t i n和待测多肽(梯度)各50μL.每组设置3个平行试验.反应2h后, S t r e p t a v i d i n-H R P用P B S T稀释100倍,向各孔加入100μL,室温振摇反应后,向各孔加入100μLT M B显色液,立即将孔板置于酶标仪中避光孵育,每5m i n读取一次650n m波长处的吸光度(A)值,共读取1h(13个检测点).2结果与讨论2.1αvβ6与现有配体的结合模式现有αvβ6多肽配体的分子对接与结合自由能计算结果列于表1.由表1可见,分子对接打分选取结合构像最佳的对接结果,所有配体分子与αvβ6的结合自由能为-147.7534~-92.8887k J/m o l.结合自由能能量分解结果表明,αvβ6与多肽配体形成的基本强相互作用:A r g R G D侧链通过氢键与αv 亚基A s p218上的羧基形成氢键;多肽的侧链残基与β6亚基种金属离子依赖的黏附位点中M g2+形成配位键;配体中-L X X L-部分可以形成两亲性α螺旋,并且该螺旋结构与β6亚基形成的特异性疏水口袋结合,该结合模式为后续设计的多肽是否合理建立标准.表1现有αvβ6配体分子对接以及结合自由能计算结果T a b l e1C a l c u l a t i o n r e s u l t s o fm o l e c u l a r d o c k i n g a n db i n d i n g f r e e e n e r g y f o r e x i s t i n gαvβ6l i g a n d多肽配体多肽序列分子对接打分结合自由能/(k J∙m o l-1) R01G C I L NM R T D L G T L L F R C R R D S D C P G A C I C R G N G Y C G10.65-136.2091S02G C R S L A R T D L D H L R G R C S S D S D C L A E C I C L E N G F C G9.36-131.6895S F L A P3G R C T G R G D L G R L C Y P D9.03-147.7534A20F M D v2N A V P N L R G D L Q V L A Q K V A R T11.39-142.9696D L X X L R T D L D S L R T7.46-92.8887T G F-β3H R G D L G R L K K8.89-139.6147P r o-T G F-β1N G F T T G R R G D L A T I H GMN R8.53-144.51582.2靶向αvβ6多肽配体的设计以R G D L X X L为核心初步筛选多肽配体的分子对接以及结合自由能计算结果列于表2.结合自由能除七肽较低为-26.3186k J/m o l外,其他3条多肽结合自由能均小于-62.8020k J/m o l.虚拟筛选的多肽与αvβ6结合模式如图1所示.由图1可见,七肽除必要的结合作用外,还与T y r178形成键长为0.29,0.39n m的氢键,并与G l u698形成键长为0.35n m的氢键;八肽除必要的结合作用外,还与A s n693和S e r602分别形成键长为0.38,0.31n m的氢键;九肽除必要的结合作用外,还与T y r178, A s p729,T h r696,A s p695分别形成键长0.32,0.32,0.40,0.38n m的氢键;十肽除必要的结合作用外,还与T y r178,S e r602,T h r696分别形成键长为0.40,0.32,0.39n m的氢键.表2多肽虚拟筛选㊁设计分子对接以及结合自由能计算结果T a b l e2V i r t u a l s c r e e n i n g o f p o l y p e p t i d e s,d e s i g no fm o l e c u l a r d o c k i n g a n d c a l c u l a t i o n r e s u l t s o f b i n d i n g f r e e e n e r g y多肽序列分子对接打分结合自由能/(k J∙m o l-1)七肽R G D L T R L8.23-26.3186八肽R G D L T R L G8.64-70.9813九肽R G D L T R L G R9.49-104.0541十肽R G D L T R L G R T9.62-116.2197图1 虚拟筛选的多肽与αv β6结合模式F i g .1 B i n d i n gp a t t e r n s o f v i r t u a l s c r e e n i n g o f p o l y p e p t i d e s a n d αv β6 考虑到现有研究中的配体均基于R G D L X X L 配体,本文对初步筛选出的九肽和十肽的氨基酸残基进行替换筛选,探索全新结构的αv β6多肽配体,最终得到2条多肽,分别为R G D V G R V G R 和R T D V G R V R G R .2条多肽配体与αv β6结合模式如图2所示.由图2可见,2条多肽与αv β6的结合模式与现有配体一致:多肽R T D V G R V G R G 除必要的结合作用外,还与A s p 150,A s p 148,I l e 147,G l u 121,S e r 656,G l u 698分别形成键长为0.29,0.34,0.35,0.31,0.32,0.32n m 的氢键;多肽R G D V G R V G R 除必要的相互作用外,还与T y r 178,A s p 695,T h r 696,L y s 119,S e r 656,P r o 654,L y s 645分别形成键长为0.38,0.27,0.23,0.33,0.29,0.29,0.28n m 的氢键.图2 氨基酸替换后多肽结合模式F i g .2 P o l y p e p t i d e b i n d i n gpa t t e r na f t e r a m i n o a c i d s ub s t i t u t i o n 2.3 多肽合成结果本文合成了4条多肽,分别为P 1(G R T D L G T L L F R ),P 2(G R R T D L A T I H G ),P 3(R T D V G R V R G R ),P 4(R G D V G R V G R ).其中P 2为p r o -T G F -β1配体的R G D 特征片段[15-17],在后续研究中作为设计多肽亲和力测定的对照,P 1为重庆理工大学靶向药物筛选与活性评价课题组原有肽库筛选出的配体,P 3和P 4为本文结构优化后的多肽配体.根据高效液相色谱(H P L C )结果,P 1纯度为98.461%,P 2纯度为95.749%,P 3纯度为95.135%,P 4纯度为97.554%.根据L C -M S 结果,P 1计算的摩尔质量为1248.42g /m o l ,测量的摩尔质量为1249.00g /m o l ;P 2计算的摩尔质量为1196.31g /m o l ,测量的摩尔质量为1196.70g /m o l ;P 3计算的摩尔质量为1171.31g /m o l ,测量的摩尔质量为1171.35g /m o l ;P 4计算的摩尔质量为971.07g /m o l ,测量的摩尔质量为971.25g /m o l .可见,经H P L C 以及L C -M S 分析,确定合成的多肽581 第1期 李跃鹏,等:新型靶向αv β6多肽的设计㊁合成与亲和力评价为本文研究的目标产物.2.4亲和力测定结果酶联免疫吸附(E L I S A)实验采用饱和浓度法进行计算[18],即体系反应达到最大响应值的50%时的浓度为结合亲和力浓度.每个孔位设置3组平行试验,测量数据结果显示标准差不超过平均值的15%.测试样品吸光度在进行计算时应减去标准品浓度为0的A平均值,即只加了T B S缓冲液的孔位吸光度.采用A20F M D v2作为阳性对照(亲和力3μm o l/L),利用O r i g i n2016对标准品吸光度及浓度绘制四参数逻辑函数曲线,计算阳性对照达到最大反应响应值的浓度约为40μm o l/L,如图3所示.利用B i o t i n标记的A20F M D v2作为检测目标分子,S t r e p t a v i d i n-H R P可与B i o t i n-A20F M D v2形成稳定复合物,当待测多肽与αvβ6结合时溶液中的B i o t i n-A20F M D v2浓度增加,T M B可与S t r e p t a v i d i n-H R P生成蓝色化合物进而检测溶液吸光度.根据4条待测多肽的A值以及浓度㊁阳性标准品的吸光度,绘制多肽-蛋白结合率和多肽浓度相关性曲线,结果如图4所示.由图4和吸光度数据可以计算出P1的亲和力为3.35μm o l/L,P3的亲和力为10.76μm o l/L .图3阳性对照A20F M D v2的浓度-吸光度曲线F i g.3C o n c e n t r a t i o n-a b s o r b a n c e c u r v eo f p o s i t i v e c o n t r o lA20F M D v 2图4待测多肽的浓度-结合率曲线F i g.4C o n c e n t r a t i o n-b i n d i n g r a t e c u r v e so f t e s t e d p o l y p e p t i d e2.5讨论本文通过分子对接与分子动力学模拟分析αvβ6与配体的结合模式,虚拟组合筛选获得符合结合模式的全新结构多肽,采用固相合成法合成多肽,通过间接E L I S A法测定多肽配体与αvβ6的结合亲和力.实验结果表明:多肽配体P3与αvβ6具有良好的结合亲和力,略大于阳性对照A20F M D v2的结合亲和力,基本符合设计预期.实验中P2和P4由于亲和力较弱,仅存在阳性结果,无法计算其亲和力值.P1虽然也有较强的αvβ6蛋白结合亲和力,但其结构仍具有R G D L X X L核心,因此在后续的研究中可作为先导化合物改造其结构.P3和P4虽都为氨基酸替换后筛选出的多肽,但其结合亲和力差距较大.综合分子对接和结合自由能相关计算结果,推测P4活性较差的原因是其-V G R V G R-部分的空间位阻较小,而αv与β6亚基形成的结合空腔较大,二者结合时存在一定的结合不稳定性.结构生物学研究表明,αvβ6不仅能识别R G D序列,还能识别-L X X L-基序,该基序折叠成两亲性α-螺旋,结合至β6亚基残基组成的疏水口袋中[15-17].研究中筛选出全新结构多肽含有的-V X X V-序列也可以形成两亲性α-螺旋,该螺旋可与β6亚基形成更多氢键相互作用,R T D序列也可以与αvβ6形成稳定的氢键,表明了本文采用氨基酸替换方法的合理性.可见,计算机辅助药物设计方法可有效提高药物开发的成功率㊁降低研发成本并缩短研发周期,是创新药物研发的重要方法.文献研究表明,整合素αvβ6的表达在许多肿瘤中显著上调,包括口腔鳞癌㊁乳腺癌㊁胃癌㊁胰腺癌㊁结直肠癌㊁肺癌等[18-20].整合素αvβ6的表达与许多癌种患者的生存率降低相关,例如结直肠癌㊁乳腺癌㊁胰腺导管腺癌㊁非小细胞肺癌和宫颈鳞状细胞癌等[21-25],并且可以通过纤维连接蛋白和胶原中的R G D序列参与肿瘤细胞与胞外基质的相互作用,该整合素的过度表达与上皮细胞向间充质样细胞转化有关,可促进肿瘤细胞的侵袭和迁移,通常在肿瘤侵袭病灶的前沿部分中整合素的表达水平较高.目前,靶向整合素成像已成为整合素研究热点之一,包括正电子发射断层扫描和单光子发射计算机断层扫描,整合素显像剂可用于对患者分层进行靶向治疗㊁评估治疗反应,并监测肿瘤生长[26-30].681吉林大学学报(理学版)第62卷由于大部分基础生物学仍未完全了解,因此αv β6为靶标的治疗肿瘤药物也尚未开发,但其仍有作为诊断靶点的潜力.综上所述,本文采用计算机辅助药物设计的方法设计了一条基于R T D V X X V 全新结构的αv β6多肽配体,研究结果表明,该配体具有良好的结合亲和力,配体分子可作为靶向αv β6多肽或肽类结合物开发的候选分子,为未来αv β6结合配体的开发提供结构基础.参考文献[1] S L A C K R J ,MA C D O N A L D SJ F ,R O P E R J A ,e ta l .E m e r g i n g T h e r a p e u t i c O p p o r t u n i t i e sf o rI n t e gr i n I n h i b i t o r s [J ].N a tR e vD r u g D i s c o v ,2022,21(1):60-78.[2] K O I V I S T OL ,B I JR ,HÄK K I N E NL ,e t a l .I n t e g r i n αv β6:S t r u c t u r e ,F u n c t i o n a n dR o l e i nH e a l t h a n dD i s e a s e [J ].I n t JB i o c h e m C e l l B i o l ,2018,99:186-196.[3] B R Z O Z OW S K A E ,D E S HMU K H S .I n t e g r i n A l p h avB e t a6(αv β6)a n dI t s I m p l i c a t i o n s i nC a n c e rT r e a t m e n t [J ].I n t JM o l S c i ,2022,23(20):12346-1-12346-18.[4] K I MU R A R H ,T E E DR ,HA C K E LBJ ,e t a l .P h a r m a c o k i n e t i c a l l y S t a b i l i z e dC y s t i n eK n o tP e pt i d e sT h a tB i n d A l p h a -v -B e t a -6I n t e g r i nw i t hS i n g l e -D i g i tN a n o m o l a rA f f i n i t i e s f o rD e t e c t i o n o f P a n c r e a t i cC a n c e r [J ].C l i nC a n c e r R e s ,2012,18(3):839-849.[5] R OE S C H S ,L I N D N E RT ,S A U T E R M ,e t a l .C o m p a r i s o no f t h eR G D M o t i f -C o n t a i n i n g αv β6I n t e g r i n -B i n d i n g P e p t i d e sS F L A P 3a n d S F I T G v 6f o r D i a g n o s t i c A p pl i c a t i o ni n HN S C C [J ].J N u c l M e d ,2018,59(11):1679-1685.[6] HA U S N E RSH ,D I C A R AD ,MA R I KJ ,e t a l .U s e o f aP e p t i d eD e r i v e d f r o mF o o t -a n d -M o u t hD i s e a s eV i r u s f o r t h eN o n i n v a s i v e I m a g i n g o fH u m a nC a n c e r :G e n e r a t i o n a n dE v a l u a t i o n o f 4-[18F ]F l u o r o b e n z o yl A 20F M D v 2f o r i n v i v o I m a g i n g o f I n t e g r i n A l p h a v b e t a 6E x p r e s s i o n w i t hP o s i t r o nE m i s s i o nT o m o g r a p h y [J ].C a n c e rR e s ,2007,67(16):7833-7840.[7] K R A F TS ,D I E F E N B A C H B ,M E H T A R ,e ta l .D e f i n i t i o no fa n U n e x p e c t e d L i g a n d R e c o gn i t i o n M o t i ff o r A l p h a vB e t a 6I n t e g r i n [J ].JB i o l C h e m ,1999,274(4):1979-1985.[8] D O N G XC ,HU D S O N N E ,L U CF ,e t a l .S t r u c t u r a lD e t e r m i n a n t s o f I n t e g r i n β-S u b u n i t S p e c i f i c i t y f o rL a t e n t T G F -β[J ].N a t S t r u c tM o l B i o l ,2014,21(12):1091-1096.[9] G A R O U S I J ,V O R O B Y E V A A ,A L T A IM.I n f l u e n c e o f S e v e r a l C o m p o u n d s a n dD r u g s o n t h eR e n a lU pt a k eo f R a d i o l a b e l e dA f f i b o d y M o l e c u l e s [J ].M o l e c u l e s ,2020,25(11):2673-1-2673-10.[10] S A R T O R O ,D E B O N O J ,C H I K N ,e ta l .L u t e t i u m -177-P S MA -617f o r M e t a s t a t i c C a s t r a t i o n -R e s i s t a n t P r o s t a t eC a n c e r [J ].N E n g l JM e d ,2021,385(12):1091-1103.[11] Z HUZ H ,M I A O W B ,L IQ W ,e ta l .99m T c -3P R G D 2f o rI n t e g r i n R e c e p t o rI m a g i n g o fL u n g C a n c e r :A M u l t i c e n t e r S t u d y [J ].JN u c lM e d ,2012,53(5):716-722.[12] C A S ED A ,C H E A T HAM ⅢT E ,D A R D E N T ,e t a l .T h eA m b e rB i o m o l e c u l a rS i m u l a t i o nP r o gr a m s [J ].J C o m p u tC h e m ,2005,26(16):1668-1688.[13] M I L L E R ⅢBR ,M C G E ETD ,J r ,S WA I L S JM ,e t a l .MM P B S A.p y :A nE f f i c i e n t P r o g r a mf o r E n d -S t a t e F r e e E n e r g y C a l c u l a t i o n s [J ].JC h e m T h e o r y C o m p u t ,2012,8(9):3314-3321.[14] N E R L I S ,S G O U R A K I SN G.C s -R o s e t t a [J ].M e t h o d sE n z y m o l ,2019,614:321-362.[15] C AM P B E L L M G ,C O R M I E R A ,I T O S ,e t a l .C r y o -E M R e v e a l s I n t e g r i n -M e d i a t e dT G F -βAc t i v a t i o nw i t h o u t R e l e a s e f r o m L a t e n tT G F -β[J ].C e l l ,2020,180(3):490-501.[16] S C H I N N E R C ,X U L F ,F R A N Z H ,e t a l .D e f e c t i v e D e s m o s o m a l Ad he s i o n C a u s e s A r r h y t h m o ge n i c C a r d i o m y o p a t h y b y I n v o l v i n g a n I n t e g r i n -αv β6/T G F -βSi g n a l i n g C a s c a d e [J ].C i r c u l a t i o n ,2022,146(21):1610-1626.[17] WA N G X ,L I US ,Y U T ,e t a l .I n h i b i t i o no f I n t e g r i n αv β6A c t i v a t i o no fT G F -βAt t e n u a t e sT e n d i n o p a t h y [J ].A d vS c i (W e i n h ),2022,9(11):e 2104469-1-e 2104469-14.[18] 陈文龙,张阳阳,张生英,等.牛病毒性腹泻病毒E 2蛋白的克隆㊁表达及多克隆抗体制备[J ].吉林大学学报(理学版),2020,58(3):711-717.(C H E N W L ,Z HA N G Y Y ,Z HA N GSY ,e t a l .C l o n i n g a n dE x pr e s s i o no f E 2P r o t e i no f B o v i n e v i r a l d i a r r h e av i r u s a n dP r e p a r a t i o no fP o l y c l o n a lA n t i b o d y [J ].J o u r n a l o f J i l i nU n i v e r s i t y 781 第1期 李跃鹏,等:新型靶向αv β6多肽的设计㊁合成与亲和力评价881吉林大学学报(理学版)第62卷(S c i e n c eE d i t i o n),2020,58(3):711-717.)[19] HU A N G H S,Y U A N M C,WU SL,e ta l.C l i n i c a lS i g n i f i c a n c eo fC-X-C M o t i fC h e m o k i n eR e c e p t o r4a n dI n t e g r i nαvβ6E x p r e s s i o n i nB r e a s tC a n c e r[J].JB r e a s tC a n c e r,2020,23(2):171-181.[20] B A A R T V M,V A N D U I J NC,V A NE GMO N DSL,e t a l.E G F Ra n dαvβ6a sP r o m i s i n g T a r g e t s f o rM o l e c u l a rI m a g i n g o fC u t a n e o u s a n d M u c o s a l S q u a m o u sC e l l C a r c i n o m a o f t h eH e a da n dN e c kR e g i o n[J].C a n c e r s,2020,12(6):1474-1-1474-13.[21] MO O R E K M,D E S A I A M,D E L G A D O B L,e ta l.I n t e g r i nαvβ6-S p e c i f i c T h e r a p y f o rP a n c r e a t i c C a n c e rD e v e l o p e d f r o m F o o t-a n d-M o u t h-D i s e a s eV i r u s[J].T h e r a n o s t i c s,2020,10(7):2930-2942.[22] L I A N GB J,L I L P,M I A O R Z,e ta l.E x p r e s s i o n o fI n t e r l e u k i n-6a n dI n t e g r i nαvβ6i n C o l o n C a n c e r:A s s o c i a t i o nw i t hC l i n i c a lO u t c o m e s a n dP r o g n o s t i c I m p l i c a t i o n s[J].C a n c e r I n v e s t,2019,37(3):174-184.[23] A HNSB,MOHAM E D A L IA,C HA N C,e ta l.C o r r e l a t i o n sb e t w e e nI n t e g r i nαvβ6E x p r e s s i o na n d C l i n i c o-P a t h o l o g i c a l F e a t u r e s i nS t a g eBa n dS t a g eCR e c t a l C a n c e r[J].P L o SO n e,2014,9(5):e97248-1-e97248-8.[24] D E S A IK,N A I R M G,P R A B HUJS,e t a l.H i g hE x p r e s s i o no f I n t e g r i nβ6i n A s s o c i a t i o n w i t ht h eR h o-R a cP a t h w a y I d e n t i f i e s aP o o rP r o g n o s t i cS u b g r o u p w i t h i n H E R2A m p l i f i e dB r e a s tC a n c e r s[J].C a n c e rM e d,2016, 5(8):2000-2011.[25] R AMO SD M,D A N GD M,S A D L E RS.T h eR o l e o f t h e I n t e g r i nA l p h a vB e t a6i nR e g u l a t i n g t h eE p i t h e l i a l t oM e s e n c h y m a lT r a n s i t i o n i nO r a l C a n c e r[J].A n t i c a n c e rR e s,2009,29(1):125-130.[26] B A T E SR C,B E L L O V I N DI,B R OWN C,e ta l.T r a n s c r i p t i o n a lA c t i v a t i o no fI n t e g r i n B e t a6d u r i n g t h eE p i t h e l i a l-M e s e n c h y m a lT r a n s i t i o n D e f i n e sa N o v e lP r o g n o s t i cI n d i c a t o ro f A g g r e s s i v e C o l o n C a r c i n o m a[J].JC l i n I n v e s t,2005,115(2):339-347.[27] L I O L I O SC,S A C H P E K I D I SC,K O L O C O U R I SA,e t a l.P E TD i a g n o s t i cM o l e c u l e sU t i l i z i n g M u l t i m e r i cC y c l i cR G DP e p t i d eA n a l o g s f o r I m a g i n g I n t e g r i nαvβ3R e c e p t o r s[J].M o l e c u l e s,2021,26(6):1792-1-1792-27. [28] D O L G I N E.R a d i o a c t i v e D r u g s E m e r g e f r o m t h e S h a d o w s t o S t o r m t h e M a r k e t[J].N a t B i o t e c h n o l,2018,36(12):1125-1127.[29] R I C E M,K R O S I N M,HA S T E P.P o s t Y t t r i u m-90I m a g i n g[J].S e m i nI n t e r v e n t R a d i o l,2021,38(4):460-465.[30] B A D I P AF,A L I R E Z A O U R B,Y O U S E F N I A H.A n O v e r v i e wo fR a d i o l a b e l e dR G D P e p t i d e s f o rT h e r a n o s t i cA p p l i c a t i o n s[J].C u r rR a d i o p h a r m,2022,16(2):107-122.(责任编辑:单凝)。

gaff力场分子模型
Gaff力场分子模型是一种计算机模拟工具，用于描述分子的
力场行为。

Gaff力场是由美国Cornell大学开发的一种力场模型，可用于描述各种有机和生物分子的分子力学行为。

Gaff
力场使用参数化方法，将分子中的原子和键与经验规则相匹配，并根据实验和计算结果进行调整。

Gaff力场分子模型基于分子力学的原理，通过计算力场势能
函数来描述分子的结构和行为。

力场势能函数由键能、角能、扭转能、非键能等多个项组成。

每个原子和键都有各自的力常数和参数，用于计算相应项的势能。

使用Gaff力场分子模型，可以模拟分子的构象、稳定性、物
理性质等多个方面。

它可以用于预测分子的结构、能量、振动模式等，对于药物研发、新材料设计等领域具有重要意义。

然而，由于力场模型是经验参数化的，所以其适用范围有限，可能存在一定的误差和不确定性。

总之，Gaff力场分子模型是一种常用的计算工具，可用于模
拟和预测分子的力场行为，对于研究和应用分子相关问题具有重要价值。

分子力场详细讲解分子力场（Molecular Force Field），简称FF，是描述分子体系中分子间相互作用和内部构型变化的数学模型。

它可以通过经验力场和基于量子化学计算的理论力场两种方式来建立。

本文将详细讲解分子力场的原理、构建方式以及在分子模拟和化学计算中的应用。

一、原理与目标分子力场的目标是通过描述原子之间的键长、键角以及相互作用力的形式，来预测分子的结构、稳定性和相对能量变化。

它的基本原理是将分子势能视为原子之间相互作用的总和，并通过参数化来拟合实验数据或量子化学计算结果。

二、常见参数分子力场的参数化包括键长、键角、二面角、扭曲能以及原子电荷等。

键长是相邻两个原子之间的距离，如C-C键长为1.54 Å。

键角是三个相邻原子所形成的夹角，如C-C-C键角为109.5°。

二面角是四个连续原子所形成的角度，如C-C-C-C二面角为180°。

扭曲能是分子内部原子之间由于旋转而引起的能量变化。

原子电荷用于描述原子之间的静电相互作用，一般可以通过量子化学计算来得到。

三、参数化方法1. 经验参数化：一种常用的方法是通过拟合实验数据来确定分子力场的参数。

例如，通过测量一系列分子的结构和能量，可得到不同键长、键角和二面角对应的能量差值。

然后采用数学方法进行拟合，从而获得各个参数的数值。

2. 理论参数化：基于量子化学计算的理论力场是另一种参数化方法。

通过量子化学软件计算分子的结构和能量，然后与实验数据进行对比，并通过优化参数得到最佳的拟合结果。

四、分子模拟与化学计算分子力场在分子模拟和化学计算中得到了广泛的应用。

1. 分子模拟（Molecular Dynamics，MD）模拟分子系统的动力学过程，通过数值求解牛顿运动方程来模拟分子的运动轨迹。

分子力场用于计算给定构型下分子的势能和受力矢量，并在模拟过程中改变原子的位置和速度。

根据分子力场的计算结果，可以得到分子的平衡构型、结构和能量的变化规律，进而研究分子的稳定性、反应动力学等性质。

分子动力学模拟研究活化剂对葡萄糖激酶活性的影响及活化机理姚雪霞【摘要】Glucokinase (GK) is a glycolic enzyme that catalyzes the phosphorylation of glucose to glucose-6-phosphate in the first step of glycolysis, and the change in GK activity plays an important role in abnormal glycemia. Thus molecular dynamics simulation was performed to study the effect of glucokinase activator (GKA) on the activity of GK, and the activation mechanism of GKA was examined. The binding free-energy analysis suggests that van der Walls forces are the main driving forces for GK-GKA complex, and glucose can enhance the binding between GK and GKA. Moreover, conformation analysis shows that GKA can constrain the conformation of GK in a period of time, but it cannot keep the active state of GK in an extended period, which is in agreement with previously reported experiment results. These observations are beneficial to understanding the activation mechanism of GKA at an atomic level and to guiding the design of novel GKA for the treatment of diabetes.%葡萄糖激酶(GK)催化葡萄糖转变为6-磷酸葡萄糖,是糖代谢的第一步；GK的活性改变对葡萄糖的代谢影响很大.通过分子动力学模拟研究了葡萄糖激酶活化剂(GKA)对GK活性的影响,探讨了GKA的活化机理.包结自由能分析表明,范德华作用是GK和GKA体系的主要驱动力；与此同时,葡萄糖使得GK和GKA的结合增强.构象分析结果表明,GKA能够在一段时间内限制GK的构象,但不能长久维持GK的活化构象不变,这和以往的实验结果相一致.相关理论分析结果不仅有利于从原子水平解释GKA的活化机理,也可为设计新的用于糖尿病治疗的GKA提供一定的理论参考.【期刊名称】《化学研究》【年(卷),期】2012(023)004【总页数】6页(P91-96)【关键词】葡萄糖激酶;活化剂;活性;活化机理;分子动力学模拟【作者】姚雪霞【作者单位】南京农业大学工学院,江苏南京210031【正文语种】中文【中图分类】O629.8葡萄糖激酶（GK）催化葡萄糖转变为6－磷酸葡萄糖，是糖代谢极为关键的一步，所以GK活性异常在糖代谢紊乱的发生发展中起着重要作用［1－2］.文献［3］对GK的晶体结构和催化过程进行了研究，认为GK存在三种构象，即：关闭型、开启型和超开启型.其中，前两种为活性构象，第3种为非活性构象.其催化过程如下：当血糖浓度较低时，GK主要以稳定的超开启构象存在，当葡萄糖分子与之结合后，其构象缓慢转变为开启型构象，然后转变为关闭型构象，在ATP存在下同时发生酶催化反应——葡萄糖磷酸化.反应结束后GK回复到开启型构象，同时释放出6－磷酸葡萄糖（G6P）和ADP.若有新的葡萄糖分子结合则继续发生磷酸化，进入快循环；若没有新的葡萄糖分子结合，GK则转变为超开启构象，进入慢循环.若GK的氨基酸残基发生失活突变，则影响相应的催化过程，导致青年型早发糖尿病（MODY 2）［4］发生.为了治疗MODY 2型糖尿病，以GK为靶点的小分子活化剂（GKA）的研究非常活跃，并取得了很大的进展［5－9］.观察GK的晶体结构［8］（如图1所示），发现GK由448个氨基酸残基构成，折叠成大小不等的两个结构域，中间凹槽为连接域，葡萄糖分子磷酸化的活性位点位于该连接域内，而GKA则位于GK的异构点，距离葡萄糖底物大约2nm左右.虽然实验已经证实GKA可以提高GK的活性，但是还是有很多基本问题没有解决.如，异构点距离葡萄糖底物的距离为2nm左右，GKA如何提高GK的活性？文献［10］虽然讨论了其中一个GKA对GK的影响，但是还有很多细节问题没有解决.如GKA是否一直对GK的构象起着约束作用？另外葡萄糖的存在反过来对GKA 有什么影响？以及其他GKA对GK的影响是否相同？为了解决这些疑问，很有必要对其他GKA对GK的活化机理做进一步的理论研究.以往我们曾分析GK的活化突变对GK活性的影响［11－12］，发现它们主要通过约束GK的构象来改变GK的活性.活化剂是否采用同样的方式来改变GK的活性呢？为了研究GKA的活化机理，作者通过分子动力学（molecular dynamical，MD）模拟和 MM－PBSA（molecular mechanics／Poisson Boltzmann surface area）方法［13－14］，对GKA、GK以及葡萄糖的相互作用进行了理论研究.1.1 蛋白质构建由于第一个包含GK、GKA和葡萄糖复合物的晶体结构［3］缺少实验数据EC50值，因此本文作者选用了一个新研制出来的晶体结构［8］作为研究对象.该体系和第一个体系相比，只有GKA不同，GKA的分子结构如图2所示.GK－GKA－glucose复合体系的初始结构取自于文献（PDB code：3FR0）［8］.目前单独的GK－glucose和单独的GK－GKA晶体结构还没有文献报道，因此单独的GK－glucose和GK－GKA只能从复合体系中去掉另一个小分子作为初始结构.另外本文中所指的GK－free是指直接从GK－GKA－glucose复合体系中去掉GKA和葡萄糖两个小分子得到的GK结构，并不是真正的晶体结构GK－free. 1.2 MD模拟MD模拟的初始结构如上所述.所有的MD模拟均采用AMBER 10程序包［15］.AMBER03力场用于蛋白质分子，GAFF力场用于葡萄糖和GKA小分子.葡萄糖和GKA小分子都采用Gaussian 03程序包［16］的HF／6－31G＊基组优化，采用RESP程序计算其电荷.程序LeaP用于中和体系，使其总电荷为零.对于所有体系，均在溶质外围加上1nm的水分子层.采用周期性边界条件和particle mesh ewald（PME）方法用于处理长程静电作用［17］.非键截距为1.2nm，使用SHAKE选项用于约束和氢原子连接的键.在运行MD之前，先进行两步能量优化：第一步约束溶质，用最速下降法优化5 000步，再用共轭梯度法优化5 000步；第二步去约束后再优化整个体系，用最速下降法优化2 500步，再用共轭梯度法优化2 500步，收敛条件为能量梯度小于4.18×10－4 kJ·mol－1·nm－1.MD模拟分为两步：首先进行200ps的约束溶质分子MD模拟，即通过200ps预热，将体系从0K加热到300K；然后进行10ns无约束恒温MD模拟，温度为300K，压力为1×105Pa（即NPT系综）.系统的温度和压力则通过Langevin piston和Hoover’s thermostat实现控制.模拟中采用的积分步长为2fs，非键连接表每50步更新一次.每4ps保存一次轨迹结构用于后期数据处理.2.1 MD轨迹稳定性的分析体系在模拟过程中产生的结构相对于晶体结构随时间变化的均方根偏差（RMSD）是衡量体系是否稳定和达到收敛的重要依据.为了检测MD轨迹的稳定性，分别对GK－free，GK－GKA，GK－glucose和 GK－GKA－glucose四个轨迹中GK的Cα原子进行RMSD分析，结果如图3所示.分析图3，发现相比于其他轨迹，GK －free的轨迹RMSD值最大，最不稳定.而GK－glucose和GK－GKA－glucose的RMSD值比较小，轨迹相对稳定一些.但是仔细分析，发现GK－GKA－glucose的轨迹最稳定.这表明GKA的存在，可以使得GK－glucose体系更加稳定.2.2 GKA和GK之间静电作用的分析分析其晶体结构（如图1所示），可以看到GKA处于GK背部的空腔内，不仅和GK有非键的范德华作用，而且还有氢键作用.晶体结构中出现的氢键，在动态的MD过程中会有什么变化？下面分别讨论两种MD情况下GK和GKA之间的氢键随时间的变化：（1）GK和GKA单配体作用.从图4（A）中可以看到其中的两个氢键非常稳定，只有第三个氢键不太稳定，在MD运行到3ns附近发生断裂，在6ns左右又重新形成；（2）GK和GKA、葡萄糖双配体作用.分析三个氢键距离随模拟时间的变化，可以发现在整个MD过程中，三个氢键都非常稳定.由于以上两种情况下MD模拟参数都完全相同，唯一的区别在于第二种情况下含有葡萄糖小分子.因此我们可以得出这样的结论，即葡萄糖的存在使GK和GKA的结合更加稳定.这和文献报道［3］一致，即GKA可以提高GK对葡萄糖的活性，同时葡萄糖的存在，反过来又可以促进GK和GKA的结合.迄今为止，还没有发现在没有葡萄糖存在的条件下，GK和GKA单配体结合的晶体结构，这也说明GK单独和GKA 在一起是不能长期稳定存在的.2.3 GKA和GK之间包结自由能的分析通过MM－PBSA方法［11］计算了两个模拟体系（GK－GKA和GK－GKA－glucose）中GKA和GK之间的相互作用能，结果如表1所示.分析计算结果，发现不论是哪一种模拟体系，GK和GKA之间的主要作用力都是范德华力.尽管GK 和GKA之间真空下的静电作用能为负值，但是极性的溶剂化能却是一个更大的正值，导致它们之间总的静电作用能为正值，不利于GK和GKA的稳定结合.通过能量计算表明，GKA主要是通过范德华力和GK结合在一起.另外分析GK－GKA和GK－GKA－glucose两个体系的能量差别，发现葡萄糖的存在，使得GK和GKA之间的能量变得更低.具体分析各个能量项，发现主要的差别不是范德华力，而是静电作用力.葡萄糖的存在使得GK－GKA之间的静电作用能变低，最终导致总的包结能降低.以上的计算结果和文献报道［3］一致，即葡萄糖可以促进GK和GKA的稳定结合.2.4 GKA对GK构象影响的分析前面我们分析了GK和GKA之间的主要作用力.下面将分析GKA对GK构象的影响.首先分析GKA对GK夹角的影响（如图5所示）.从图5中可以看到对于GK－free，大小域的夹角在开始的1ns变化很大，从15°迅速增加到30°左右.随后，GK－free夹角的值趋于稳定，说明GK达到一种相对较稳定的构象.GK－free夹角的这种变化趋势和其RMSD的变化趋势相一致，从图3中可以看到，其RMSD 值在开始的1 ns上升非常快，迅速达到3nm后趋于稳定.和GK－free相比，加了GKA以后，GK大小域的夹角和初始结构夹角相比，变化较小.从图5中可以看到，在整个5ns时间内，夹角的值比较稳定，在18°左右波动.因此可以得出GKA 对GK的构象变化确实起到了一种门闩的作用，约束着GK大小域的开合.但是5ns 以后，我们发现GKGKA体系的夹角从18°增加到22°.以上结果表明GKA可以在一段时间内约束GK的构象，使之处于活化状态，但不能一直保持GK的构象不变，这和以往的实验结果是一致的［3］.为了进一步研究GKA对GK构象的影响，分别用动态相关性矩阵（Dynamic cross－correlation matrices，DCCM）［18－19］分析了四种MD情况下GK 448个氨基酸残基内部的相关性，如图6所示.DCCM最吸引人的地方在于预测一些非相邻的氨基酸之间的动态相关性，也就是我们通常所说的蝴蝶效应，这是静态的晶体结构无法提供的.因此DCCM已被广泛用于蛋白质动态相关性分析［18－19］，并取得了很好的预测效果.图6 （A）为GK－free的DCCM 图，从图6（A）可以看到，大部分氨基酸残基的相关性很小，颜色以浅蓝色和绿色为主，很小的黄色区域表示该区域的氨基酸残基之间存在正的相关性.另外还存在一些很小的蓝色区域，表示该部分残基之间的负相关性较强.和GK－free相比，GK－GKA的氨基酸残基负相关性非常明显，部分区域的颜色显示为深蓝色（如图6（B）所示）.仔细分析，我们发现GK 300～400的氨基酸残基和100～150的氨基酸残基之间有很明显的负相关性.进一步分析晶体结构（如图7所示），发现300～400的氨基酸残基位于GK的大域，而100～150的氨基酸残基位于GK的小域，它们之间相距甚远.在没有GKA存在的情况下，它们之间并没有负相关性运动，但是在GKA存在的情况下，这种负相关运动则非常明显.表明GKA可以影响GK大小域的开合运动.这和我们前面分析GKA对GK夹角的影响是一致的.GK－glucose的DCCM如图6（C）所示.分析图6（C）发现，在葡萄糖存在的情况下，GK的正负相关性运动都非常不明显，可忽略不计.我们推测是因为GK和葡萄糖的强烈静电作用，导致GK的构象非常稳定，所以观察不到GK的相关性运动.最后一个DCCM图用来说明GK－GKA－glucose的氨基酸残基的相关性运动（见图6（D））.分析图6（D）发现，GK－GKA－glucose的负相关性介于GK －GKA和GK－glucose之间.我们推测是因为葡萄糖和GK的相互作用，导致GKA对GK的负相关性影响减小.以上DCCM分析和作者以前所研究的Y214C［11］和M197V［12］活性突变的DCCM结果基本相同.表明该活化剂的作用类似于GK本身的一些活性突变，都是通过约束GK的构象来提高GK的活性.通过分子动力学模拟分析了GKA和GK的相互作用.包结自由能分析表明，范德华作用是GK和GKA体系的主要驱动力.同时葡萄糖的存在使得GK和GKA结合得更加稳定.构象分析表明，GKA可以在一段时间约束GK的构象，但是不能一直维持GK的活化构象不变，这和以往的实验结果相一致.本研究从原子水平解释了GKA的活化机制，为今后设计新的GKA小分子提供了一定的理论参考.致谢：南京大学理论与计算化学研究所提供计算机服务器用于计算.【相关文献】［1］FERRÉT，RIU E，BOSCH F，et al.Evidence from transgenic mice that glucokinase is rate limiting for glucose utilization in the liver［J］.FASEB J，1996，10（10）：1213－1218.［2］NOGUEIRA F N，SANTOS M F，NICOLAU J.Influence of streptozotocin－induced diabetes on hexokinase activity of rat salivary glands［J］.J Physiol Biochem，2005，61（3）：421－428.［3］KAMATA K，MITSUYA M，NISHIMURA T，et al.Structural basis for allosteric regulation of themonomeric allosteric enzyme human glucokinase［J］.Structure，2004，12（3）：429－438.［4］VIONNET N，STOFFEL M，TAKEDA J，et al.Nonsense mutation in the glucokinase gene causes early－onset non－insulin－dependent diabetes mellitus［J］.Nature，1992，356（6371）：721－722.［5］GRIMSBY J，SARABU R，CORBETT W L，et al.Allosteric activators of glucokinase：potential role in diabetes therapy［J］.Science，2003，301（5631）：370－373.［6］FUTAMURA M，HOSAKA H，KADOTANI A，et al.An allosteric activator of glucokinase impairs the interaction of glucokinase and glucokinase regulatory protein and regulates glucose metabolism［J］.J Biol Chem，2006，281（49）：37668－37674.［7］IINO T，TSUKAHARA D，KAMATA K，et al.Discovery of potent and orally active 3－alkoxy－5－phenoxy－N－thiazolyl benzamides as novel allosteric glucokinase activators ［J］.Bioorg Med Chem，2009，17（7）：2733－2743.［8］NISHIMURA T，IINO T，MITSUYA M，et al.Identification of novel and potent 2－amino benzamide derivatives as allosteric glucokinase activators［J］.J Bioorg Med Chem Lett，2009，19（5）：1357－1360.［9］MATSCHINSKY F M.Assessing the potential of glucokinase activators in diabetes therapy：nature reviews drug discovery［J］.Nat Rev Drug Discov，2009，8（5）：399－416.［10］ZHANG Jian，LI Chen Jing，CHEN Kai Xian，et al.Conformational transition pathway in the allosteric process of human glucokinase［J］.Proc Natl Acad Sci USA，2006，103（36）：13368－13373.［11］姚雪霞.葡萄糖激酶Y214C活性突变的理论研究［J］.化学研究，2010，21（6）：71－76.［12］姚雪霞，徐素明.葡萄糖激酶 M197V突变的活化机理研究［J］.化学学报，2011，69（4）：405－410.［13］WANG Jun Mei，MORIN P，WANG Wei，et e of MM－PBSA in reproducing the binding free Energies to HIV－1 RT of TIBO derivatives and predicting the binding mode to HIV－1RT of efavirenz by docking and MM－PBSA ［J］.J Am Chem Soc，2001，123（22）：5221－5230.［14］GOUDA H，KUNTZ I D，CASE D A，et al.Free energy calculations for theophylline binding to an RNA aptamer：MMPBSA and comparison of thcrmodynamic integration methods［J］.Biopolymers，2003，68（1）：16－34.［15］CASE D A，DARDEN T A，CHEATHAM T E，et al.AMBER 10［CP］.University of California：San Francisco，2008.［16］FRISCH M J，TRUCKS G W，SCHLEGEL H B，et al.Gaussian 03，Revision C.02［CP］.Gaussian，Inc.：Wallingford CT，2004.［17］DARDEN T，YORK D，PEDERSEN L.Particle mesh Ewald：An N.log（N）method for Ewald sums in large systems［J］.J Chem Phy，1993，98（12）：10089－10092.［18］LI Li Wei，UVERSKY V N，DUNKER A K，et al.A computational investigation of allostery in the catabolite activator protein［J］.J Am Chem Soc，2007，129（50）：15668－15676.［19］GHOSH A，VISHVESHWARA S.A study of communication pathways in methionylt RNA synthetase by molecular dynamics simulations and structure network analysis ［J］.Proc Natl Acad Sci USA，2007，104（40）：15711－15716.。

gap 反应力场分子动力学-概述说明以及解释1.引言1.1 概述在当代材料科学和化学研究中，理解和预测材料性质是至关重要的。

然而，由于原子尺度的复杂性和大量的可能性，实验方法难以满足对材料行为的全面分析和预测的需求。

因此，计算机模拟成为研究者们实现这一目标的有力工具之一。

GAP 反应力场分子动力学是一种计算模拟方法，结合了反应力场和分子动力学模拟的优势。

反应力场是一种理论模型，它通过描述原子间相互作用来模拟并预测材料的结构和性质。

而分子动力学模拟则是一种仿真方法，通过数值求解牛顿运动方程来模拟原子和分子在时间和空间上的运动。

GAP 反应力场分子动力学模拟方法的主要优势之一是其高保真度和高效率。

通过基于大量实验数据训练得到的反应力场，该方法可以准确地描述各种材料的结构和性质。

同时，分子动力学模拟的高效率使得研究者们可以在相对较短的时间内获得大量的数据，从而更全面地理解材料的行为。

在本文中，我们将首先介绍GAP 反应力场的基本原理和方法。

然后，我们将详细解释分子动力学模拟的原理，并阐述GAP 反应力场分子动力学方法的应用场景和方法。

接下来，我们将展示一些具体的研究结果，并进行相应的讨论。

最后，我们将总结本文的研究成果，并展望未来该方法在实际应用中的潜力和价值。

通过本文的阅读，读者将能够全面了解GAP 反应力场分子动力学方法的基本原理和应用，以及其在材料科学和化学研究中的实际应用价值。

同时，本文也将为读者提供一些启示，鼓励他们在自己的研究中尝试并发展这一方法，以推动材料科学和化学领域的进步。

文章结构部分主要用于介绍整个文章的组织架构，让读者清楚地了解文章的内容安排。

本文的结构如下：1. 引言1.1 概述1.2 文章结构1.3 目的1.4 总结2. 正文2.1 GAP 反应力场介绍2.2 分子动力学模拟原理2.3 应用场景与方法2.4 结果与讨论3. 结论3.1 研究总结3.2 展望未来3.3 实际应用价值3.4 结束语在引言部分，我们将概括性地介绍本篇文章的内容，并提供文章结构的概览，同时明确文章的目的和总结。

amber中生成小分子模板第一步：生成小分子模板蛋白质中各氨基酸残基的力参数是预先存在的，但是很多模拟过程会涉及配体分子，这些有机小分子有很高的多样性，他们的力参数和静电信息不可能预存在库文件中，需要根据需要自己计算生成模板。

amber中的antechamber 程序就是生成小分子模板的。

生成模板要进行量子化学计算，这一步可以由antechamber中附带的mopac完成，也可以由gaussian完成，这里介绍用gaussian计算的过程。

建议在计算前用sybyl软件将小分子预先优化，不要用gaussian优化，大基组从头计算进行几何优化花费的计算时间太长。

gaussian计算的输入文件可以用antechamber程序直接生成，生成后去掉其中关于几何优化的参数即可将小分子优化后的结构存储为mol2各式，上传到工作目录，用antechamber程序生成gaussian输入文件，命令如下：antechamber -i 49.mol2 -fi mol2 -o 49.in -fo gzmat这样可以生成49.in文件，下载到windows环境，运行gaussian计算这个文件，如果发现计算时间过长或者内存不足计算中断，可以修改文件选择小一些的基组。

获得输出文件49.out之后将它上传到工作目录，再用antechamber生成模板，命令如下：antechamber -i 49.out -fi gout -o 49mod.mol2 -fo mol2 -c resp运行之后就会生成一个新的mol2文件，如果用看图软件打开这个文件会发现，原子的颜色很怪异，这是因为mol2的原子名称不是标准的原子名称，看图软件无法识别。

下面一步是检查参数，因为可能会有一些特殊的参数在gaff中不存在需要程序注入，命令如下：parmchk -i 49mod.mol2 -f mol2 -o 49mod这样那些特殊的参数就存在49mod这个文件中了第二步：处理蛋白质文件amber自带的leap程序是处理蛋白质文件的，他可以读入PDB各式的蛋白质文件，根据已有的力场模板为蛋白质赋予键参数和静电参数。

第33卷 第2期原 子 与 分 子 物 理 学 报Vol.33No.22016年4月JOURNAL OF ATOMIC AND MOLECULAR PHYSICSApr.2016doi:103969/j.issn.1000-0364．2016．04．027抑制剂8CA 与脂肪细胞脂肪酸结合蛋白(A -FABP )结合模式的分子动力学研究尹妍妍1,梁志强1,王 伟1,伊长虹1,李洪云1,赵 娟1,张庆刚2(1．山东交通学院理学院,济南250357;2．山东师范大学物理与电子科学学院,济南250014)摘 要:脂肪细胞脂肪酸结合蛋白A -FABP (Adipocyte fatty -acid binding protein )是治疗脂质调节生物过程相关疾病的重要靶标．采用分子动力学模拟和MM -PBSA 方法研究抑制剂8CA 与A -FABP 结合模式,结果表明静电相互作用和范德华作用驱动了抑制剂8CA 与A -FABP 的结合．基于残基的能量分解表明抑制剂8CA 与R126间的极性相互作用为抑制剂与A -FABP 的结合提供了重要贡献,该残基与8CA 的相互作用较好地稳定了抑制剂与A -FABP 复合物的稳定性．期望该研究可为治疗炎症㊁动脉硬化和代谢病药物设计提供一定的理论指导．关键词:分子动力学模拟;MM -PBSA 方法;A -FABP ;抑制剂-残基相互作用中图分类号:Q641．12 文献标识码:A 文章编号:1000-0364(2016)02-0339-06Insight into binding mode of inhibitor 8CA to A -FABPbased on molecular dynamics simulationYIN Yan-Yan 1,LIANG Zhi-Qiang 1,WANG Wei 1,YI Chang-Hong 1,LI Hong-Yun 1,ZHAO Juan 1,ZHANG Qing-Gang 2(1．School of Science,Shandong Jiaotong University,Jinan 250357,China;2．College of Physics and Electronics,Shandong Normal University,Jinan 250014,China)Abstract :Adipocyte fatty -acid binding protein (A -FABP)is an important target of drug designs treating some diseases related to lipid -mediated biology．Molecular dynamics (MD)simulations coupled with molecularmechanics Poisson -Boltzmann surface area (MM -PBSA)calculation were carried out to study the binding mode of 8CA to A -FABP．The results show that electrostatics and van der Waals interactions drive the binding of 8CA to A -FABP．The calculation from residue -based free energy decomposition suggests that the polar in-teraction of 8CA with the residue R126provides an important contribution to the 8CA binding．This polar interac-tion plays a key role in the stabilization of 8CA /A -FABP complex．We expect that this study can contribute some theoretical guidance for design of potent inhibitors within the fields of metabolic disease,inflammation andatherosclerosis．Key words :Molecular dynamics simulation;MM -PBSA;A -FABP;Inhibitor -residue interaction收稿日期:2015-09-06基金项目:国家自然科学基金(11504206,11274206);山东省高等学校科技计划项目(J14LJ07)作者简介:尹妍妍(1980 ),女,山东高唐人,讲师,生物分子的模拟计算.E-mail:yinyanyansdu@126．com1 引 言脂肪酸结合蛋白(FABP)是小型的细胞质蛋白,存在于多种组织中,被定义为脂肪酸转运蛋白．该蛋白周围包绕了大量的相关蛋白,一些除结合脂肪酸外,还结合了疏水的配体．最近几年,原子与分子物理学报第33卷对FABP的组织分布㊁配体亲和力和特异性以及其结构特性进行了集中研究,结果均表明FABPs 参与细胞内脂质代谢[1,2]．众所周知,脂肪细胞脂肪酸结合蛋白(A-FABP)是9个已知脂肪酸结合蛋白亚型之一,该蛋白多出现于脂肪组织和巨噬细胞中[3,4]．先前研究表明A-FABP在诸如代谢综合征和心脏病的特定专一性方面起到重要作用[5]．研究也显示A -FABP的移除和其功能紊乱能减少心脏病和糖尿病的患病几率[6]．因此A-FABP成为治疗与脂质调节生物相关疾病药物设计的重要靶标．A-FABP功能的药物干扰的确在糖尿病和心脏病的治疗中起到重要作用[7]．Ringom等人设计几个不同结构的抑制剂,这些抑制剂能较好抑制A-FABP的活性功能[7,8]．A-FABP抑制剂的问世向倍受心脏病和糖尿病煎熬的病人显现光明．但是目前抑制剂与A-FABP的结合模式仍有待于进一步阐明．故原子层次上研究抑制剂与A-FABP的作用机制将有助于高效的A-FABP抑制剂的研发．分子动力学模拟和自由能计算已成为研究抑制剂与蛋白㊁蛋白与蛋白相互作用的重要工具[9-17]．为达到研究目的,选取由Ringom研发的抑制剂8CA,研究其与A-FABP的作用机制．A -FABP及抑制剂8CA的结构显示在图1中．抑制剂8CA与A-FABP的结合能力达到纳摩尔数量级[7]．本文采用基于分子动力学轨迹的MM-PBSA方法[18-20]计算抑制剂8CA与A-FABP的结合自由能,从原子层次上探索8CA与A-FABP 的作用机制．2 模型和计算2．1 分子动力学模拟蛋白质库的晶体结构(3FR2)被提取且用作动力学模拟的初始结构．通过从3FR2中直接删除8CA获得未结合的A-FABP结构．晶体结构中所有结晶水分子保留在初始模型中,晶体结构中缺失的氢原子由Amber[21]中的leap模块添加．A-FABP和水分子的力场参数由Amber中的FF03[22]力场产生．小分子抑制剂8CA的力场参数由GAFF力场产生[23],且BCC电荷分配给抑制剂的每一个原子．复合体8CA/A-FAB溶解在由显性的TIP3P水分子组成的水盒子里,水盒子边缘与复合体最近原子的距离是12．0Å．两个钠离子添图1 分子结构:(A)8CA/A-FABP复合物的结构;(B)抑制剂8CA的结构Fig.1 Structures of molecules:(A)for the structureof8CA/A-FABP complex;(B)for thestructure of8CA加到水和复合体8CA/A-FAB组成的系统,以确保整个系统呈电中性．为了消除原子间一些不合理的接触,采用Amber中的Sander工具执行两步的系统优化以便对整个复合体系统进行弛豫,获得合理的初始优化结构:首先约束溶质,优化溶剂和中和离子,采用的约束力常数为100kcal/(mol㊃Å2);然后无约束地弛豫整个系统,使整个系统沿着低能量结构的方向优化．每一步优化先执行2000步的最陡下降优化,接着执行5000步的共轭梯度优化．随后在500ps内把系统从0K缓慢地加热到300 K,随后在常温300K,常压1标准大气压下进行500ps的动力学平衡．最后是两个30ns的无约束分子动力学模拟,一个是8CA/A-FABP复合体的模拟,另一个为未结合A-FABP的模拟．模拟期间采用SHAKE方法约束与氢原子相连化学键的伸缩[24]．模拟积分步长设定为2fs．PME方法用来计算长程的静电相互作用．应用周期性边界条件以消除溶剂盒子的边缘效应[25]．非成键相互作用的截断值为12．0Å．同时,监测和评估主链原子相对于晶体结构的RMSD随时间的函数,以判断系统平衡是否可靠．043第2期尹妍妍,等:抑制剂8CA与脂肪细胞脂肪酸结合蛋白(A-FABP)结合模式的分子动力学研究2．2 MM-PBSA方法计算结合自由能采用MM-PBSA方法计算抑制剂8CA与A-FABP的结合自由能．从动力学模拟轨迹中每隔一定的时间间隔取出8CA/A-FABP复合体的结构,删掉构象中的水分子和钠离子,使用下面方程计算结合自由能．ΔG=ΔE MM+ΔG sol-TΔS(1)该方程中ΔE MM表示气相中的分子力学能,可以分解为两部分:ΔE MM=ΔE ele+ΔE vdw(2)式中ΔE ele和ΔE vdw分别表示气相中的8CA与A-FABP静电相互作用和范德华能．方程(1)中ΔG sol 表示溶解自由能对8CA与A-FABP结合的贡献,也可分解成如下两个成分:ΔG sol=ΔG pb+ΔG surf(3)式中ΔG pb和ΔG surf分别代表极性溶解自由能和非极性溶解自由能对抑制剂8CA与A-FABP结合的贡献．前者使用有限差分法求解泊松-玻尔兹曼方程获得结果,计算中溶质和溶剂的电介常数分别设定为1．0和80．0．后者使用如下经验方程求解:ΔG surf=γSASA+β(4)式中γ和β的值设定为0．00542kcal/(mol㊃Å2)和0．92kcal/mol．方程(1)中最后一项TΔS是与8CA结合诱导的自由度变化形成的熵变对结合自由能的贡献,该项使用简振模和传统的热力学进行计算．2．3 抑制剂-残基相互作用为详细地评估每个残基对结合自由能的贡献,基于残基的能量分解法被采用计算抑制剂与A-FABP中每一个残基的相互作用能．所采用的计算方程表示如下:ΔG=ΔE ele+ΔE vdw+ΔG pb+ΔG surf(5)方程中ΔE ele和ΔE vdw分别表示气相中静电相互作用和范德华作用．ΔG pb和ΔG surf分别表示极性溶解自由能和非极性溶解自由能．这四项通常采用GB 模型进行计算．3 结果和讨论3．1 分子动力学平衡的稳定性为了评估MD模拟的质量,使用Amber中的工具Ptraj程序计算分析了A-FABP主链原子在动力学模拟过程中相对于晶体结构的均方根偏差(RMSD)随时间的变化情况(图2)．图2表明8CA/A-FABP复合体在模拟进行5ns后,系统达到了稳定平衡,而未结合的A-FABP在模拟进行大约8ns才趋近平衡．由图2观察到RMSD的平均值为1．21Å,涨落范围低于0．55Å,该结果表明系统的平衡是可靠的,能为系统的后加工分析进行可靠的构象取样．此外由图2还观察到复合体的RMSD的平均值要比未结合A-FABP的RMSD的平均值低大约0．45Å,这个结果表明抑制剂8CA的结合对A-FABP产生一定的约束作用,与8CA的结合有利于稳定A-FABP的结构．图2 MD模拟中A-FABP主链原子的均方根偏差Fig.2 Root-mean-square deviation of backbone atoms in A-FABP during molecular dynamics simula-tion3．2 结构的柔性为考察抑制剂8CA的结合对A-FABP结构柔性的影响,A-FABP中各残基Cα原子涨落(RMSF)显示在图3中．图3表明8CA结合对A -FABP内部的柔性产生重要影响．未结合的A-FABP具有更大的柔性．在残基30㊁45㊁57和106附近区域的RMSF有了明显下降,而在残基90和100附近RMSF的值轻微升高．该结果表明抑制剂8CA可能与该区域的残基产生相互作用,进而影响了A-FABP内部结构的变化．3．3 结合自由能计算该工作采用基于单轨迹的MM-PBSA方法评估了抑制剂8CA与A-FABP的结合自由能．分离的相互作用成分有助于抑制剂与A-FABP结合的驱动力和结合模式的阐明．各成分及其结合自由能均呈现在表1中．由表1观察到:抑制剂8CA与A-FABP的结合自由能是-11．89kcal/mol．这表明抑制剂143原 子 与 分 子 物 理 学 报第33卷表1 采用MM -PBSA 方法计算所得到的自由能成分(kcal /mol)aTable 1 Results of free energies (in kcal /mol)calculated by MM -PBSA methodComponent ComplexA -FABP8CADeltaMeanStdMeanStdMean StdMeanStdE ele -3806．6171．42-3684．0571．273．990．99-126．557．91E vdw -511．6217．76-488．4317．2110．361．98-33．552．72E int 2718．9333．452660．2332．4558．704．550．000．00E gas -1599．3071．74-1512．2571．2673．054．36-160．107．72G surf 41．340．7342．440．733．760．02-4．860．08G pb 1782．8862．00-1816．4562．71-79．230．98112．815．93G sol -1741．5462．10-1774．0262．36-75．470．97107．955．93G pbele -5589．4825．67-5500．5025．60-75．250．88-13．745．56G pbtot -3340．8334．02-3286．2633．62-2．424．15-52．155．31-TS2252．515．042225．426．4070．350．5640．265．78ΔG bind-11．89aE ele 和E vdw 分别表示静电作用能和范德华作用能,E int 表示由键伸㊁键角弯折和二面角扭转贡献的内能,E gas =E ele +E vdw +E int ;G surf和G pb 分别是非极性溶剂化能和极性溶剂化能,G sol 是溶解自由能且G sol =G surf +G pb ;G pbele =E ele +G pb ;G pbtot =G sol +E gas ;TS 自由度变化对自由能的贡献图3 MD 模拟中残基Cα原子的原子涨落Fig.3 Root -mean -square fluctuation of Cαatoms in A-FABP during molecular dynamics simulation8CA 与A -FABP 的结合能力较强．8CA 与A -FABP 间的范德华作用为-33．5kcal /mol,该相互作用为8CA 的结合提供有利因素．8CA 与A -AFABP 的静电相互作用非常强,达到-126．55kcal /mol,在强度上是范德华作用的4倍．但是这个有利的因素被不利的极性溶解自由能抵消一部分,二者共产生-13．74kcal /mol 的有利因素．与溶剂可及表面积相对应的非极性溶解自由能为8CA 结合提供了较小有利贡献,其值为-4．86kcal /mol．综合比较而言,非极性溶解自由能远弱于范德华相互作用和极性相互作用,因此范德华作用和静电相互作用驱动抑制剂8CA 与A -FABP 的结合．此外,表一还表明自由度变化诱导的熵变对8CA 结合的贡献是40．26kcal /mol．该结果表明A -FABP 在与8CA 结合后,运动的自由度数减少,反映了8CA 结合对A -FABP 的运动产生一定的约束作用．当然该相互作用严重消弱了8CA 与A -FABP 的结合．3．4 抑制剂-残基相互作用为定量的理解每一个残基对8CA 结合的贡献及阐明8CA /A -FABP 复合体的结构和亲和能关系,采用基于残基的能量分解方法计算8CA 与A -FABP 中每一个残基的相互作用,其相互作用谱显示在图4中．基于动力学模拟获得最小能量结构在图5中给出了抑制剂8CA 相对于重要残基的相对位置．依据图4,抑制剂8CA 与A -FABP 中9个残基产生较强的相互作用．这9个残基包括F16㊁M20㊁S53㊁S55㊁F57㊁A75㊁R78㊁R106和R126．这9个残基与8CA 的作用共分为3类:F16㊁M20㊁S53㊁S55㊁F57和A75与8CA 产生疏水性相互作用,R126与8CA 形成氢键相互作用,R78㊁R106和R126与8CA 的带电基团产生强烈的极性相互作用．243第2期尹妍妍,等:抑制剂8CA 与脂肪细胞脂肪酸结合蛋白(A -FABP )结合模式的分子动力学研究图4 抑制剂8CA 与A -FABP 中各个分离残基的相互作用谱Fig.4 Inhibitor -residue interaction spectrum betweenthe inhibitor 8CA and A -FABP图5 A -FABP 和8CA 复合体中主要残基的相对位置及形成氢键Fig.5 Relative geometries of the key residues in 8CA /A -FABP complex and hydrogen bondsR126与8CA 的作用最强,其作用强度为-6．8kcal /mol．图5表明该相互作用主要来自两方面的贡献:(1)R126与8CA 形成两个氢键相互作用;(2)R126带电基团与8CA 的羰基产生强烈的相互作用．同理R78和R106与8CA 的相互作用能分别为-2．96和-2．66kcal /mol．这两个相互作用也是源自R78和R106的带电基团与8CA 羰基的静电相互作用．残基F16和M20与8CA 的相互作用能分别是-1．45和-1．26kcal /mol,从结构上看(图5),这两个相互作用可能源自F16和M20疏水性侧链与8CA 疏水性基团间的π-π相互作用和CH -π相互作用．残基S53和S55的CH 基团临近8CA 的疏水性苯环(图5),有利于与8CA 产生CH -π相互作用,这两个相互作用分别为8CA 的结合提供了-1．3和-1．61kcal /mol 的能量贡献(图4)．图5表明F57的苯环和A75的烷基位于8CA 疏水性环的附近,分别与8CA 产生π-π相互作用和CH -π相互作用,该相互作用的强度分别是-1．84和-1．85kcal /mol．本工作的研究吻合了先前的实验和理论计算研究[7,8]．基于上述分析8CA 与残基间的氢键作用㊁静电相互作用以及π-π相互作用和CH -π相互作用驱动了8CA 与A -FABP 的结合．在以A -FABP 为靶标的药物设计过程中,如何调节小分子结构,优化这些相互作用成为成功研发高效A -FABP 抑制剂的重要一环．4 结 论本文采用分子动力学模拟和MM -PBSA 方法研究了抑制剂8CA 与A -FABP 结合模式．动力学分析披露8CA 的结合约束了A -FABP 内部一些区域的运动,降低了该区域的原子涨落,稳定了复合物结构．结合自由能结果证明8CA 能与A -FABP 产生较强相互作用,且极性相互作用和范德华作用驱动8CA 与A -FABP 的结合．基于残基的能量分解计算表明8CA 能与A -FABP 的9个残基产生较强相互作用,且氢键作用㊁静电相互作用以及π-π相互作用和CH -π相互作用是抑制剂结合的主导力量．我们期望这个研究能为以A -FABP 为靶标的药物设计提供合理的理论指导．参考文献:[1] Hotamisligil G S．Inflammation and metabolic disorders [J]．Nature ,2006,444:860．[2]Banaszak L,Winter N,Xu Z,et al ．Lipid -bindingproteins:a family of fatty acid and retinoid transport proteins [J]．Adv ．Protein Chem ．,1994,45:89．[3] Hertzel A V,Bernlohr D A．The mammalian fatty acid -binding protein multigene family:molecular and ge-netic insights into function [J]．Trends Endocrinol Metab ,2000,11:175．[4]Hotamisligil G S,Johnson R S,Distel R J,et al ．Un-coupling of obesity from insulin resistance through a targeted mutation in aP2,the adipocyte fatty acid bind-ing protein [J]．Science ,1996,274:1377．[5] Uysal K T,Scheja L,Wiesbrock S M,et al ．Improved glucose and lipid metabolism in genetically obese mice343原子与分子物理学报第33卷lacking aP2[J]．Endocrinology,2000,141:3388．[6] Makowski L,Boord J B,Maeda K,et al．Lack ofmacrophage fatty-acid–binding protein aP2protectsmice deficient in apolipoprotein E against atherosclero-sis[J]．Nat．Med．,2001,7:699．[7] Ringom R,Axen E,Uppenberg J,et al．Substitutedbenzylamino-6-(trifluoromethyl)pyrimidin-4(1H)-ones:a novel class of selective human A-FABP inhibitors[J]．Bioorg．Med．Chem．Lett．,2004,14:4449．[8] Chen J,Wang J,Zhu W．Binding modes of three in-hibitors8CA,F8A and I4A to A-FABP studied basedon molecular dynamics simulation[J]．PloS one,2014,9:e99862．[9] Cheng W Y,Liang Z Q,Zhang Q G,et al．Insight in-to p53-MDM2interaction based on molecular dynam-ics simulation and molecular mechanics[J]．J．At．Mol．Phys．,2012,29:393(in Chinese)[程伟渊,梁志强,张庆刚,等.p53-MDM2相互作用的分子力学和动力学研究[J]．原子与分子物理学报,2012,29:393][10] Yi C H,Liang Z Q,Wang W,et al．Insight into inter-action mode of HIV-1protease with inhibitor BEGbased molecular docking method[J]．J．At．Mol．Phys．,2014,31:718(in Chinese)[伊长虹,梁志强,王伟,等.HIV-1蛋白酶与抑制剂BEG相互作用的分子对接计算研究[J]．原子与分子物理学报,2014,31:718][11] Chen J,Wang J,Zhu W,et al．A computational anal-ysis of binding modes and conformation changes ofMDM2induced by p53and inhibitor bindings[J]．J．Comput．Aided Mol．Des．,2013,27:965．[12] Wu E L,Han K L,Zhang J Z．Selectivity of neutral/weakly basic P1group inhibitors of thrombin and tryp-sin by a molecular dynamics study[J]．Chem．Eur．J．,2008,14:8704．[13] Wang J,Yang H,Zuo Z,et al．Molecular dynamicssimulations on the mechanism of transporting methyla-mine and ammonia by ammonium transporter amtB[J].J．Phys．Chem．B,2010,114:15172．[14] Zhu T,Xiao X,Ji C,et al．A new quantum calibratedforce field for zinc-protein complex[J]．J．Chem．Theory Comput．,2013,9:1788．[15] Shi S H,Zhang S L,Zhang Q G．Binding free energyinsight into interaction mechanism of inhibitor PDIQwith MDM2[J]．J．At．Mol．Phys．,2014,31:618(in Chinese)[时术华,张少龙,张庆刚.MDM2与抑制剂PDIQ作用机制的结合自由能研究[J]．原子与分子物理学报,2014,31:618] [16] Wu E L,Mei Y,Han K L,et al．Quantum and mo-lecular dynamics study for binding of macrocyclic in-hibitors to human[alpha]-thrombin[J]．Biophys．J．,2007,92:4244．[17] Yi C H,Liang Z Q,Wang W,et al．A computationalinsight into interaction of inhibitor3G3with HIV-1prease based on molecular dynamics simulation[J]．J．At．Mol．Phys．,2013,30:994(in Chinese)[伊长虹,梁志强,王伟,等.HIV-1蛋白酶与抑制剂3G3相互作用机制的分子动力学研究[J]．原子与分子物理学报,2013,30:994][18] Wang J,Morin P,Wang W,et al．Use of MM-PBSAin reproducing the binding free energies to HIV-1RTof TIBO derivatives and predicting the binding mode toHIV-1RT of efavirenz by docking and MM-PBSA[J]．Journal of the American Chemical Society,2001,123:5221．[19] Chen J,Liang Z,Wang W,et al．Revealing origin ofdecrease in potency of darunavir and amprenavir a-gainst HIV-2relative to HIV-1protease by molecu-lar dynamics simulations[J]．Sci．Rep．,2014,4:6872．[20] Hu G,Wang J．Ligand selectivity of estrogen receptorsby a molecular dynamics study[J]．Eur．J．Med．Chem．,2014,74:726．[21] Case D A,Darden T A,Cheatham III T E,et al．AM-BER12,San Francisco:University of California,2012．[22] Duan Y,Wu C,Chowdhury S,et al．A point-chargeforce field for molecular mechanics simulations of pro-teins based on condensed-phase quantum mechanicalcalculations[J]．J．Comp．Chem．,2003,24:1999．[23] Wang J,Wolf R M,Caldwell J W,et al．Developmentand testing of a general amber force field[J]．J．Comput．Chem．,2004,25:1157．[24] Coleman T G,Mesick H C,Darby R L．Numerical in-tegration[J]．Ann．Biomed．Eng．,1977,5:322．[25] Darden T,York D,Pedersen L．Particle mesh Ewald:An N㊃log(N)method for Ewald sums in large sys-tems[J]．J．Chem．Phys．,1993,98:10089．443。