CN111161810B - 一种基于约束概率分布函数优化的自由能微扰方法 - Google Patents

一种基于约束概率分布函数优化的自由能微扰方法 Download PDF

Info

Publication number
CN111161810B
CN111161810B CN201911412608.4A CN201911412608A CN111161810B CN 111161810 B CN111161810 B CN 111161810B CN 201911412608 A CN201911412608 A CN 201911412608A CN 111161810 B CN111161810 B CN 111161810B
Authority
CN
China
Prior art keywords
probability distribution
constraint
free energy
calculation
delta
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.)
Active
Application number
CN201911412608.4A
Other languages
English (en)
Other versions
CN111161810A (zh
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.)
Sun Yat Sen University
Original Assignee
Sun Yat Sen University
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 Sun Yat Sen University filed Critical Sun Yat Sen University
Priority to CN201911412608.4A priority Critical patent/CN111161810B/zh
Priority to PCT/CN2020/071663 priority patent/WO2020147668A1/zh
Publication of CN111161810A publication Critical patent/CN111161810A/zh
Application granted granted Critical
Publication of CN111161810B publication Critical patent/CN111161810B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
    • G16C—COMPUTATIONAL CHEMISTRY; CHEMOINFORMATICS; COMPUTATIONAL MATERIALS SCIENCE
    • G16C20/00—Chemoinformatics, i.e. ICT specially adapted for the handling of physicochemical or structural data of chemical particles, elements, compounds or mixtures
    • G16C20/50—Molecular design, e.g. of drugs

