一种基于MM/PB（GB）SA的蛋白质-药物结合自由能预测方法及预测系统

技术领域

本发明涉及药物筛选技术领域，尤其涉及一种基于MM/PB(GB)SA的蛋白质-药物结合自由能预测方法及预测系统。

背景技术

预测药物活性是当今新药研发中非常具有挑战性的一个课题。在生物体系中，结合自由能决定了许多生物进程的方向，比如：蛋白质折叠、酶催化以及药物靶标结合。因此，结合自由能预测在相关领域扮演着不可或缺的角色。在药物发现的过程中，药物与生物靶标的结合是药物发挥药效的基础，其亲和力的强弱直接决定了药物的生物学活性，而结合自由能又是药物与靶标亲和力大小的量化指标。因此，在先导化合物优化阶段，使用计算机预测候选药物与生物靶标亲和力的强弱，可以为先导化合物的优化提供理论指导，加速药物发现的进程。

目前，药物发现领域常常用来预测受体-配体相互作用的方法有自由能微扰(FEP)、热力学积分(TI)以及分子力学泊松-玻尔兹曼(广义伯恩)表面积(MM/PB(GB)SA)三种方法，相比于另外两种方法，MM/PB(GB)SA由于速度快、鲁棒性高以及计算资源消耗较少的优点，被广泛用于配体-受体自由能预测。

发明内容

针对上述背景技术，本发明解决的问题是：利用自由能预测技术在低计算资源消耗的前提下提升蛋白质药物靶标与药物配体结合自由能预测的准确度，从而节省药物研发的时间和经济成本。本研究在MM/PB(GB)SA理论的基础上，通过预先寻找合适的计算参数，来进一步提升MM/PB(GB)SA计算结合自由能的准确度。在测试阶段，本方法使用较低的计算资源消耗，取得了和FEP相媲美的准确度。

为实现上述目的，本发明采取的技术方案为：

一方面，本发明提供一种基于MM/PB(GB)SA的蛋白质-药物结合自由能预测方法，包括以下步骤：

(1)使用参数对所述蛋白质及其配体进行结合自由能预测，所述参数包括所述蛋白质的力场、所述配体的力场和所述配体的电荷数中的至少一种；

(2)将步骤(1)中预测的多组结合自由能分别与实验测得的结合自由能进行数据拟合；

(3)以步骤(2)中拟合最好的一组参数和方法分析所述蛋白质与药物的结合自由能；

其中，步骤(1)中，所述结合自由能预测的方法为MM/PBSA和/或MM/GBSA。

在本发明的技术方案中，所述MM/PB(GB)SA计算结合自由能的过程如下：

ΔG

＝ΔE

＝ΔE

其中，ΔG

作为优选地实施方式，步骤(1)中，所述配体与所述药物具有相同的化合物骨架以及相似的晶体结构，即所述配体与所述药物是由同一个先导化合物优化得到，或所述药物由所述配体优化得到，或所述配体由所述药物优化得到；其中，先导化合物即需要进行改良优化的化合物分子，由化合物骨架和至少一个官能团组成，化合物骨架是先导化合物分子的核心部分，而官能团是决定先导化合物分子的化学性质的原子或原子团，常见官能团包括羟基、羧基、醚键、醛基、羰基等。

在某些具体的实施方式中，MM/GBSA方法中，常用的GB模型包括GBHCT、GBOBC、GBOBC2、GBneck、GBneck2；

优选地，所述配体的电荷数的计算方法选自经验法、半经验法和量子化学计算法中的任一种；其中，基于经验的配体电荷计算方法，如：基于CHARMM普适力场的电荷计算方法；经验与量化相结合的半经验配体电荷计算方法，如：AM1-BCC原子电荷计算方法；基于量子化学的原子电荷计算方法，如基于密度泛函理论密度泛函理论(density functionaltheory，DFT)的静电势计算联用RESP(Restrained Electro Static Potential，有限制的静电电位)来拟合计算原子电荷数；

