CN111755064A - 基于cmap势函数的耦合二面角参数优化方法及蛋白质力场 - Google Patents

基于cmap势函数的耦合二面角参数优化方法及蛋白质力场 Download PDF

Info

Publication number
CN111755064A
CN111755064A CN202010599142.XA CN202010599142A CN111755064A CN 111755064 A CN111755064 A CN 111755064A CN 202010599142 A CN202010599142 A CN 202010599142A CN 111755064 A CN111755064 A CN 111755064A
Authority
CN
China
Prior art keywords
cmap
force field
distribution
phi
psi
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Pending
Application number
CN202010599142.XA
Other languages
English (en)
Inventor
康微
蒋帆
吴云东
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Peking University Shenzhen Graduate School
Original Assignee
Peking University Shenzhen Graduate School
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Peking University Shenzhen Graduate School filed Critical Peking University Shenzhen Graduate School
Priority to CN202010599142.XA priority Critical patent/CN111755064A/zh
Publication of CN111755064A publication Critical patent/CN111755064A/zh
Pending legal-status Critical Current

Links

Images

Classifications

    • GPHYSICS
    • G16INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
    • G16BBIOINFORMATICS, i.e. INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR GENETIC OR PROTEIN-RELATED DATA PROCESSING IN COMPUTATIONAL MOLECULAR BIOLOGY
    • G16B5/00ICT specially adapted for modelling or simulations in systems biology, e.g. gene-regulatory networks, protein interaction networks or metabolic networks
    • GPHYSICS
    • G16INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
    • G16BBIOINFORMATICS, i.e. INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR GENETIC OR PROTEIN-RELATED DATA PROCESSING IN COMPUTATIONAL MOLECULAR BIOLOGY
    • G16B15/00ICT specially adapted for analysing two-dimensional or three-dimensional molecular structures, e.g. structural or functional relations or structure alignment
    • G16B15/20Protein or domain folding
    • GPHYSICS
    • G16INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
    • G16BBIOINFORMATICS, i.e. INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR GENETIC OR PROTEIN-RELATED DATA PROCESSING IN COMPUTATIONAL MOLECULAR BIOLOGY
    • G16B35/00ICT specially adapted for in silico combinatorial libraries of nucleic acids, proteins or peptides
    • G16B35/20Screening of libraries