Landscapes

  • Engineering & Computer Science (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Chemical & Material Sciences (AREA)
  • Medicinal Chemistry (AREA)
  • Health & Medical Sciences (AREA)
  • General Health & Medical Sciences (AREA)
  • Physics & Mathematics (AREA)
  • Pharmacology & Pharmacy (AREA)
  • Spectroscopy & Molecular Physics (AREA)
  • Crystallography & Structural Chemistry (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Computing Systems (AREA)
  • Theoretical Computer Science (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

本发明提供一种基于约束概率分布函数优化的自由能微扰方法，包括：根据热力学循环构建自由能计算分子模拟窗口，自动化添加约束，对于每个窗口进行分子动力学模拟，统计得出相邻窗口之间的势能差，对于添加约束步骤的势能差统计得到概率分布，基于约束概率分布函数针对添加约束过程的概率分布的拟合，根据拟合后的概率分布重新生成势能差，最终统计得出药物‑靶标绝对结合自由能。本发明提供的自由能微扰方法，实现全自动得到添加约束所需的复合物拓扑文件，约束概率分布函数用于拟合实验中的概率分布，使FEP得到的约束窗口内的自由能变化自由值更精准，在保证计算精度的前提下，极大减少了约束步骤的窗口数，减少了计算量，提高了计算效率。

Description

一种基于约束概率分布函数优化的自由能微扰方法
技术领域
本发明涉及药物研发技术领域，更具体的，涉及一种基于约束概率分布函数优化的自由能微扰方法。
背景技术
计算药物小分子与靶标蛋白的结合能力一直是基于结构的药物设计的核心目的。一种能准确并快速预测药物-靶标结合能力的方法，可以高效经济地指导药物化学家对先导化合物的结构优化，进而大幅度降低新药研发的周期和成本。
迄今为止，有许多计算技术尝试估计或精准计算药物-靶标结合自由能的技术，包括如对接打分函数、基于隐性溶剂模型的MM/PBSA(分子力学/泊松玻尔兹曼表面)、MM/GBSA(分子力学/广义伯恩表面)等方法。但由于这些方法通常是基于实验数据，添加相应的经验常数，导致其计算结果与实际的实验数据相差较大。著名药物设计软件公司薛定谔Schrodinger发展了一套基于自由能微扰的相对结合自由能预测方法(J.Am.Chem.Soc.2015,137,2695-2703.)。然而该方法只能用于预测结构差异较小的同骨架化合物与靶标蛋白的的相对结合自由能变化(差异需小于10个原子)。
适用于薛定谔公司的相对结合自由能计算方法的计算体系中药物小分子间的结构差异较小(差异小于10个原子)，不同小分子在自由能微扰过程中产生的误差相差较小。但当配体小分子的结构差异性较大时，差异性大的小分子在自由能微扰过程中产生的误差较大，会导致计算结果收敛性变差，不能准确预测不同骨架结构的药物分子与靶标蛋白结合能力。Martin Karplus提出新的热力学循环(J.Phys.Chem.B 2003,107,9535-9551)，通过引入配体与受体之间的约束步骤，约束小分子与靶标蛋白的相对位置，提高自由能微扰计算方法的收敛性。加入约束步骤后的完整热力学循环如图2。图中Rec代表标靶受体；Lig代表药物小分子配体；Dumm代表虚原子，所有的虚原子与周围的环境没有任何的相互作用；Restrain代表对药物小分子施加简谐势约束。
Figure GDA0003466482030000011
表示将溶液中药物小分子全部微扰至湮灭状态的能量；
Figure GDA0003466482030000012
表示将药物-靶标复合物中药物小分子全部微扰至湮灭状态的能量；
Figure GDA0003466482030000013
表示对药物-靶标复合物中药物小分子添加简谐势约束的能量变化，
Figure GDA0003466482030000021
表示对溶液中的药物小分子添加相应约束的能量变化。其中，湮灭状态定义为药物小分子的电荷以及范德华力L-J势能为0，周围原子无法感知到药物小分子的存在。
但现有的基于该热力学循环的绝对结合自由能计算方法(Chem.Sci.,2016,7,207-218；J.Am.Chem.Soc.2017,139,946-957等)，在计算
Figure GDA0003466482030000022
这一项过程中需设置11个扰动窗口，能量计算收敛差，计算耗时长，计算效率低。
发明内容
本发明为克服上述现有基于热力学循环的绝对结合自由能计算方法存在能量计算收敛差，计算耗时长，计算效率低的技术缺陷，提供一种基于约束概率分布函数优化的自由能微扰方法。
为解决上述技术问题，本发明的技术方案如下：
一种基于约束概率分布函数优化的自由能微扰方法，包括以下步骤：
S1：对受体蛋白-配体体系，即Rec-Lig体系进行常规的分子动力学模拟，并基于模拟获得的轨迹得到添加约束所需的距离、夹角、二面角数据；
S2：构建热力学循环，构建
Figure GDA0003466482030000023
和
Figure GDA0003466482030000024
所需的分子动力学模拟参数，
Figure GDA0003466482030000025
所需的距离、夹角、二面角数据添加约束得到，并对
Figure GDA0003466482030000026
的自由能进行解析计算；
其中：
Figure GDA0003466482030000027
表示将溶液中药物小分子全部微扰至湮灭状态的能量；
Figure GDA0003466482030000028
表示将药物-靶标复合物中药物小分子全部微扰至湮灭状态的能量；
Figure GDA0003466482030000029
表示对药物-靶标复合物中药物小分子添加简谐势约束的能量变化；
Figure GDA00034664820300000210
表示对溶液中的药物小分子添加相应约束的能量变化；
S3：利用分子动力学模拟程序，对
Figure GDA00034664820300000211
步骤中的每一个窗口进行分子动力学模拟采样，得到
Figure GDA00034664820300000212
Figure GDA00034664820300000213
步骤中每一个窗口的分子动力学轨迹；
S4：根据
Figure GDA00034664820300000214
步骤中的每一个窗口的分子动力学轨迹，通过相邻窗口的力场文件进行计算，得到每个相邻窗口之间的势能差ΔU；
S5：对于
Figure GDA00034664820300000215
这两步，根据对应窗口的ΔU，分别计算得到
Figure GDA00034664820300000216
步骤的自由能；
S6：对于
Figure GDA00034664820300000217
这一步每个相邻窗口之间的势能差ΔU进行概率分布统计，利用约束概率分布函数进行拟合，得到平滑的概率分布，并根据拟合结果检查收敛性；
S7：根据拟合后的概率分布，重新生成
Figure GDA0003466482030000031
这一步中每一窗口的势能差，通过计算得到
Figure GDA0003466482030000032
的能量，并根据之前步骤得到的
Figure GDA0003466482030000033
Figure GDA0003466482030000034
以及
Figure GDA0003466482030000035
计算得到药-靶结合自由能ΔAbinding。
其中，所述步骤S1具体包括以下步骤：
S11：准备受体蛋白Rec-配体Lig复合物坐标文件；
S12：准备小分子配体的电荷，并对配体指定力场；
S13：将受体-配体复合物中的受体指定力场，并将该复合物在水箱中进行浸泡，准备好复合物进行分子模拟的初始体系，记为Rec-Lig体系，并导出力场文件作为Rec-Lig体系初始状态参数；
S14：将单独小分子配体浸泡在水箱中，准备好单独的小分子的分子模拟初始体系，记为Lig体系，并导出力场文件作为Lig体系初始状态参数；
S15：根据Rec-Lig体系初始状态的力场参数，进行短时间的分子动力学模拟，得到运动轨迹；
S16：基于模拟得到的轨迹，统计与小分子配体相对位置变化较小的受体蛋白的氨基酸残基，取该残基中与小分子配体相对距离最稳定的三个原子作为添加约束所需的属于受体蛋白的三个原子，选取小分子配体中最稳定的三个原子作为添加约束所需的属于配体小分子的三个原子，得到添加约束所需的距离、夹角、二面角数据。
其中，所述步骤S2具体包括以下步骤：
S21：根据准备好的Rec-Lig体系的力场参数，根据步骤S16得到的6个原子，添加一个长度、两个夹角、三个二面角约束，得到计算
Figure GDA0003466482030000036
所需窗口的分子动力学模拟参数；
S22：进一步将Rec-Lig体系中的小分子的电荷以及范德华参数依次分别减小至0，得到计算
Figure GDA0003466482030000037
所需窗口的分子动力学模拟参数；
S23：对于Rec-Lig体系中小分子的电荷以及范德华参数依次分别减小至0，得到计算
Figure GDA0003466482030000038
所需窗口的分子动力学模拟参数；
S24：对
Figure GDA0003466482030000039
的自由能进行解析计算。
其中，在所述步骤S24中，所述
Figure GDA00034664820300000310
自由能的计算公式具体为：
Figure GDA0003466482030000041
其中，kB为玻尔兹曼常数，V为系统的体积，raA为一个两原子间的距离用于距离简谐势约束，θa和θA为两个三原子间的夹角用于夹角简谐势约束，诸多k值为用于约束一个距离raA、两个夹角θa,θA、三个二面角
Figure GDA0003466482030000042
的简谐势力常数。
其中，所述步骤S5具体为：
根据对应窗口的ΔU，采用BenneetAcceptance Ratio计算方法，即BAR计算方法进行
Figure GDA0003466482030000043
自由能的计算，具体计算公式为：
Figure GDA0003466482030000044
Figure GDA0003466482030000045
其中，ΔAi,i+1为两相邻窗口i与i+1的结合自由能差，ΔUi,i+1为基于拟合后的概率分布新生成的势能差，β＝1/kBT，kB为玻尔兹曼常数，ni与ni+1分别为两相邻状态i与i+1采样所包含的样本量，C为方程组之间的隐含变量；
将
Figure GDA0003466482030000046
和
Figure GDA0003466482030000047
对应的计算窗口得到的ΔAi,i+1进行加和，得到
Figure GDA0003466482030000048
和
Figure GDA0003466482030000049
的数值。
其中，所述步骤S6具体包括以下步骤：
S61：对计算
Figure GDA00034664820300000410
的窗口中每个窗口之间的势能差进行概率分布统计，得到概率分布P(ΔU)data；
S62：基于约束概率分布函数对概率分布P(ΔU)data进行拟合，具体表达式为：
Figure GDA00034664820300000411
其中，P(ΔU)fit为经过约束概率分布函数拟合后的平滑的概率分布，其中ai,bi,ci,di,hi,μi,σi为第i个“约束概率分布函数”的参数，n为选取拟合函数的数目，n1和n2为两个正整数的修正因子。
其中，所述步骤S7具体包括以下步骤：
S71：根据概率分布P(ΔU)fit，重新生成每一状态的势能差ΔU′；
S72：基于新的势能差ΔU′，利用步骤S5的计算方法重新计算得到
Figure GDA0003466482030000051
的计算数值；
S73：根据步骤S24、步骤S5和步骤S72得到的
Figure GDA0003466482030000052
Figure GDA0003466482030000053
和
Figure GDA0003466482030000054
计算出药-靶结合自由能ΔAbinding。
其中，在所述步骤S73中，所述药-靶结合自由能ΔAbinding的计算过程具体为：
Figure GDA0003466482030000055
与现有技术相比，本发明技术方案的有益效果是：
本发明提供的一种基于约束概率分布函数优化的自由能微扰方法，实现全自动得到添加约束所需的复合物拓扑文件，约束概率分布函数用于拟合实验中的概率分布，使自由能微扰(FEP)得到的约束窗口内的自由能变化自由值更精准，在保证计算精度的前提下，极大减少了约束步骤的窗口数，减少了计算量，提高了计算效率。
附图说明
图1为本发明所述方法流程示意图；
图2为对添加约束过程的状态采样所得概率被约束概率分布函数拟合的结果示意图。
图3为计算绝对结合自由能所涉及的热力学循环示意图。
具体实施方式
附图仅用于示例性说明，不能理解为对本专利的限制；
为了更好说明本实施例，附图某些部件会有省略、放大或缩小，并不代表实际产品的尺寸；
对于本领域技术人员来说，附图中某些公知结构及其说明可能省略是可以理解的。
下面结合附图和实施例对本发明的技术方案做进一步的说明。
实施例1
如图1所示，一种基于约束概率分布函数优化的自由能微扰方法，包括以下步骤：
S1：对受体蛋白-配体体系，即Rec-Lig体系进行常规的分子动力学模拟，并基于模拟获得的轨迹得到添加约束所需的距离、夹角、二面角数据；
S2：构建热力学循环，构建
Figure GDA0003466482030000056
和
Figure GDA0003466482030000057
所需的分子动力学模拟参数；
Figure GDA0003466482030000061
根据所需的距离、夹角、二面角数据添加约束得到，并对
Figure GDA0003466482030000062
的自由能进行解析计算；
其中：
Figure GDA0003466482030000063
表示将溶液中药物小分子全部微扰至湮灭状态的能量；
Figure GDA0003466482030000064
表示将药物-靶标复合物中药物小分子全部微扰至湮灭状态的能量；
Figure GDA0003466482030000065
表示对药物-靶标复合物中药物小分子添加简谐势约束的能量变化；
Figure GDA0003466482030000066
表示对溶液中的药物小分子添加相应约束的能量变化；
S3：利用分子动力学模拟程序，对
Figure GDA0003466482030000067
步骤中的每一个窗口进行分子动力学模拟采样，得到
Figure GDA0003466482030000068
Figure GDA0003466482030000069
步骤中每一个窗口的分子动力学轨迹；
S4：根据
Figure GDA00034664820300000610
步骤中的每一个窗口的分子动力学轨迹，通过相邻窗口的力场文件进行计算，得到每个相邻窗口之间的势能差ΔU；
S5：对于
Figure GDA00034664820300000611
这两步，根据对应窗口的ΔU，分别计算得到
Figure GDA00034664820300000612
步骤的自由能；
S6：对于
Figure GDA00034664820300000613
这一步每个相邻窗口之间的势能差ΔU进行概率分布统计，利用约束概率分布函数进行拟合，得到平滑的概率分布，并根据拟合结果检查收敛性；
S7：根据拟合后的概率分布，重新生成
Figure GDA00034664820300000614
这一步中每一窗口的势能差，通过计算得到
Figure GDA00034664820300000615
的能量，并根据之前步骤得到的
Figure GDA00034664820300000616
Figure GDA00034664820300000617
以及
Figure GDA00034664820300000618
计算得到药-靶结合自由能ΔAbinding。
在具体实施过程中，本发明提供的一种基于约束概率分布函数优化的自由能微扰方法，实现全自动得到添加约束所需的复合物拓扑文件，约束概率分布函数用于拟合实验中的概率分布，使自由能微扰(FEP)得到的约束窗口内的自由能变化自由值更精准，在保证计算精度的前提下，极大减少了约束步骤的窗口数，减少了计算量，提高了计算效率。
在具体实施过程中，根据以上步骤，基于本发明提供的自动添加约束方法和约束概率分布函数拟合约束步骤的概率分布，以人类磷酸二酯酶4D及其小分子抑制剂罗氟司特的复合物体系(PDB代码为1XOQ)和人类BRD4蛋白及其小分子配体阿普唑仑(PDB代码为3U5J)为例，将约束窗口由现有文献中使用的10个(Chem.Sci.,2016,7,207-218；J.Am.Chem.Soc.2017,139,946-957等)减少到1个，通过对比考察10步方案计算所得的
Figure GDA0003466482030000071
与1步方案差距，获得的测试结果如表1所示。
表1.选择不同数目的窗口计算得出的
Figure GDA0003466482030000072
Figure GDA0003466482030000073
从上表可以看出，在都以本发明提供的自动添加约束程序进行简谐势约束的前提下，设置不同的约束步骤扰动窗口数，经过“约束概率分布函数”拟合其概率分布后，计算得到
Figure GDA0003466482030000074
相差较小，差值仅为0.2kcal/mol以内。该对比说明在使用本发明提供的自动添加约束程序及“约束概率分布函数”拟合基础上，只设置1个约束扰动窗口的快速约束方案是可行的。
更具体的，所述步骤S1具体包括以下步骤：
S11：准备受体蛋白Rec-配体Lig复合物坐标文件；
S12：准备小分子配体的电荷，并对配体指定力场；
S13：将受体-配体复合物中的受体指定力场，并将该复合物在水箱中进行浸泡，准备好复合物进行分子模拟的初始体系，记为Rec-Lig体系，并导出力场文件作为Rec-Lig体系初始状态参数；
S14：将单独小分子配体浸泡在水箱中，准备好单独的小分子的分子模拟初始体系，记为Lig体系，并导出力场文件作为Lig体系初始状态参数；
S15：根据Rec-Lig体系初始状态的力场参数，进行短时间的分子动力学模拟，得到运动轨迹；
S16：基于模拟得到的轨迹，统计与小分子配体相对位置变化较小的受体蛋白的氨基酸残基，取该残基中与小分子配体相对距离最稳定的三个原子作为添加约束所需的属于受体蛋白的三个原子，选取小分子配体中最稳定的三个原子作为添加约束所需的属于配体小分子的三个原子，得到添加约束所需的距离、夹角、二面角数据。
更具体的，所述步骤S2具体包括以下步骤：
S21：根据准备好的Rec-Lig体系的力场参数，根据步骤S16得到的6个原子，添加一个长度、两个夹角、三个二面角约束，得到计算
Figure GDA0003466482030000081
所需窗口的分子动力学模拟参数；
S22：进一步将Rec-Lig体系中的小分子的电荷以及范德华参数依次分别减小至0，得到计算
Figure GDA0003466482030000082
所需窗口的分子动力学模拟参数；
S23：对于Rec-Lig体系中小分子的电荷以及范德华参数依次分别减小至0，得到计算
Figure GDA0003466482030000083
所需窗口的分子动力学模拟参数；
S24：对
Figure GDA0003466482030000084
的自由能进行解析计算。
其中，在所述步骤S24中，所述
Figure GDA0003466482030000085
自由能的计算公式具体为：
Figure GDA0003466482030000086
其中，kB为玻尔兹曼常数，V为系统的体积，raA为一个两原子间的距离用于距离简谐势约束，θa和θA为两个三原子间的夹角用于夹角简谐势约束，诸多k值为用于约束一个距离raA、两个夹角θa,θA、三个二面角
Figure GDA0003466482030000087
的简谐势力常数。
更具体的，所述步骤S5具体为：
根据对应窗口的ΔU，采用BenneetAcceptance Ratio计算方法，即BAR计算方法进行
Figure GDA0003466482030000088
自由能的计算，具体计算公式为：
Figure GDA0003466482030000089
Figure GDA00034664820300000810
其中，ΔAi,i+1为两相邻窗口i与i+1的结合自由能差，ΔUi,i+1为基于拟合后的概率分布新生成的势能差，β＝1/kBT，kB为玻尔兹曼常数，ni与ni+1分别为两相邻状态i与i+1采样所包含的样本量，C为方程组之间的隐含变量；
将
Figure GDA00034664820300000811
和
Figure GDA00034664820300000812
对应的计算窗口得到的ΔAi,i+1进行加和，得到
Figure GDA00034664820300000813
和
Figure GDA00034664820300000814
的数值。
更具体的，所述步骤S6具体包括以下步骤：
S61：对计算
Figure GDA00034664820300000815
的窗口中每个窗口之间的势能差进行概率分布统计，得到概率分布P(ΔU)data；
S62：基于约束概率分布函数对概率分布P(ΔU)data进行拟合，具体表达式为：
Figure GDA0003466482030000091
其中，P(ΔU)fit为经过约束概率分布函数拟合后的平滑的概率分布，其中ai,bi,ci,di,hi,μi,σi为第i个“约束概率分布函数”的参数，n为选取拟合函数的数目，n1和n2为两个正整数的修正因子，在本申请给出的计算案例中，取n，n1和n2分别为1，10和4。为说明拟合效果，图2为对添加约束过程的状态采样所得概率分布P(ΔU)被一系列约束概率分布函数拟合的结果示意图。图2中灰点为采样所得概率分布点，黑线为最终拟合所得的分布曲线，灰实线和灰虚线分别为公式中的前后两项。
更具体的，所述步骤S7具体包括以下步骤：
S71：根据概率分布P(ΔU)fit，重新生成每一状态的势能差ΔU′；
S72：基于新的势能差ΔU′，利用步骤S5的计算方法重新计算得到
Figure GDA0003466482030000092
的计算数值；
S73：根据步骤S24、步骤S5和步骤S72得到的
Figure GDA0003466482030000093
Figure GDA0003466482030000094
和
Figure GDA0003466482030000095
计算出药-靶结合自由能ΔAbinding。
其中，在所述步骤S73中，所述药-靶结合自由能ΔAbinding的计算过程具体为：
Figure GDA0003466482030000096
实施例2
实施例1的结果说明自动添加约束的方法结合约束概率分布函数进行拟合，可以使得
Figure GDA0003466482030000097
项能量结算结果达到非常好的收敛性。基于实施例1的基础上，以两个针对BRD4蛋白靶标的药物小分子为例。采用自动添加约束并在约束过程中只设置一个扰动窗口的约束方案，在处理约束过程中的概率分布数据时，采用本发明中的约束概率分布函数对其概率分布散点图进行拟合。并基于图3提供的完整版热力学循环，计算绝对结合自由能ΔAbinding。具体步骤如下：
S1：下载BRD4蛋白晶体结构3MXF和4HBV，准备BRD4的力场参数，选定Rec-Lig体系以及Lig体系，对Rec-Lig体系进行一小段常规的分子动力学模拟，并基于模拟获得的轨迹得到添加约束所需距离、夹角、二面角数据；
S2：根据图2所示构建热力学循环，构建
Figure GDA0003466482030000098
步骤所需的分子动力学模拟参数，其中对
Figure GDA0003466482030000099
这一步添加约束使用S1步骤得到的距离、夹角、二面角数据，并对
Figure GDA0003466482030000101
这一步的自由能进行解析解计算；
S3：利用分子动力学模拟程序，对
Figure GDA0003466482030000102
步骤中的每一个窗口进行分子动力学模拟采样，得到
Figure GDA0003466482030000103
Figure GDA0003466482030000104
步骤中的每一个窗口的分子动力学轨迹；
S4：根据
Figure GDA0003466482030000105
步骤中的每一个窗口的分子动力学轨迹，通过相邻窗口的力场文件进行计算，得到每个相邻窗口之间的势能差ΔU；
S5：对于
Figure GDA0003466482030000106
这两步，根据对应窗口的ΔU，分别计算得到
Figure GDA0003466482030000107
步骤的自由能。
S6：对于
Figure GDA0003466482030000108
这一步每个相邻窗口之间的势能差ΔU进行概率分布统计，利用本专利给出的“约束概率分布函数”进行拟合，得到平滑的概率分布，并根据拟合结果检查收敛性；
S7：根据拟合后的概率分布，重新生成
Figure GDA0003466482030000109
这一步中每一窗口的势能差，通过计算得到
Figure GDA00034664820300001010
的能量，并根据之前步骤得到的
Figure GDA00034664820300001011
Figure GDA00034664820300001012
以及
Figure GDA00034664820300001013
计算得到药-靶结合自由能ΔAbinding。
在具体实施过程中，根据以上步骤，考察针对BRD4蛋白靶标与两个两个药物小分子的计算结合自由能ΔAbinding与实验值的差别，获得的结果如表2所示：
表2.药-靶结合自由能计算结果与实验结果对比a
Figure GDA00034664820300001014
a计算过程中使用自动添加约束方法
bΔAcal为计算所得药靶结合自由能
cΔAexp为实验所得药靶结合自由能
从计算与实验结果的对比可知，计算结果与实验结果非常接近，差距在2kcal/mol之内。
在最新文献中，计算
Figure GDA00034664820300001015
至少需10个扰动窗口，且需要手动添加约束(Chem.Sci.，2016，7，207-218；J.Am.Chem.Soc.2017，139，946-957)。基于本发明的技术方案，可实现全自动得到添加简谐势约束所需的复合物拓扑文件，约束概率分布函数可用于拟合实验所得的ΔU与其概率的散点图，使得FEP计算得到的约束窗口内的自由能变化值更精准。FEP方法是目前准确度较高的方法，但该方法需要非常大的计算量，而使用本发明的技术方案，可以在保证计算精确度的前提下，极大减少约束步骤的窗口数，使得计算窗口数目从10降至1，这同时也极大减少了计算量，计算量仅为之前的1/10。两个实施例也说明了结果的准确性。
显然，本发明的上述实施例仅仅是为清楚地说明本发明所作的举例，而并非是对本发明的实施方式的限定。对于所属领域的普通技术人员来说，在上述说明的基础上还可以做出其它不同形式的变化或变动。这里无需也无法对所有的实施方式予以穷举。凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明权利要求的保护范围之内。

Claims (8)

1.一种基于约束概率分布函数优化的自由能微扰方法，其特征在于，包括以下步骤：
S1：对受体蛋白-配体体系，即Rec-Lig体系进行常规的分子动力学模拟，并基于模拟获得的轨迹得到添加约束所需的距离、夹角、二面角数据；
S2：构建热力学循环，构建
Figure FDA0003466482020000011
和
Figure FDA0003466482020000012
所需的分子动力学模拟参数，
Figure FDA0003466482020000013
根据所需的距离、夹角、二面角数据添加约束得到，并对
Figure FDA0003466482020000014
的自由能进行解析计算；
其中：
Figure FDA0003466482020000015
表示将溶液中药物小分子全部微扰至湮灭状态的能量；
Figure FDA0003466482020000016
表示将药物-靶标复合物中药物小分子全部微扰至湮灭状态的能量；
Figure FDA0003466482020000017
表示对药物-靶标复合物中药物小分子添加简谐势约束的能量变化；
Figure FDA0003466482020000018
表示对溶液中的药物小分子添加相应约束的能量变化；
S3：利用分子动力学模拟程序，对
Figure FDA0003466482020000019
步骤中的每一个窗口进行分子动力学模拟采样，得到
Figure FDA00034664820200000110
Figure FDA00034664820200000111
步骤中每一个窗口的分子动力学轨迹；
S4：根据
Figure FDA00034664820200000112
步骤中的每一个窗口的分子动力学轨迹，通过相邻窗口的力场文件进行计算，得到每个相邻窗口之间的势能差ΔU；
S5：对于
Figure FDA00034664820200000113
这两步，根据对应窗口的ΔU，分别计算得到
Figure FDA00034664820200000114
步骤的自由能；
S6：对于
Figure FDA00034664820200000115
这一步每个相邻窗口之间的势能差ΔU进行概率分布统计，利用约束概率分布函数进行拟合，得到平滑的概率分布，并根据拟合结果检查收敛性；
S7：根据拟合后的概率分布，重新生成
Figure FDA00034664820200000116
这一步中每一窗口的势能差，通过计算得到
Figure FDA00034664820200000117
的能量，并根据之前步骤得到的
Figure FDA00034664820200000118
Figure FDA00034664820200000119
以及
Figure FDA00034664820200000120
计算得到药-靶结合自由能ΔAbinding。
2.根据权利要求1所述的一种基于约束概率分布函数优化的自由能微扰方法，其特征在于，所述步骤S1具体包括以下步骤：
S11：准备受体蛋白Rec-配体Lig复合物坐标文件；
S12：准备小分子配体的电荷，并对配体指定力场；
S13：将受体-配体复合物中的受体指定力场，并将该复合物在水箱中进行浸泡，准备好复合物进行分子模拟的初始体系，记为Rec-Lig体系，并导出力场文件作为Rec-Lig体系初始状态参数；
S14：将单独小分子配体浸泡在水箱中，准备好单独的小分子的分子模拟初始体系，记为Lig体系，并导出力场文件作为Lig体系初始状态参数；
S15：根据Rec-Lig体系初始状态的力场参数，进行短时间的分子动力学模拟，得到运动轨迹；
S16：基于模拟得到的轨迹，统计与小分子配体相对位置变化符合要求的受体蛋白的氨基酸残基，取该残基中与小分子配体相对距离最稳定的三个原子作为添加约束所需的属于受体蛋白的三个原子，选取小分子配体中最稳定的三个原子作为添加约束所需的属于配体小分子的三个原子，得到添加约束所需的距离、夹角、二面角数据。
3.根据权利要求2所述的一种基于约束概率分布函数优化的自由能微扰方法，其特征在于，所述步骤S2具体包括以下步骤：
S21：根据准备好的Rec-Lig体系的力场参数，根据步骤S16得到的6个原子，添加一个长度、两个夹角、三个二面角约束，得到计算
Figure FDA0003466482020000021
所需窗口的分子动力学模拟参数；
S22：进一步将Rec-Lig体系中的小分子的电荷以及范德华参数依次分别减小至0，得到计算
Figure FDA0003466482020000022
所需窗口的分子动力学模拟参数；
S23：对于Rec-Lig体系中小分子的电荷以及范德华参数依次分别减小至0，得到计算
Figure FDA0003466482020000023
所需窗口的分子动力学模拟参数；
S24：对
Figure FDA0003466482020000024
的自由能进行解析计算。
4.根据权利要求3所述的一种基于约束概率分布函数优化的自由能微扰方法，其特征在于，在所述步骤S24中，所述
Figure FDA0003466482020000025
自由能的计算公式具体为：
Figure FDA0003466482020000026
其中，kB为玻尔兹曼常数，V为系统的体积，raA为一个两原子间的距离用于距离简谐势约束，θa和θA为两个三原子间的夹角用于夹角简谐势约束，诸多k值为用于约束一个距离raA、两个夹角θa,θA、三个二面角
Figure FDA0003466482020000027
的简谐势力常数。
5.根据权利要求3所述的一种基于约束概率分布函数优化的自由能微扰方法，其特征在于，所述步骤S5具体为：
根据对应窗口的ΔU，采用BenneetAcceptance Ratio计算方法，即BAR计算方法进行
Figure FDA0003466482020000031
自由能的计算，具体计算公式为：
Figure FDA0003466482020000032
Figure FDA0003466482020000033
其中，ΔAi,i+1为两相邻窗口i与i+1的结合自由能差，ΔUi,i+1为基于拟合后的概率分布新生成的势能差，β＝1/kBT，kB为玻尔兹曼常数，ni与ni+1分别为两相邻状态i与i+1采样所包含的样本量，C为方程组之间的隐含变量；
将
Figure FDA0003466482020000034
和
Figure FDA0003466482020000035
对应的计算窗口得到的ΔAi,i+1进行加和，得到
Figure FDA0003466482020000036
和
Figure FDA0003466482020000037
的数值。
6.根据权利要求5所述的一种基于约束概率分布函数优化的自由能微扰方法，其特征在于，所述步骤S6具体包括以下步骤：
S61：对计算
Figure FDA0003466482020000038
的窗口中每个窗口之间的势能差进行概率分布统计，得到概率分布P(ΔU)data；
S62：基于约束概率分布函数对概率分布P(ΔU)data进行拟合，具体表达式为：
Figure FDA0003466482020000039
其中，P(ΔU)fit为经过约束概率分布函数拟合后的平滑的概率分布，其中ai,bi,ci,di,hi,μi,σi为第i个“约束概率分布函数”的参数，n为选取拟合函数的数目，n1和n2为两个正整数的修正因子。
7.根据权利要求6所述的一种基于约束概率分布函数优化的自由能微扰方法，其特征在于，所述步骤S7具体包括以下步骤：
S71：根据概率分布P(ΔU)fit，重新生成每一状态的势能差ΔU′；
S72：基于新的势能差ΔU′，利用步骤S5的计算方法重新计算得到
Figure FDA00034664820200000310
的计算数值；
S73：根据步骤S24、步骤S5和步骤S72得到的
Figure FDA0003466482020000041
Figure FDA0003466482020000042
和
Figure FDA0003466482020000043
计算出药-靶结合自由能ΔAbinding。
8.根据权利要求7所述的一种基于约束概率分布函数优化的自由能微扰方法，其特征在于，在所述步骤S73中，所述药-靶结合自由能ΔAbinding的计算过程具体为：
Figure FDA0003466482020000044
CN201911412608.4A 2019-01-17 2019-12-31 一种基于约束概率分布函数优化的自由能微扰方法 Active CN111161810B (zh)

Priority Applications (2)

Application Number Priority Date Filing Date Title
CN201911412608.4A CN111161810B (zh) 2019-12-31 2019-12-31 一种基于约束概率分布函数优化的自由能微扰方法
PCT/CN2020/071663 WO2020147668A1 (zh) 2019-01-17 2020-01-13 预测药-靶结合强度的绝对结合自由能微扰方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201911412608.4A CN111161810B (zh) 2019-12-31 2019-12-31 一种基于约束概率分布函数优化的自由能微扰方法

Publications (2)

Publication Number Publication Date
CN111161810A CN111161810A (zh) 2020-05-15
CN111161810B true CN111161810B (zh) 2022-03-22

Family

ID=70560036

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201911412608.4A Active CN111161810B (zh) 2019-01-17 2019-12-31 一种基于约束概率分布函数优化的自由能微扰方法

Country Status (1)

Country Link
CN (1) CN111161810B (zh)

Families Citing this family (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112199909A (zh) * 2020-10-22 2021-01-08 深圳晶泰科技有限公司 一种准确计算气体分子绝对自由能的方法
CN112216350B (zh) * 2020-11-05 2022-09-13 深圳晶泰科技有限公司 物理严格且相空间重叠最大化的相对自由能计算方法
CN114187971B (zh) * 2021-12-10 2023-03-28 上海智药科技有限公司 分子自由能计算、稳定性分析方法、装置、设备及存储介质
WO2023123396A1 (zh) * 2021-12-31 2023-07-06 深圳晶泰科技有限公司 增强采样方法及计算复合物的结合自由能的方法

Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2004102159A2 (en) * 2003-05-13 2004-11-25 The Penn State Research Foundation Quantum mechanics based method for scoring protein-ligand interactions
CN102930152A (zh) * 2012-10-26 2013-02-13 中国科学院上海药物研究所 一种模拟配体分子与靶标受体反应并计算预测该反应的热力学与动力学参数的方法和系统
CN107423570A (zh) * 2017-08-02 2017-12-01 南昌立德生物技术有限公司 快速准确计算蛋白酶与药物分子之间亲和自由能的算法
CN109859806A (zh) * 2019-01-17 2019-06-07 中山大学 一种预测药物-靶标结合强度的绝对自由能微扰方法
CN110047559A (zh) * 2019-03-06 2019-07-23 山东师范大学 蛋白质与药物结合自由能的计算方法、系统、设备及介质
CN110400598A (zh) * 2019-07-03 2019-11-01 江苏理工学院 基于mm/pbsa模型的蛋白质-配体结合自由能计算方法

Family Cites Families (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2015002860A1 (en) * 2013-07-02 2015-01-08 Epigenetx, Llc Structure-based modeling and target-selectivity prediction
WO2015066415A1 (en) * 2013-11-01 2015-05-07 University Of Florida Research Foundation, Inc. Movable type method applied to protein-ligand binding

Patent Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2004102159A2 (en) * 2003-05-13 2004-11-25 The Penn State Research Foundation Quantum mechanics based method for scoring protein-ligand interactions
CN102930152A (zh) * 2012-10-26 2013-02-13 中国科学院上海药物研究所 一种模拟配体分子与靶标受体反应并计算预测该反应的热力学与动力学参数的方法和系统
CN107423570A (zh) * 2017-08-02 2017-12-01 南昌立德生物技术有限公司 快速准确计算蛋白酶与药物分子之间亲和自由能的算法
CN109859806A (zh) * 2019-01-17 2019-06-07 中山大学 一种预测药物-靶标结合强度的绝对自由能微扰方法
CN110047559A (zh) * 2019-03-06 2019-07-23 山东师范大学 蛋白质与药物结合自由能的计算方法、系统、设备及介质
CN110400598A (zh) * 2019-07-03 2019-11-01 江苏理工学院 基于mm/pbsa模型的蛋白质-配体结合自由能计算方法

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
Absolute Binding Free Energy Calculation and Design of a Subnanomolar Inhibitor of Phosphodiesterase-10;ZHE LI ,et al;《 JOURNAL OF MEDICINAL CHEMISTRY》;20190128;第2099-2111页 *
高斯函数增强的自由能微扰方法（GA-FEP）的开发及其在药物设计中的应用;李哲;《2019中国化学会第十五届全国计算（机）化学学术会议》;20191114;第1页 *

Also Published As

Publication number Publication date
CN111161810A (zh) 2020-05-15

Similar Documents

Publication Publication Date Title
CN111161810B (zh) 一种基于约束概率分布函数优化的自由能微扰方法
Savelyev et al. Competition among Li+, Na+, K+, and Rb+ monovalent ions for DNA in molecular dynamics simulations using the additive CHARMM36 and Drude polarizable force fields
Lan et al. Predicting drug–target interaction using positive-unlabeled learning
Jo et al. CHARMM-GUI Ligand Binder for absolute binding free energy calculations and its application
Liu et al. Combined MEDV-GA-MLR method for QSAR of three panels of steroids, dipeptides, and COX-2 inhibitors
Bajaj et al. Halide ion microhydration: Structure, energetics, and spectroscopy of small halide–water clusters
CN109859806B (zh) 一种预测药物-靶标结合强度的绝对自由能微扰方法
Hofer et al. Hydration of highly charged ions
Shi et al. Thermodynamic coupling of protonation and conformational equilibria in proteins: theory and simulation
Cole et al. Protein-protein interactions from linear-scaling first-principles quantum-mechanical calculations
González-Díaz et al. ANN-QSAR model for selection of anticancer leads from structurally heterogeneous series of compounds
CN110400598B (zh) 基于mm/pbsa模型的蛋白质-配体结合自由能计算方法
Kruse et al. QM computations on complete nucleic acids building blocks: analysis of the Sarcin–Ricin RNA motif using DFT-D3, HF-3c, PM6-D3H, and MM approaches
CN107423570B (zh) 快速准确计算蛋白酶与药物分子之间亲和自由能的算法
Flores et al. Multiscale modeling of macromolecular biosystems
Sfriso et al. Exploration of conformational transition pathways from coarse-grained simulations
Izairi et al. Comparison study of polar and nonpolar contributions to solvation free energy
Kantonen et al. Data-driven mapping of gas-phase quantum calculations to general force field Lennard-Jones parameters
Pan et al. A simplified charge projection scheme for long-range electrostatics in ab initio QM/MM calculations
Giese et al. Combined QM/MM, machine learning path integral approach to compute free energy profiles and kinetic isotope effects in RNA cleavage reactions
Chiappino-Pepe et al. Integration of metabolic, regulatory and signaling networks towards analysis of perturbation and dynamic responses
Gallegos et al. NNAIMQ: A neural network model for predicting QTAIM charges
Prado-Prado et al. Multi-target spectral moment: QSAR for antiviral drugs vs. different viral species
Suruzhon et al. Enhancing ligand and protein sampling using sequential monte carlo
Deng et al. Connecting free energy surfaces in implicit and explicit solvent: An efficient method to compute conformational and solvation free energies

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
GR01 Patent grant
GR01 Patent grant