在本发明的技术方案中，所述配体的力场由配体的电荷数的计算方法决定，即，电荷的计算方法确定之后，配体力场随之确定。在本发明中，可选的配体力场为普适AMBER力场2(Generation Amber Force Field 2)或者CHARMM普适力场(CHARMM General ForceField)等。

在某些具体的实施方式中，所述蛋白质为膜蛋白，所述参数还包括所述膜蛋白的介电常数(membrane dielectric constant)。

作为优选地实施方式，步骤(2)中，所述实验为湿实验；

在本发明的技术方案中，采用湿实验即通过实验室中采用分子、细胞、生理学等试验方法测得的结合自由能而非通过计算机模拟、生物信息学方法得到的结合自由能与步骤(1)预测的数据进行拟合。

又一方面，本发明提供一种基于MM/PB(GB)SA的蛋白质-药物结合自由能预测系统，包括：

探索数据模块：用于收集并存储所述蛋白质及其配体的结构信息数据、所述药物的结构信息数据以及实验测得的所述蛋白质-配体的结合自由能数据；

训练数据模块：用于对所述蛋白质及其配体进行结合自由能预测并将预测得到的结合自由能与实验结果拟合筛选出拟合最好的一组参数；

预测数据模块：以筛选出的拟合最好的一组预测参数分析所述蛋白质和所述药物的结合自由能；

其中，所述结合自由能预测的方法为MM/PBSA和/或MM/GBSA。

优选地，所述训练数据模块包括预处理单元，用于所述蛋白质晶体结构的预处理。

优选地，所述预处理包括加氢、质子化和能量最小化。

又一方面，本发明提供一种计算机可读存储介质，其上存储有计算机程序，该程序被处理器执行上述基于MM/PB(GB)SA的蛋白质-药物结合自由能预测方法。

上述技术方案具有如下优点或者有益效果：

本发明公开了一种基于MM/PB(GB)SA的蛋白质-配体结合自由能计算预测方法，相对于FEP及TI，具有以下优点：(1)本发明提供的预测方法的计算资源消耗较低，且能够节省计算时间；(2)本发明提供的预测方法具有高鲁棒性；(3)本发明提供的预测方法的准确率与自由能微扰技术(FEP，free energy perturbation)及热动力学积分(TI，thermodynamicintegration)相当，且速度要快40～50倍；(4)FEP/TI方法多用于水溶性蛋白，对膜蛋白的效果未知，而本发明提供的方法可以用于预测膜蛋白体系的药物分子结合自由能，且预测结果具有一定的精准性。

附图说明

图1是本发明基于MM/PB(GB)SA的蛋白质-配体结合自由能预测方法的流程图。

图2是本发明以FEP预测水溶性蛋白与其配体结合自由能与实验测得的结合自由能的相关性图。

图3是本发明以MM/PB(GB)SA测试水溶性蛋白与其配体结合自由能与实验测得的结合自由能的相关性图。

图4是本发明以MM/PB(GB)SA测试膜蛋白与配体结合自由能与实验测得的结合自由能的相关性图。

图5是本发明CDK2蛋白测试体系的配体结构图。

图6-1和图6-2是本发明P38蛋白测试体系的配体结构图。

图7是本发明Thrombin蛋白测试体系的配体结构图。

图8是本发明Tyk2蛋白测试体系的配体结构图。

图9是本发明mPGES蛋白测试体系的配体结构图。

图10是本发明GPBAR蛋白测试体系的配体结构图。

图11是本发明OX1蛋白测试体系的配体结构图。

具体实施方式

下述实施例仅仅是本发明的一部分实施例，而不是全部的实施例。因此，以下提供的本发明实施例中的详细描述并非旨在限制要求保护的本发明的范围，而是仅仅表示本发明的选定实施例。基于本发明的实施例，本领域技术人员在没有作出创造性劳动的前提下所获得的所有其他实施例，都属于本发明的保护范围。