Landscapes

  • Physics & Mathematics (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Engineering & Computer Science (AREA)
  • Health & Medical Sciences (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Spectroscopy & Molecular Physics (AREA)
  • Evolutionary Biology (AREA)
  • Biophysics (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Biotechnology (AREA)
  • Theoretical Computer Science (AREA)
  • General Health & Medical Sciences (AREA)
  • Medical Informatics (AREA)
  • Molecular Biology (AREA)
  • Chemical & Material Sciences (AREA)
  • Library & Information Science (AREA)
  • Physiology (AREA)
  • Crystallography & Structural Chemistry (AREA)
  • Biochemistry (AREA)
  • Peptides Or Proteins (AREA)

Abstract

本发明涉及分子动力学模拟技术领域,具体涉及蛋白质力场,更具体涉及一种基于CMAP势函数的耦合二面角参数优化方法。该优化方法以现有的蛋白质力场为母力场,以卷曲库中的局部固有构象偏好为拟合目标,对母力场中氨基酸残基的二面角参数进行优化,使得分子动力学模拟获得的局部构象分布与卷曲库中的分布达到预定的相似度。本发明还涉及通过该优化方法获得的蛋白质力场。本发明可以简单而高效地对蛋白质分子中耦合的二面角进行同时优化,二面角参数优化过程简单、高效、可自动化,获得了兼容性强的残基特异性力场RSFF2C,可广泛用于多种MD软件。

Description

基于CMAP势函数的耦合二面角参数优化方法及蛋白质力场
技术领域
本发明涉及分子动力学模拟技术领域,具体涉及蛋白质力场,更具体涉及一种基于CMAP势函数的耦合二面角参数优化方法及所得的蛋白质力场。
背景技术
作为细胞中生命活动的主要承担者,蛋白质一直是生命科学及其相关学科中的重要研究对象。关于蛋白质的研究主要围绕“序列-结构-功能”这个基本范式开展,其中蛋白质的结构尤其是关注重点,因为蛋白质的空间结构能帮助人们更好地从微观上理解其功能,进而为蛋白质工程和药物设计与筛选等研究提供重要的基础。目前,得益于基因测序技术的发展和普及,蛋白质序列数据库快速增长,然而,已知实验结构的蛋白质相对比较少。尽管蛋白质结构测定实验技术(包括X射线晶体衍射和核磁共振技术等)已经逐渐发展成熟,通过实验手段来确定蛋白质的空间结构仍然是一个成本高且周期长的复杂过程。作为一种辅助手段,通过计算方法来研究蛋白质的空间结构、蛋白质的构象转化、以及蛋白质-蛋白质相互作用等成为了一个很重要的研究手段。
在众多计算方法中,分子动力学(Molecular Dynamics,MD)模拟是一个很重要的蛋白质研究工具。MD模拟的基本思想是把物质看成由遵循经典力学的粒子构成,通过求解牛顿运动方程而得到体系中各粒子的运动轨迹从而模拟系统随时间推进的微观过程,然后通过统计分析得到系统的结构、热力学性质、动力学性质等。自1959年由Alder和Wainwright提出以来,经过几十年的蓬勃发展,MD模拟如今已成为物理、化学、生物和材料等多学科的重要研究工具。尤其是在生物大分子的研究中,MD模拟可以帮助人们在很广的时空尺度进行研究,深入理解生物大分子的结构和功能。一个成功的MD模拟依赖于准确的力场和充分的构象采样。用于描述粒子间相互作用的势能函数形式和具体参数被称为力场。近年来,计算能力的不断提升和多种增强抽样方法的发展大大增强了构象采样的效率,同时也对力场的准确度提出了更高的要求。
蛋白质力场的发展始于20世纪80年代初,在过去的三十多年间,研究者们发展出了多个系列的蛋白质力场,主要包括AMBER、CHARMM、GROMOS和OPLS等几个系列。在势能函数形式方面,各个系列经典的蛋白质力场的函数表达大同小异,形式相对简单,主要由键伸缩项Ebond、键角弯折项Eangle、二面角扭曲项Etorsion(包括用来保证原子共平面的impropertorsion)、Lennard-Jones(L-J)ELJ势和库伦势ECoulomb,如图1和公式1所示。前两项由谐振子模型来近似描述,二面角扭曲项通常由一系列傅里叶展开式来描述,L-J势用来描述由三个或三个以上的键相隔的电中性原子对之间的能量(色散作用与交换排斥作用之和)随距离的变化,通常称为“范德华”项,库伦势则用来描述带电粒子对之间的库伦吸引或者排斥能,ELJ和ECoulomb又合称为非键相互作用项。此外,在CHARMM系列力场中还增加一个Urey-Bradley项作为对标准的键伸缩项和键角弯折项的补充,以更好地重现振动光谱实验的结果。
E=Ebond+Eangle+Etorsion+ELI+ECoulomb 式(1)
纵观各系列力场的发展历史,可发现两个突出的特点:其一,近年来大部分力场参数改进的方向集中在骨架二面角或侧链二面角,尤其是骨架二面角的参数,由于其直接影响了各种二级结构间的平衡;而且大部分力场的优化策略是对各个二面角进行单独优化,基本不考虑相邻二面角之间的耦合,直到2003年,Brooks等人首次向CHARMM力场中引入2D格点形式的修正函数(Correction MAP,CMAP)对骨架二面角
Figure BDA0002558010610000021
ψ进行修正,以改善CHARMM22力场存在的π螺旋含量过高的问题。其二,二面角参数主要通过拟合小分子模型(最常用的是二肽模型)在气相中的QM计算得到的
Figure BDA0002558010610000022
ψ势能面获得。
大部分蛋白质力场的参数化依赖于小分子模型(最常用的是二肽模型)在气相中的QM计算结果。近年来,随着实验手段的进步和实验数据的积累,利用实验数据对力场参数进一步优化成为一种新趋势。其中,利用PDB结构中卷曲区域(coil region)的残基(即卷曲库,coil library)的构象统计结果来优化力场取得了不错的成效。所谓卷曲库,就是从PDB中选择高分辨率的蛋白质晶体结构,除去位于各种二级结构中的残基,只留下卷曲区域,这些残基便构成卷曲库。通过对卷曲库进行统计分析,人们发现不同氨基酸残基具有不同的局部构象分布(即二面角分布)特征。由于这些局部构象分布实际反映了残基的固有特征,因此又被称为固有构象偏好。参照这些统计结果,可以对现有蛋白质力场进行进行残基特异性的优化,使得每种残基都可以被准确描述,这样的力场被称为残基特异性力场(Residue-Specific Force Field,RSFF)。
此前,蒋帆和周晨阳等人先后分别对OPLS-AA/L力场和AMBER 99SB力场进行了优化得到RSFF1力场(Jiang,F.;Zhou,C.Y.;Wu,Y.D.,Residue-specific force field basedon the protein coil library.RSFF1:modification of OPLS-AA/L.J.Phys.Chem.B2014,118(25),6983-98)和RSFF2力场(Zhou,C.Y.;Jiang,F.;Wu,Y.D.,Residue-specificforce field based on protein coil library.RSFF2:modification of AMBERff99SB.J.Phys.Chem.B 2015,119(3),1035-47)。他们主要优化了二面角参数,但发现仅用传统的傅里叶展开式无法准确描述相邻二面角之间的耦合,因此引入了一些特殊的1-4/1-5/1-6相互作用(如图2所示)以更好地重现卷曲库中各种氨基酸残基的固有构象分布。经测试,RSFF1和RSFF2力场比其它力场能更好地重现氨基酸二肽的NMR 3JHNHα耦合常数。另外,RSFF1和RSFF2力场是分别针对TIP4P-Ew和TIP3P水模型进行的参数化,因此,它们也通常搭配对应的水模型一起使用。然而,由于TIP3P和TIP4P-Ew这些水模型低估了水分子与蛋白质的短程色散作用,RSFF1和RSFF2力场与这些水模型搭配模拟固有无序蛋白质(Intrinsically Disordered Proteins,IDP)时,常给出过于塌缩的构象系综。将它们与更准确的水模型TIP4P-D搭配使用时,能给出更合理的构象系综,但这种搭配低估A14和AQ15等体系的螺旋度。为了兼顾折叠态和无序态蛋白,吴昊南等人对RSFF2力场进行了修正,发展了新的RSFF2+力场(Wu,H.-N.;Jiang,F.;Wu,Y.-D.,Significantly improved proteinfolding thermodynamics using a dispersion-corrected water model and a newresidue-specific force field.J.Phys.Chem.Lett.2017,8(14),3199-3205)。RSFF2+力场在蛋白质分子骨架Oi和Hi+4原子之间加入了微弱的吸引作用,用来补偿经典力场在α螺旋结构中缺失的多体效应(电荷极化)。
这些RSFF力场后来被成功用于蛋白质从头折叠模拟、蛋白质结构优化、环肽结构预测、IDP的系综研究等工作,且通常比相应的母力场表现更佳,这说明以卷曲库的统计结果为参照的二面角参数优化策略是可行的。但是这些力场存在以下两方面的缺点:
(1)参数优化复杂,不是一个通用的好办法。由于残基的固有构象分布同时受到二面角参数和局部的1-4/1-5/1-6作用等多方面的影响,参数优化往往需要经验丰富的研究者反复尝试对多种作用进行精细调节,会耗费大量的人力资源和计算资源。目前的RSFF系列力场是针对特定的显式水模型优化得到的,如要将它们与其它水模型搭配使用,原则上需要对力场参数重新优化。另外,如果要对其它力场进行类似的二面角优化或者需要更改拟合目标时,也需要重复此过程。重新优化绝不是一个简单的重复过程,而是要将前述的耗时耗力的过程每个环节再重新做一遍。因此,我们需要一种更简单、高效、最好是可自动化的参数优化方法。
(2)目前RSFF系列力场的软件兼容性差,无法被用于多种MD软件。由于包含局部的1-4/1-5/1-6项,除了GROMACS软件外,其它的MD软件并不支持或很难实现RSFF系列力场,这使得其它软件使用者不能使用这些力场,从而无法实现对力场的广泛验证、应用以及进一步发展。实际上,不同MD软件各具优点:GROMACS的CPU计算效率很高,这得益于它的代码实现结合了硬件特点,且该软件免费;AMBER软件更新迭代快,包含很多模拟和分析的工具,使用方便,是一个主流的MD软件;NAMD支持大规模并行计算;OpenMM灵活性很大,既可以独立使用,也可以作为第三方函数库被调用;LAMMPS支持大规模原子/分子并行模拟计算等。因此,理想的情况是采用一种更普适的参数优化方法,使得优化得到的力场能被广泛用于多个MD软件,充分利用各MD软件的优势。
发明内容
针对现有RSFF力场存在的参数化复杂、兼容性差的问题,受CHARMM力场启发,本发明提出一种基于CMAP势函数的耦合二面角参数优化方法,进而获得了优化的蛋白质力场。
因此,在第一方面,本发明提供一种基于CMAP势函数的耦合二面角参数优化方法,该优化方法以现有的蛋白质力场为母力场,以卷曲库中的局部固有构象偏好为拟合目标,对母力场中氨基酸残基的二面角参数进行优化,使得分子动力学模拟获得的局部构象分布与卷曲库中的分布达到预定的相似度。
进一步地,对于Ala和Gly残基,以卷曲库中的局部构象(φ,ψ)分布作为拟合目标,得到最终的骨架CMAP。
进一步地,对于除Ala和Gly残基以外的其它氨基酸残基,均使用Ala的骨架CMAP项,并以卷曲库中的3D(χ1,φ,ψ)分布作为拟合目标,得到残基特异性的侧链CMAP,并参考卷曲库中的相应分布对χi(i>1)二面角叠加了傅里叶展开式修正项。
进一步地,对Ala/Gly二肽进行REMD模拟,获得二面角(φ,ψ)分布,并与卷曲库中的(φ,ψ)分布进行相似度比较,如果所得的相似度没有达到预定的相似度,则更新(φ,ψ)CMAP进行迭代优化,直到达到预定的相似度,得到最终的骨架CMAP。
进一步地,对除Ala/Gly二肽以外的其他二肽进行REMD模拟,获得二面角(φ,ψ,χ1)分布,并与卷曲库中的(φ,ψ,χ1)分布进行相似度比较,如果所得的相似度没有达到预定的相似度,则分解两个(φ,ψ,χ1)分布的差异,更新(χ1,φ)CMAP和更新更新(χ1,ψ)CMAP进行迭代优化,直到达到预定的相似度,得到最终的侧链CMAP。
还进一步地,在与卷曲库中的(φ,ψ,χ1)分布进行相似度比较中,还比较三种χ1旋转异构体g+/t/g-的比例是否合适。
在本发明的一个具体实施方案中,预定的相似度大于0.99。
在本发明的一个具体实施方案中,母力场为ff14SB力场,并采用TIP3P水模型。
在第二方面,本发明涉及一种蛋白质力场,该蛋白质力场通过根据本发明第一方面所述的优化方法获得。
在本发明的一个具体实施方案中,蛋白质力场为残基特异性力场RSFF2C。
本发明的有益效果:
(1)简单而高效地对蛋白质分子中耦合的二面角进行同时优化,二面角参数优化过程简单、高效、可自动化。由于没有使用局部L-J势来描述相邻二面角的耦合,而是使用更简单的2D CMAP势函数,参数优化过程中将不需要对参数进行手动调整和验证,可以大大简化参数优化流程,实现简单、高效的参数优化,甚至可实现完全自动化。在此基础上,可搭配不同的水模型,使用不同的拟合目标,对很多现有的力场进行快速的二面角参数优化,大大丰富现有蛋白质力场家族,为蛋白质模拟提供更多的备选力场。
(2)获得兼容性强的残基特异性力场RSFF2C,可广泛用于多种MD软件。目前,多种MD软件包括GROMACS、AMBER、CHARMM、OpenMM等均支持CMAP势函数,因此这种基于CMAP势函数的残基特异性力场具有很好的软件兼容性。另外,这种参数优化方法底层的自由能分解策略具有普遍的物理基础,因此也可被用于其它分子力场的优化。
附图说明
图1显示了经典蛋白质力场的能量项示意图,从左至右依次为键伸缩项、键角弯折项、二面角扭曲项、L-J势和库伦势;
图2显示了RSFF1/2力场中的特殊的局部1-4/1-5/1-6相互作用示例;
图3显示了本发明的参数优化方法所采用的自由能分解策略;
图4显示了原始的卷曲库中的局部构象分布与重构的分布非常相似,证明本发明的参数优化方法所采用的自由能分解策略是合适的;
图5显示了与TIP3P水模型搭配的RSFF2C力场的CMAP参数优化总体流程;
图6显示了RSFF2C力场优化过程中骨架二面角(φ,ψ)分布随优化周期的变化,目标分布(卷曲库中的分布)与RSFF2力场下的模拟轨迹中的(φ,ψ)分布也示于图中方便对比,其中(A)对应Ala的(φ,ψ)分布,而(B)为天冬氨酸残基(简称Asp)在三种χ1旋转异构体下骨架二面角(φ,ψ)的分布;
图7显示了使用不同力场模拟所得的局部构象分布与卷曲库中的分布的相似度,其中(A)骨架二面角(φ,ψ)分布的相似度,包括总体(φ,ψ)分布以及不同χ1旋转异构体下的(φ,ψ)分布,(B)三种χ1旋转异构体含量的对比,横坐标为其在卷曲库中的含量(+代表g+,-代表g-,菱形代表t);
图8显示了从二肽(无Pro和His)的REMD模拟的300K副本轨迹中计算得到的3JHNHα耦合常数与实验值之间的相关性;
图9显示了RSFF2C力场中最终优化得到的CMAP修正项,图中相邻等高线能量差为0.5kcal/mol;
图10显示了Glu、Gln和Asn的侧链二面角χ2和χ3的PMF优化前(即ff14SB+CMAPs,叉号标记的曲线)后(即RSFF2C,三角标记的曲线)与卷曲库(无标记的曲线)对比图。
具体实施方式
为了使本发明所解决的技术问题、所采用的技术方案及所获得的有益效果更加清楚明白,以下结合附图及具体实施例对本发明进行进一步详细说明。应当理解,此处所描述的具体实施例仅仅用于解释本发明,并不意在限定本发明。
本发明的基于CMAP势函数的耦合二面角参数优化方法采用一种自由能分解策略(参见图3)。如图3所示,在除了Gly和Ala之外的氨基酸残基中,侧链原子团R与骨架两个相邻肽键平面之间存在相互作用,这些侧链-骨架相互作用直接决定了氨基酸残基的局部固有构象偏好,表现在不同的氨基酸残基具有不同的(φ,ψ,χ1)分布特征。从能量上来看(图3上),在除丙氨酸残基(简称Ala)和甘氨酸残基(简称Gly)之外的其它氨基酸残基中,其侧链原子团R与两个相邻的肽平面(椭圆形)的相互作用可以被分解为两部分:R与左侧肽平面的相互作用和与右侧肽平面的相互作用。从局部构象分布上来看(图3下),氨基酸残基的3D(φ,ψ,χ1)分布可以表示为三个2D分布的乘积:(φ,ψ)分布及(χ1,φ)分布和(χ1,ψ)分布。2D分布图的z坐标为频数,且相邻等高线的自由能量差为RT(其中R为理想气体常数,T为温度)。如无特殊说明,本发明中所有的二面角2D分布图均采用此梯度。
首先假定R对骨架的作用可以分解为两部分:R与左侧肽平面的相互作用和与右侧肽平面的相互作用,其中R与骨架作用强度主要由侧链二面角x1以及骨架二面角φ,ψ决定(图3上)。按照这个思路,这个氨基酸残基X的自由能GX(x1,φ,ψ)可以分解为三个量(公式2):Ala的自由能GA(φ,ψ),两个与x1相关的分量GX(x1,φ)和GX1,ψ)。根据波尔兹曼关系式,公式2可以改写为公式3,即在局部构象分布上来看,三维的x1,φ,ψ分布可以分解为三个二维的分布的乘积:第一项是Ala的φ,ψ分布,后两项分别为χ1,φ分布和χ1,ψ分布。在已知pX1,φ,ψ)和pA(φ,ψ)的情况下,后两项可以通过公式4多次迭代求得,其中pX1,φ,ψ)和pA(φ,ψ)可以通过高斯核密度估计方法从蛋白质卷曲库或MD模拟轨迹中提取。
GX(x1,φ,ψ)=GA(φ,ψ)+GX(x1,φ)+GX(x1,ψ) 式(2)
pX1,φ,ψ)=pA(φ,ψ)·pX(x1,φ)·pX(x1,ψ) 式(3)
Figure BDA0002558010610000071
Figure BDA0002558010610000072
为了证明以上提出的自由能分解策略是合适的,对于每一个氨基酸残基,在对其卷曲库分布进行如上分解后,又从获得的三个二维分布出发重构了三维分布。通过比较发现,原始的卷曲库中的局部构象分布与重构的分布非常相似(参见图4)。图4A显示了自由能分解策略的示意图,来自卷曲库的初始的3D分布(original 3D distr.from coil lib.)可以分解为三个2D分布,而且可以从所得的2D分布重构出一个与原来分布非常像的3D分布(reconstructed 3D distr.),以谷氨酸残基(简称Glu)为例,S=0.996,相邻等高线的自由能量差为RT。图4B显示了除Ala和Gly外的所有氨基酸残基的原始卷曲库分布与重构得到的分布之间存在非常高的相似度,其中左侧为不同旋转异构体(g+、t和g-)下(φ,ψ)分布的相似度,右侧为三种旋转异构体的含量的比较。
这个测试说明了以上提出的分解策略并不会损失重要的信息,其对信息的保留度足够高。值得注意的是,这里的能量分解是对自由能的分解而不是经典力场框架下的势能的成对分解,因为这里的能量已经包含了溶剂化作用:蛋白质晶体结构中的卷曲部分是蛋白质晶体中高度溶剂化的区域,从MD模拟轨迹中提取得到的自由能面同样也包含了溶剂的影响。
在验证了能量分解公式2后,可分别对来自卷曲库的3D自由能面与来自MD的3D自由能面进行分解,而两者的差也可以分解成三个2D项,每项分别用格点形式的CMAP势函数来表示,对母力场施加这些CMAP项可弥补该自由能差,使MD自由能与卷曲库自由能相等(如以下式5、6和7所示),从而实现对耦合二面角的同时优化。
GX,coil1,φ,ψ)=GA,coil(φ,ψ)+GX,coil1,φ)+GX,coil1,ψ) 式5
GX,MD1,φ,ψ)=GA,MD(φ,ψ)+GX,MD1,φ)+GX,MD1,ψ) 式6
GX,coil1,φ,ψ)-GX,MD1,φ,ψ)=EA,CMAP(φ,ψ)+EX,CMAP1,φ)+EX,CMAP1,ψ) 式7
理论上来说,本发明所采用的二面角参数优化策略可被用于多种蛋白质分子力场的优化,而不仅仅限于经典分子力场。以下仅以经典力场为例,来阐述该参数优化流程。
以AMBER ff14SB为母力场,搭配TIP3P水模型,以卷曲库中的局部固有构象偏好为拟合目标,对母力场中每一个残基的二面角参数进行优化,使得MD模拟获得的局部构象分布与卷曲库中的分布达到预定的相似度。具体拟合目标为:对于Ala和Gly而言,拟合Coil-6卷曲库中的局部构象(φ,ψ)分布;对于其它残基,使用Coil-3+卷曲库中的3D(χ1,φ,ψ)分布以及其它χi(i>1)分布作为拟合目标。该过程优化得到的力场被命名为RSFF2C。
1.参数优化流程:
图5显示了与TIP3P水模型搭配的RSFF2C力场的CMAP参数优化总体流程。图5左边是Ala和Gly的骨架二面角(φ,ψ)CMAP的优化。当Ala的(φ,ψ)CMAP优化完成后,对其它氨基酸残基的骨架二面角使用Ala的(φ,ψ)CMAP,在此基础上优化其它氨基酸残基的侧链二面角χ1相关的CMAP,包括(χ1,φ)CMAP和(χ1,ψ)CMAP,如图5右边所示。
(1)骨架CMAP项的优化流程
对于Ala和Gly,因为没有侧链基团的存在,骨架CMAP项的参数优化相对简单(如图5左边所示)。以Ala为例,首先使用REMD模拟对二肽(Ace-Ala-Nme)进行充分采样,获得300K时的二面角分布pMD(φ,ψ)。基于玻尔兹曼分布,卷曲库中的分布pcoil(φ,ψ)和模拟所得的分布pMD(φ,ψ)的差异可以通过对现有力场(AMBER ff14SB)增加一个CMAP势函数来弥补,即Ala骨架二面角的能量修正项ECMAP(φ,ψ)为:
Figure BDA0002558010610000081
其中pcoil(φ,ψ)和pMD(φ,ψ)分别为卷曲库和MD轨迹中Ala骨架二面角处于φ,ψ的概率,通过高斯核密度估计获得。这种参数优化策略与之前RSFF1/2的优化策略存在显著不同,虽然均以卷曲库为拟合目标,但RSFF1/2是通过用一维的傅里叶展开去拟合各个二面角的平均力势(Potential of Mean Force,PMF),并同时对多个局部1-4/1-5/1-6作用进行经验性调整而获得的。
实际上,ECMAP(φ,ψ)的优化过程是迭代的。首先,使ECMAP(φ,ψ)≡0,即没有对母力场进行任何修改,然后在每一次优化循环中,通过公式9a计算目标自由能面Gcoil(φ,ψ)与当前模拟所得的自由能面
Figure BDA0002558010610000094
之间的差
Figure BDA0002558010610000095
然后通过公式9b将这项与上一轮优化所得的修正能量项
Figure BDA0002558010610000096
相加得到新的修正项
Figure BDA0002558010610000097
Figure BDA0002558010610000098
Figure BDA0002558010610000099
循环进行该优化过程直到卷曲库分布与模拟所得分布的相似度系数S>0.99为止。
(2)侧链CMAP项的优化流程
对于非Gly非Ala氨基酸残基,其拟合过程如图6右边所示。在以上的自由能分解策略的基础上,可以将残基X在卷曲库中的自由能与模拟中的自由能的差分解为三项:
GX,coil(x1,φ,ψ)-GX,ff14SB1,φ,ψ)
=[GA,coil(φ,ψ)+GX,coil(x1,φ)+GX,coil(x1,ψ)]
-[GA,MD(φ,ψ)+GX,MD(x1,φ)+GX,MD(x1,ψ)]
=EA,CMAP(φ,ψ)+EX,CMAP(x1,φ)+EX,CMAP(x1,ψ) 式(10)
当对X的骨架二面角施加与Ala同样的修正EA,CMAP(φ,ψ)后,目标自由能与模拟自由能差变为:
Figure BDA00025580106100000910
基于波尔兹曼关系式,公式11变为:
Figure BDA0002558010610000091
与公式4类似,δpX,CMAP1,φ)和δpX,CMAP1,ψ)可通过以下两个方程从已知的pX,coil1,φ,ψ)和pX,ff14SB+CMAP_A1,φ,ψ)获得:
Figure BDA0002558010610000092
Figure BDA0002558010610000093
首先使δpX,CMAP1,ψ)≡1,然后通过循环迭代公式13,直到δpX,CMAP1,ψ)和δpX,CMAP1,φ)均收敛(一般在10次迭代以内可以实现)。然后δpX,CMAP1,φ)和δpX,CMAP1,ψ)可以分别被转化为CMAP势修正项:
EX,CMAP(x1,φ)=-RTlnδpX,CMAP(x1,φ) 式(14a)
EX,CMAP(x1,ψ)=-RTlnδpX,CMAP(x1,ψ) 式(14b)
与Ala和Gly的骨架CMAP ECMAP(φ,ψ)的优化类似,最终的EX,CMAP1,φ)和EX,CMAP1,ψ)也通过几个优化循环获得。
以上只对(φ,ψ,χ1)分布进行优化,对于χi(i>1),原有的AMBER ff14SB力场对大部分氨基酸残基的χi(i>1)描述相对准确。对于少数氨基酸残基,ff14SB模拟得到的侧链二面角χ2和χ3的PMF与卷曲库中的PMF存在较大差别,本发明采用经典的傅里叶展开式拟合两者的差异,如公式15,与之前的RSFF1/2类似。
Figure BDA0002558010610000101
2.参数优化结果
CMAP势显著提高了模拟中的局部构象分布与卷曲库中的分布的相似度。图6给出了参数优化过程中拉氏图的变化情况,模拟所得的(φ,ψ)分布经过两轮参数优化便可与卷曲库分布很像。具体地,如图7所示,在使用母力场AMBER ff14SB时,模拟得到的分布已经与卷曲库很像(S>0.9),跟之前的AMBER ff99SB的结果(对于Ala,S=0.855;对于Gly,S=0.923)相比,ff14SB中Ala的骨架参数有了改进,在骨架二面角描述方面较之前的版本更准确。然而,仔细对比AMBER ff14SB和卷曲库的拉氏图,存在明显的差异:(1)对于Gly二肽,AMBER ff14SB给出的αR和PPII区域大致呈圆形,而卷曲库呈现倾斜的椭圆形。(2)对于Ala二肽的αR区域,也存在类似的现象,而且在φ=-160°时,卷曲库中Ala更倾向于处于ψ=160°(β区域)而不是在ψ=0°(α区域),而ff14SB下的模拟结果并没有体现这个倾向。这些问题源自母力场无法准确地描述φ/ψ耦合,无法通过继续优化原有的二面角参数来解决(母力场只有一维的傅里叶展开,没有显式描述二面角耦合的函数项)。这里使用的CMAP修正项很好地修正了母力场在这方面的不足,在经过一轮优化后,S就能提高至0.99,第二轮后S达到0.995以上,第二轮获得的CMAP修正项就是RSFF2C力场的最终参数。另外,与原来的RSFF2的结果(S=0.97)相比,RSFF2C的结果与卷曲库分布的相似度更高,说明使用CMAP修正项与RSFF2中所用的特殊1-4/1-5/1-6L-J势相比具有明显优势。即,不仅CMAP项优化流程简单,而且对二面角耦合的描述更准确。
对于除了Ala和Gly外的其它残基,使用侧链CMAP可以改善总体(包含三个χ1旋转异构体)和不同χ1旋转异构体下的(φ,ψ)图(详见图7的(A))。在RSFF2C中,18个氨基酸残基中有14个残基的总体相似度S>0.990,最低S=0.985,而ff14SB所得的S值约为0.894,最低为0.782(Asn),最高为0.979(Pro)。RSFF2的二肽模拟结果的平均S为0.975,大大低于RSFF2C(0.993)。特别是,在RSFF2模拟结果中,Asp在旋转异构体t下以及Leu在旋转异构体g-下给出的S值约为0.80,这两者在RSFF2C中已提高到0.98以上。
当分别考虑每个旋转异构体时,g+旋转异构体下,ff14SB模拟所得的分布的S均大于0.80,而Asp和Asn的t旋转异构体的S值低至<0.6,另外有7个氨基酸残基的g-旋转异构体的S值<0.80。这说明,在母力场中,使用同一组骨架二面角参数可能不足以同时准确描述同一残基的不同旋转异构体的(φ,ψ)分布特征。另一方面,RSFF2C对所有旋转异构体均保持较高的S值,优于RSFF2。以天冬氨酸残基Asp为例,优化过程中不同旋转异构体下的(φ,ψ)分布以及卷曲库分布和相应的RSFF2分布如图6的(B)所示。各力场下的模拟所得的分布之间最明显的差异在于旋转异构体t。在旋转异构体t下,ff14SB和RSFF2均偏爱PPII构象,但卷曲库中的统计分布明显偏向于C7构象。
RSFF2C中的侧链相关的CMAP项可以解决RSFF2中的这个问题。在RSFF2C模拟结果中,仅Val和Ile处于g-时(φ,ψ)分布的S<0.97(图7的(A))。这表明尽管不同旋转异构体下的(φ,ψ)分布存在显着差异,但RSFF2C中的侧链CMAP仍能准确捕获与旋转异构体有关的自由能分布。对于Val和Ile,由于每个Cβ原子上都有两个Cγ基团,这些侧链具有β分支的残基可能不像其它残基一样适合之前提出的自由能分解方案(图4)。尽管如此,与RSFF2相比,RSFF2C下的模拟分布仍然与卷曲库分布高度一致。
另外,RSFF2C和RSFF2给出的三个侧链χ1旋转异构体(g+、t和g-)的相对偏好与卷曲库中的相应分布非常相似(图7的(B))。而ff14SB给出的结果与卷曲库统计显著不同,平均而言,t旋转异构体倾向于占有更大的比例。
此外,使用CMAP修正项后,二肽的3JHNHα耦合常数与实验值符合更好,如图8所示。对比ff14SB和RSFF2C的结果,模拟结果与实验值之间的相关系数从0.59提升为0.93。RSFF2C的结果与之前RSFF2的结果类似。
3.RSFF2C的最终参数
综上,RSFF2C力场中的所有CMAP修正项包括:(1)分别对Ala和Gly的骨架(φ,ψ)CMAP修正项;(2)其它氨基酸残基除了共用Ala的骨架(φ,ψ)CMAP项,还有与侧链二面角χ1相关的(χ1,φ)CMAP修正项和(χ1,ψ)CMAP修正项。具体如图9所示。
除了CMAP项之外,也对一些残基的其它侧链二面角χi(i>1)进行了修正,主要是Glu、Gln和Asn这三个残基。以Glu为例,从ff14SB+CMAPs模拟结果得到的χ2的PMF曲线上290°处能量最低,而180°处最高,这与卷曲库统计结果相差较大;而χ3的PMF曲线上则存在多个极小值点,也与统计结果不符。对这些二面角叠加了傅里叶展开形式的修正。图10给出了修正前后这些二面角的PMF与卷曲库结果的对比。
以上应用了具体实例对本发明进行了阐述,只是用于帮助理解本发明,并不用以限制本发明。本发明所属技术领域的技术人员依据本发明的构思,还可以做出若干简单推演、变形或替换。这些推演、变形或替换方案也落入本发明的权利要求范围内。