在本发明中，若非特指，所有的设备和原料等均可从市场购得或是本行业常用的。下述实施例中的方法，如无特别说明，均为本领域的常规方法。

在一个实施例中，基于MM/PB(GB)SA的蛋白质-药物结合自由能预测方法的流程图如图1所示。

在下文中，结合具体操作对其进行描述。

下述实施中的测试方法和测试对象：水溶性体系中，采用自由能微扰技术FEP测试了4个蛋白，包括：CDK2、P38、Thrombin、Tyk2蛋白；膜蛋白对象包括3个：mPGES、GPBAR、OX1蛋白。

测试体系的构建：在水溶液体系中，首先对已知的蛋白质晶体的三维结构进行加氢、预测氨基酸的质子化状态、能量最小化等处理；然后将其配体分子和所要预测的药物分子通过限制性分子对接技术获得其与蛋白质的结合模式；然后再通过分子动力学工具Gromacs构建12×12×12埃

预测已知的蛋白质药物靶标与配体的结合自由能：上述测试体系建立后，通过分子动力学模拟(molecular dynamics simulations，MD)模拟相关蛋白质药物靶标与配体分子在5ns时间尺度下的多重动态信息；通过MM/PB(GB)SA的热力学循环方程来预测蛋白质药物靶标与配体分子结合自由能，如下所示：

ΔG

＝ΔE

＝ΔE

上式中，ΔG

在通过上式预测结合自由能过程中，对计算过程中的不同参数及获取该参数的方法进行系统化测试，包括：由不同方法计算得到的配体的电荷数、蛋白质药物靶标的分子力场、药物配体的分子力场、是否添加额外的单点卤素模型处理卤素、对膜蛋白体系采用不同的隐式GB模型等。

在一个实施例中，测试了不同的配体电荷计算方法得到的电荷数与不同的蛋白质的分子力场对上述水溶性蛋白质以及膜蛋白与其配体分子结合自由能的影响，其中，配体电荷计算方法包括：RESP-DFT、RESP-HF、AM1-BCC、CGenFF；分子力场包括：CHARMM36m(简写为CHARMM)、Amber FF99SB(简写为FF99SB)。以上述参数和计算方法对蛋白质及其配体预测结合自由能，与实验测得的结合自由能一一拟合后得到的相关系数(obs.R)以及平均绝对误差(MAE)见表1(表中，obs.R后括号中的higher is better意为该指标的数值越大越好，MAE后lower is better的意为该指标的数值越小越好，下述表格同理)。

在一个实施例中，测试了是否添加额外的单点卤素模型处理卤素计算配体的电荷数和不同的蛋白质分子力场对上述水溶性蛋白以及膜蛋白与其配体分子结合自由能的影响，其中，添加单点卤素处理模型的配体电荷计算方法包括RESP_DFT_EP和RESP_HF_EP，未添加单点卤素处理模型的的是RESP_DFT和RESP_HF(其中EP意为额外单点)；分子力场包括：CHARMM36m(简写为CHARMM)、Amber FF99SB(简写为FF99SB)。以上述参数和计算方法对蛋白质及其配体预测结合自由能，与实验测得的结合自由能一一拟合后得到的相关系数(obs.R)以及平均绝对误差(MAE)见表2。

在一个实施例中，测试了不同的隐式GB模型和不同的蛋白质分子力场对上述膜蛋白与其配体分子结合自由能的影响，其中，隐式GB模型包括：GBHCT(igb＝1)、GBOBC(igb＝2)、GBOBC2(igb＝5)、GBneck(igb＝7)、GBneck2(igb＝8)；分子力场包括：CHARMM36m(简写为CHARMM)、Amber FF99SB(简写为FF99SB)，所用到的配体电荷数计算方法包括：RESP_DFT、RESP_HF、AM1_BCC、CGenFF。以上述参数和计算方法对蛋白质及其配体预测结合自由能，与实验测得的结合自由能一一拟合后得到的相关系数(obs.R)以及平均绝对误差(MAE)见表3。