Claims (10)

1.一种基于CMAP势函数的耦合二面角参数优化方法,其特征在于,以现有的蛋白质力场为母力场,以卷曲库中的局部固有构象偏好为拟合目标,对母力场中氨基酸残基的二面角参数进行优化,使得分子动力学模拟获得的局部构象分布与卷曲库中的分布达到预定的相似度。
2.根据权利要求1所述的优化方法,其特征在于,对于Ala和Gly残基,以卷曲库中的局部构象(φ,ψ)分布作为拟合目标,得到骨架CMAP项。
3.根据权利要求1所述的优化方法,其特征在于,对于除Ala和Gly残基以外的其它氨基酸残基,均使用Ala的骨架CMAP项,并以卷曲库中的3D(χ1,φ,ψ)分布作为拟合目标,得到残基特异性的侧链CMAP项,并参考卷曲库中的相应分布对χi(i>1)二面角叠加了傅里叶展开式修正项。
4.根据权利要求2所述的优化方法,其特征在于,对Ala/Gly二肽进行REMD模拟,获得二面角(φ,ψ)分布,并与卷曲库中的(φ,ψ)分布进行相似度比较,如果所得的相似度没有达到所述预定的相似度,则更新(φ,ψ)CMAP进行迭代优化,直到达到所述预定的相似度,得到最终的骨架CMAP。
5.根据权利要求3所述的优化方法,其特征在于,对除Ala/Gly二肽以外的其他二肽进行REMD模拟,获得二面角(φ,ψ,χ1)分布,并与卷曲库中的(φ,ψ,χ1)分布进行相似度比较,如果所得的相似度没有达到所述预定的相似度,则分解两个(φ,ψ,χ1)分布的差异,更新(χ1,φ)CMAP和(χ1,ψ)CMAP进行迭代优化,直到达到所述预定的相似度,得到最终的侧链CMAP。
6.根据权利要求5所述的优化方法,其特征在于,在与卷曲库中的(φ,ψ,χ1)分布进行相似度比较中,还比较三种χ1旋转异构体g+/t/g-的比例是否合适。
7.根据权利要求1-6中任一项所述的优化方法,其特征在于,所述预定的相似度大于0.99。
8.根据权利要求1-6中任一项所述的优化方法,其特征在于,所述母力场为ff14SB力场,并采用TIP3P水模型。
9.一种蛋白质力场,其特征在于,所述蛋白质力场通过根据权利要求1-8中任一项所述的优化方法获得。
10.根据权利要求9所述的蛋白质力场,其特征在于,所述蛋白质力场为残基特异性力场RSFF2C。
CN202010599142.XA 2020-06-28 2020-06-28 基于cmap势函数的耦合二面角参数优化方法及蛋白质力场 Pending CN111755064A (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202010599142.XA CN111755064A (zh) 2020-06-28 2020-06-28 基于cmap势函数的耦合二面角参数优化方法及蛋白质力场

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202010599142.XA CN111755064A (zh) 2020-06-28 2020-06-28 基于cmap势函数的耦合二面角参数优化方法及蛋白质力场

Publications (1)

Publication Number Publication Date
CN111755064A true CN111755064A (zh) 2020-10-09

Family

ID=72677836

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202010599142.XA Pending CN111755064A (zh) 2020-06-28 2020-06-28 基于cmap势函数的耦合二面角参数优化方法及蛋白质力场

Country Status (1)

Country Link
CN (1) CN111755064A (zh)

Cited By (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112216350A (zh) * 2020-11-05 2021-01-12 深圳晶泰科技有限公司 物理严格且相空间重叠最大化的相对自由能计算方法
CN112233733A (zh) * 2020-11-05 2021-01-15 深圳晶泰科技有限公司 分子力场质量控制系统及其控制方法
CN114300035A (zh) * 2021-12-21 2022-04-08 上海交通大学 一种用于蛋白质力场模拟的个性化参数生成方法
WO2022094870A1 (zh) * 2020-11-05 2022-05-12 深圳晶泰科技有限公司 物理严格且相空间重叠最大化的相对自由能计算方法
CN116453587A (zh) * 2023-06-15 2023-07-18 之江实验室 一种任务执行方法、装置、存储介质及电子设备

Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN105787292A (zh) * 2014-12-18 2016-07-20 中国科学院大连化学物理研究所 蛋白质折叠的并行预测方法
CN107203702A (zh) * 2017-01-13 2017-09-26 北京理工大学 一种分析蛋白质侧链构象含时动力学演化的方法

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN105787292A (zh) * 2014-12-18 2016-07-20 中国科学院大连化学物理研究所 蛋白质折叠的并行预测方法
CN107203702A (zh) * 2017-01-13 2017-09-26 北京理工大学 一种分析蛋白质侧链构象含时动力学演化的方法

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
FAN JIANG 等: "Developments and Applications of Coil-Library-Based Residue-specific Force Fields for Molecular Dynamics Simulations of Peptides and Proteins", J. CHEM. THEORY COMPUT., pages 2761 - 2773 *
WEI KANG 等: "Universal Implementation of a Residue-Specific Force Field Based on CMAP Potentials and Free Energy Decomposition", J. CHEM. THEORY COMPUT., pages 4474 - 4486 *

Cited By (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112216350A (zh) * 2020-11-05 2021-01-12 深圳晶泰科技有限公司 物理严格且相空间重叠最大化的相对自由能计算方法
CN112233733A (zh) * 2020-11-05 2021-01-15 深圳晶泰科技有限公司 分子力场质量控制系统及其控制方法
WO2022094870A1 (zh) * 2020-11-05 2022-05-12 深圳晶泰科技有限公司 物理严格且相空间重叠最大化的相对自由能计算方法
CN112216350B (zh) * 2020-11-05 2022-09-13 深圳晶泰科技有限公司 物理严格且相空间重叠最大化的相对自由能计算方法
CN112233733B (zh) * 2020-11-05 2023-04-07 深圳晶泰科技有限公司 分子力场质量控制系统及其控制方法
CN114300035A (zh) * 2021-12-21 2022-04-08 上海交通大学 一种用于蛋白质力场模拟的个性化参数生成方法
CN116453587A (zh) * 2023-06-15 2023-07-18 之江实验室 一种任务执行方法、装置、存储介质及电子设备
CN116453587B (zh) * 2023-06-15 2023-08-29 之江实验室 一种基于分子动力学模型预测配体亲和力的任务执行方法

Similar Documents

Publication Publication Date Title
CN111755064A (zh) 基于cmap势函数的耦合二面角参数优化方法及蛋白质力场
Lesage et al. Smoothed biasing forces yield unbiased free energies with the extended-system adaptive biasing force method
Song et al. AIMOES: Archive information assisted multi-objective evolutionary strategy for ab initio protein structure prediction
Dawson et al. Coarse-grained modeling of RNA 3D structure
Cardamone et al. Multipolar electrostatics
Tom et al. Genarris 2.0: A random structure generator for molecular crystals
Klenin et al. Modelling proteins: Conformational sampling and reconstruction of folding kinetics
US20210383897A1 (en) Structure search method, structure search apparatus, and non-transitory computer-readable storage medium for storing structure search program
Kim et al. Constructing an interpolated potential energy surface of a large molecule: A case study with bacteriochlorophyll a model in the Fenna–Matthews–Olson complex
Héry et al. X-ray diffuse scattering and rigid-body motion in crystalline lysozyme probed by molecular dynamics simulation
Peng et al. Simulating large-scale conformational changes of proteins by accelerating collective motions obtained from principal component analysis
Chen et al. ELF: an extended-Lagrangian free energy calculation module for multiple molecular dynamics engines
Simonson Free energy calculations
Skolnick et al. Computational studies of protein folding
Marinelli et al. Force-correction analysis method for derivation of multidimensional free-energy landscapes from adaptively biased replica simulations
Rahimi et al. Comparison of on-the-fly probability enhanced sampling and parallel tempering combined with metadynamics for atomistic simulations of RNA tetraloop folding
Miao et al. Assessing the Performance of Peptide Force Fields for Modeling the Solution Structural Ensembles of Cyclic Peptides
Rick et al. An implementation of replica exchange with dynamical scaling for efficient large-scale simulations
Brunger et al. Annealing in crystallography: a powerful optimization tool
Maruyama et al. Analysis of molecular dynamics simulations of 10-residue peptide, chignolin, using statistical mechanics: Relaxation mode analysis and three-dimensional reference interaction site model theory
Gao et al. NMR assignments of sparsely labeled proteins using a genetic algorithm
Álvarez et al. Predicting protein tertiary structure and its uncertainty analysis via particle swarm sampling
Meshkin et al. Toward convergence in free energy calculations for protein conformational changes: A case study on the thin gate of Mhp1 transporter
Chen et al. Free-energy calculations along a high-dimensional fragmented path with constrained dynamics
Anile et al. Determination of protein structure and dynamics combining immune algorithms and pattern search methods

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
WD01 Invention patent application deemed withdrawn after publication

Application publication date: 20201009

WD01 Invention patent application deemed withdrawn after publication