在一个实施例中，测试了在考虑膜蛋白的介电常数(membrane dielectricconstant)的影响下以及不同的蛋白质分子力场对上述膜蛋白与其配体分子结合自由能的影响，其中，膜蛋白的介电常数(membrane dielectric constant)分别设置为emem＝1～9；分子力场包括：CHARMM36m(简写为CHARMM)、Amber FF99SB(简写为FF99SB)，所用到的配体电荷数计算方法包括：RESP_DFT、RESP_HF、AM1_BCC、CGenFF。以上述参数和计算方法对蛋白质及其配体预测结合自由能，与实验测得的结合自由能一一拟合后得到的相关系数(obs.R)以及平均绝对误差(MAE)见表4(其中，inp为非极性溶剂化自由能的计算优化方法)。

在一个实施例中，上述四个水溶性蛋白体系中，以自由能微扰技术(FEP)方法测得的蛋白与其配体的结合自由能与湿实验测得的结合自由能的相关性见图2。

在一个实施例中，上述四个水溶性蛋白体系中，以MM/PB(GB)SA测试蛋白与其配体的结合自由能与实验测得的结合自由能的相关性图见图3。

在一个实施例中，上述三个膜蛋白体系中，以MM/PB(GB)SA测试蛋白与其配体的结合自由能与实验测得的结合自由能的相关性图见图4。

在一个实施例中，以分子对接(docking)、表1-表4中优化后的MM/PB(GB)SA方法与自由能微扰技术(FEP)方法测得的蛋白质与配体结合自由能与通过湿实验测得的结合自由能的相关系数见表5(加粗的为MM/GBSA、MM/PBSA中的较佳数据)：表5中，数值的绝对值越大，预测的精度越高。

表5

因此，本发明提供的预测方法，系统性优化参数后的MM/PB(GB)SA结果普遍优于现有“黄金标准”的自由能微扰技术FEP。

此外，由于目前的FEP/TI技术多采用OPLS分子力场，该力场不能很好的重现真实细胞膜的物理性质。甚至在长时间尺度的分子动力学模拟过程中细胞膜会发生破裂现象，因此基于该力场的FEP/TI不适合膜蛋白体系的药物分子结合自由能预测。由于本发明中采用的CHARMM36m力场和AMBER99SB力场均有非常可靠的膜蛋白力场，可以很好的适用于膜蛋白体系的药物分子结合自由能预测。从预测结果可以看出，三个不同测试体系对于膜蛋白体系的预测结果均与实验值高度一致。

在速度方面，对FEP/TI计算而言，从公开报道来看，一个GTX2080ti显卡(GPU)上，FEP/TI一天只能预测0.4-0.5个分子。而优化后的MM/PB(GB)SA一天可以预测20个分子。速度比FEP/TI的方法要快40-50倍。

在一个实施例中，本发明提供一种基于MM/PB(GB)SA的蛋白质-药物结合自由能预测系统，包括：

探索数据模块：用于收集并存储蛋白质及其配体的结构信息数据、药物的结构信息数据以及实验测得的蛋白质-配体的结合自由能数据；

训练数据模块：用于对蛋白质及其配体进行结合自由能预测并将预测得到的结合自由能与实验结果拟合筛选出拟合最好的一组参数；

预测数据模块：以筛选出的拟合最好的一组预测参数分析蛋白质和药物的结合自由能；

其中，结合自由能预测的方法为MM/PBSA和/或MM/GBSA。

优选地，训练数据模块包括预处理单元，用于蛋白质晶体结构的预处理。

优选地，预处理包括加氢、质子化和能量最小化。

以上所述仅是本发明的优选实施方式，应当指出：对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。

/>