WO2020147668A1 - 预测药-靶结合强度的绝对结合自由能微扰方法 - Google Patents

预测药-靶结合强度的绝对结合自由能微扰方法 Download PDF

Info

Publication number
WO2020147668A1
WO2020147668A1 PCT/CN2020/071663 CN2020071663W WO2020147668A1 WO 2020147668 A1 WO2020147668 A1 WO 2020147668A1 CN 2020071663 W CN2020071663 W CN 2020071663W WO 2020147668 A1 WO2020147668 A1 WO 2020147668A1
Authority
WO
WIPO (PCT)
Prior art keywords
probability distribution
free energy
drug
calculation
window
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.)
Ceased
Application number
PCT/CN2020/071663
Other languages
English (en)
French (fr)
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
Priority claimed from CN201910045003.XA external-priority patent/CN109859806B/zh
Priority claimed from CN201911412608.4A external-priority patent/CN111161810B/zh
Application filed by Sun Yat Sen University filed Critical Sun Yat Sen University
Publication of WO2020147668A1 publication Critical patent/WO2020147668A1/zh
Anticipated expiration legal-status Critical
Ceased legal-status Critical Current

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

Definitions

  • the present invention relates to the technical field of drug research and development, and more specifically, to an absolute binding free energy perturbation method for predicting drug-target binding strength.
  • the present invention provides a perturbation method of absolute binding free energy for predicting the strength of drug-target binding.
  • the technical scheme of the present invention is as follows:
  • the absolute binding free energy perturbation method for predicting drug-target binding strength includes the following steps:
  • S1 Perform conventional molecular dynamics simulation on the receptor protein-ligand system, namely the Rec-Lig system, and obtain the distance, angle, and dihedral angle data required to add constraints based on the trajectory obtained by the simulation;
  • step S1 specifically includes the following steps:
  • step S2 specifically includes the following steps:
  • step S21 According to the prepared force field parameters of the Rec-Lig system and the 6 atoms obtained in step S16, add a length, two included angles, and three dihedral angle constraints to obtain the calculation Molecular dynamics simulation parameters of the required window;
  • k B is Boltzmann's constant
  • V is the volume of the system
  • raA is the distance between two atoms for the restriction of distance harmonic potential
  • ⁇ a and ⁇ A are the angles between two three atoms for Simple harmonic potential constraint of the included angle
  • many k values are used to constrain a distance raA , two included angles ⁇ a , ⁇ A , three dihedral angles
  • the simple harmonic force constant.
  • any molecular dynamics simulation software including AMBER, Gromacs, Namd is used for calculation.
  • step S5 specifically includes the following steps:
  • S51 Perform probability distribution statistics for the potential energy difference between each adjacent state to obtain the probability distribution P( ⁇ U) data ;
  • P( ⁇ U) fit is the probability distribution after Gaussian function fitting
  • c i is the parameter of the i-th Gaussian function
  • ⁇ i is the average of the i-th Gaussian function
  • ⁇ i is the i-th Gaussian function Standard deviation
  • ⁇ A i,i+1 is the binding free energy difference between two adjacent windows i and i+1
  • ⁇ U i,i+1 is the potential energy difference newly generated based on the fitted probability distribution
  • ⁇ 1/k B T
  • k B are Boltzmann's constants
  • n i and n i+1 are the sample sizes included in the samples of two adjacent states i and i+1 respectively
  • C is the implicit variable between the equations
  • step S6 specifically includes the following steps:
  • P ( ⁇ U) fit through the constraint of the probability distribution function fitting the smoothed probability distribution, where a i, b i, c i , d i, h i, ⁇ i, ⁇ i is the i th "probability constraint “Distribution function” parameter, n is the number of selected fitting functions, n 1 and n 2 are correction factors of two positive integers.
  • step S7 specifically includes the following steps:
  • step S72 Based on the new potential energy difference ⁇ U′, use the calculation method of step S5 to recalculate Calculated value;
  • step S73 obtained according to step S24, step S5 and step S72 with Calculate the drug-target binding free energy ⁇ A binding .
  • the calculation process of the drug-target binding free energy ⁇ A binding is specifically:
  • the invention provides an absolute binding free energy perturbation method for predicting drug-target binding strength, which is based on the fitting of the constraint probability distribution function to the probability distribution of the constraint addition process, and the probability distribution obtained by statistics of the windows of the remaining thermodynamic cycle steps,
  • the linear combination of multiple Gaussian functions is used to fit the probability distribution of the remaining windows, and the potential energy difference is regenerated according to the fitted probability distribution, and finally the drug-target absolute binding free energy is obtained by statistics; by perturbing the intermediate state of the free energy
  • the probability distribution corresponding to the kinetic simulation trajectory is fitted to make the probability distribution of sampling more reasonable, which greatly improves the convergence and accuracy of the energy calculation results.
  • Figure 1 is a schematic diagram of the thermodynamic cycle involved when the mass gap between different small molecules of the same target is relatively small without adding constraints;
  • Figure 2 is a schematic diagram of the probability distribution P( ⁇ U) of the conversion of ⁇ U obtained by sampling
  • Figure 3 is a schematic diagram of the result of fitting the probability distribution P( ⁇ U) obtained by sampling by a series of Gaussian functions
  • Figure 4 is a schematic diagram of the accuracy of the method before and after fitting the sampling probability distribution by the Gaussian function
  • Figure 5 is a schematic diagram of the relationship between exp(- ⁇ U)P 0 ( ⁇ U) and P 0 ( ⁇ U) in the process of calculating energy;
  • FIG. 6 is a schematic diagram of accuracy comparison between the absolute binding free energy calculation method enhanced by the Gaussian function in the present invention and the relative binding free energy calculation method of the Schrodinger formula;
  • Figure 7 is a schematic flow diagram of the method
  • Figure 8 is a thermodynamic cycle diagram involved in calculating absolute binding free energy
  • Figure 9 is a schematic diagram of the result of fitting the probability obtained by the state sampling of the constraint process by the constraint probability distribution function.
  • thermodynamic cycle shown in Figure 1 is used without considering the binding force, where Rec stands for Receptor, which is drug-target receptor; Lig stands for Ligand, which is drug small molecule ligand ; Dumm stands for Dummy, that is, virtual atoms. All virtual atoms have no interaction with the surrounding environment.
  • the binding free energy ⁇ A binding of the Lig system and the Rec-Lig system is obtained by calculation. According to the thermodynamic cycle shown in Figure 1, the with After the size, the binding free energy ⁇ A binding can be obtained.
  • S11 Prepare the structure of the receptor-ligand complex, and use the molecular docking method to dock 16 small drug molecule ligands into the binding pocket of the large molecule CDK2 receptor;
  • ⁇ ij is the LJ potential energy parameter between atom i and atom j
  • r ij is the distance between atom i and atom j
  • the van der Waals interaction between atoms is weakened by reducing ⁇ ij ;
  • the van der Waals interaction is annihilated, it becomes more difficult to sample these unreal states. Therefore, when the van der Waals effect is reduced to a relatively small level, the different states need to be closer to improve the convergence of sampling. Therefore, the ⁇ of each atom decreases in the form of (1-0.2i) 6 . After five steps of modification of van der Waals interaction in this way, all the effects of small molecules with the outside world disappeared, achieving the effect of annihilating small molecules. Save all modified ⁇ parameter files as independent force field files as the parameters of the intermediate perturbation state;
  • step S2 Use the molecular dynamics simulation program AMBER to perform molecular dynamics simulation on the Lig system and the Rec-Lig system using the original force field parameter file and the force field file of the intermediate perturbation state in step S1.
  • the SHAKE function cannot be used to restrict the system during simulation, because the characteristics of the SHAKE function will make the intermediate state unstable, so the simulation step size is chosen as 1fs.
  • the periodic boundary conditions of the system need to be calculated using the PME method for the long-range electrostatic interaction, and the cut-off value of 8A is selected for the short-range van der Waals interaction; for each state of each system, it is performed under the condition of 300K 4ns molecular dynamics simulation. In the last 2ns simulation process, one frame of constellation is extracted every 100fs, and the trajectory of each state contains 20,000 frames of constellation;
  • the system average of the potential energy difference between the state i+1 and the state i is solved, and the result is the forward perturbative potential energy difference; solve the state i-1 state and state i
  • the system average of the potential energy difference of, the backward perturbation potential energy difference is solved from this; the results of the forward perturbation potential energy difference and the backward perturbation potential energy difference are combined, and the BAR method is used for calculation, which effectively improves the final energy calculation Accuracy.
  • S4 Perform probability distribution statistics for the potential energy difference between each adjacent state, and use the linear combination of multiple Gaussian functions to fit to obtain a smooth probability distribution, which is specifically:
  • step S3 the trajectory of each state contains 20,000 conformations. Therefore, the trajectory of each state is calculated to obtain 20000 data points of ⁇ U; as shown in Figure 2, the specifics are:
  • P( ⁇ U) fit is the probability distribution after Gaussian function fitting
  • c i is the parameter of the i-th Gaussian function
  • ⁇ i is the average of the i-th Gaussian function
  • ⁇ i is the i-th Gaussian function Standard deviation
  • c i 1, 2, 3, 4, 5
  • P( ⁇ U) data is the probability distribution obtained by dynamic sampling
  • P( ⁇ U) fit is the probability distribution after Gaussian function fitting
  • c i is the parameter of the i-th Gaussian function
  • ⁇ i is the i-th Gaussian function
  • the average value of the function, ⁇ i is the standard deviation of the i-th Gaussian function
  • (c i , ⁇ i , ⁇ i ) t is the parameter value after t iterations; through continuous iteration until the tth iteration (c i , ⁇ i , ⁇ i ) t is the fitted function The most important parameters. After fitting, the probability distribution of the result is shown in Figure 3.
  • ⁇ A i,i+1 is the binding free energy difference between two adjacent states i and i+1
  • ⁇ U i,i+1 is the potential energy difference newly generated based on the fitted probability distribution
  • ⁇ 1/k B T
  • k B are Boltzmann's constants
  • n i and n i+1 are the sample sizes included in the samples of two adjacent states i and i+1 respectively
  • C is the implicit variable between the equations
  • the 16 small molecules of the CDK2 system are calculated.
  • the structures of these 16 small molecules are known, and the IC 50 can be determined experimentally, combined with the calculation formula of the experimental value of the free energy for:
  • Fig. 4 The result obtained is shown in Fig. 4, which can be obtained by comparison.
  • the absolute binding degree of freedom calculated by the method of the present invention has a very high correlation with the experimental value, so the result obtained by the method has high accuracy.
  • the free energy perturbation method enhanced by the Gaussian function for predicting the absolute binding free energy between ligands and receptors has high accuracy, and the existing free energy perturbation method can only be used to calculate the relative With the free energy of binding, the present invention can calculate the absolute free energy of binding, so the method of the present invention can be more widely used in the work of rational drug design.
  • the method of fitting sampling probability distribution P( ⁇ U) data to improve the accuracy and convergence of free energy perturbation is as follows:
  • ⁇ A - ⁇ ln ⁇ exp(- ⁇ U)P 0 ( ⁇ U)d ⁇ U;
  • ⁇ is a constant
  • ⁇ A represents the free energy difference between the two states
  • the sampled P( ⁇ U) needs to be multiplied by exp(- ⁇ U), and the result is shown in Fig. 5.
  • the part of P( ⁇ U) in the smaller ⁇ U is amplified by the product of exp(- ⁇ U), which causes the sampling accuracy of this part to have a very large impact on the calculation of energy; and this part is due to the data points
  • the probability of occurrence is small, which often results in insufficient sampling. Therefore, it will lead to insufficient convergence and accuracy when faced with large perturbations.
  • the fitting is performed with the sum of a series of Gaussian functions.
  • the fitting function fully takes into account the overall P( ⁇ U) In nature, the part where P( ⁇ U) is at a smaller ⁇ U is also optimized by extrapolation, and the accuracy and convergence of the optimized energy calculation are greatly improved. As shown, the product is integrated, and the curve area obtained is proportional to the free
  • the accuracy of absolute free energy calculation using free energy perturbation by the method of the present invention is greatly improved.
  • the existing free energy perturbation method does not fit the probability distribution P( ⁇ U), and the linear regression results of the experimental and calculated values are shown in Figure 4.
  • the regression coefficient between the experiment and the calculated value is 0.40; and after the probability distribution is fitted and calculated using the method of the present invention, the regression coefficient between the obtained experiment and the calculated value is 0.88.
  • the comparison shows that the invention patent Compared with the existing free energy perturbation method, the free energy perturbation method has greatly improved the accuracy when calculating the absolute binding free energy of the drug and the target.
  • the free energy perturbation method enhanced by the Gaussian function is used to calculate the absolute free energy, and the results are compared with the existing binding free energy calculation methods.
  • the comparison results are shown in Figure 6. It can be seen from the figure that the accuracy of the absolute free energy calculation using the free energy perturbation method enhanced by the Gaussian function is equivalent to the relative binding free energy calculation method, but Absolute free energy calculation does not require a known reference drug molecule, and its application range is far greater than the existing calculation methods of relative binding free energy.
  • S1 Perform conventional molecular dynamics simulation on the receptor protein-ligand system, namely the Rec-Lig system, and obtain the distance, angle, and dihedral angle data required to add constraints based on the trajectory obtained by the simulation;
  • the present invention provides a free energy perturbation method based on optimization of the constraint probability distribution function, which realizes the fully automatic acquisition of the composite topology file required for adding constraints, and the constraint probability distribution function is used to fit the experimental results
  • the probability distribution makes the free energy change free value in the free energy perturbation (FEP) constrained window more accurate.
  • FEP free energy perturbation
  • the probability distribution of the constraint step is fitted based on the automatic constraint addition method and constraint probability distribution function provided by the present invention, based on human phosphodiesterase 4D and its small molecule inhibitor roflumilast
  • the complex system (PDB code is 1XOQ) and human BRD4 protein and its small molecule ligand alprazolam (PDB code 3U5J) are taken as examples, and the constraint window is changed from 10 used in the existing literature (Chem.Sci., 2016, 7, 207-218; J. Am. Chem. Soc. 2017, 139, 946-957, etc.) reduced to one, calculated by comparing and examining the 10-step plan Compared with the 1-step program, the test results obtained are shown in Table 1.
  • step S1 specifically includes the following steps:
  • step S2 specifically includes the following steps:
  • step S21 According to the prepared force field parameters of the Rec-Lig system and the 6 atoms obtained in step S16, add a length, two included angles, and three dihedral angle constraints to obtain the calculation Molecular dynamics simulation parameters of the required window;
  • k B is Boltzmann's constant
  • V is the volume of the system
  • raA is the distance between two atoms for the restriction of distance harmonic potential
  • ⁇ a and ⁇ A are the angles between two three atoms for Simple harmonic potential constraint of the included angle
  • many k values are used to constrain a distance raA , two included angles ⁇ a , ⁇ A , three dihedral angles
  • the simple harmonic force constant.
  • step S5 is specifically:
  • ⁇ A i,i+1 is the binding free energy difference between two adjacent windows i and i+1
  • ⁇ U i,i+1 is the potential energy difference newly generated based on the fitted probability distribution
  • ⁇ 1/k B T
  • k B are Boltzmann's constants
  • n i and n i+1 are the sample sizes included in the samples of two adjacent states i and i+1 respectively
  • C is the implicit variable between the equations
  • step S6 specifically includes the following steps:
  • FIG. 9 is a schematic diagram of the result of fitting the probability distribution P( ⁇ U) obtained by sampling the state of the constraint adding process by a series of constraint probability distribution functions.
  • the gray point in the figure is the probability distribution point obtained by sampling
  • the black line is the distribution curve obtained by the final fitting
  • the gray solid line and gray dashed line are the first and last two items in the formula.
  • step S7 specifically includes the following steps:
  • step S72 Based on the new potential energy difference ⁇ U′, use the calculation method of step S5 to recalculate Calculated value;
  • the calculation process of the drug-target binding free energy ⁇ A binding is specifically:
  • Example 1 show that the Gaussian function fitting can increase with The convergence of the calculation.
  • Example 2 show that the method of automatically adding constraints combined with the constraint probability distribution function for fitting can make The term energy settlement result achieves very good convergence.
  • two small drug molecules targeting BRD4 protein targets are taken as examples.
  • the constraint scheme of automatically adding constraints and setting only one disturbance window during the constraint process is adopted.
  • the constraint probability distribution function of the present invention is used to fit its probability distribution scatter plot.
  • the absolute binding free energy ⁇ A binding is calculated. Specific steps are as follows:
  • thermodynamic cycle as shown in Figure 8, construct The molecular dynamics simulation parameters required for the step, where the This step adds constraints using the distance, included angle, and dihedral angle data obtained in step S1, and The free energy in this step is calculated analytically;
  • the number of windows in the constraint step can be greatly reduced on the premise of ensuring the calculation accuracy, so that the calculation
  • the number of windows is reduced from 10 to 1, which also greatly reduces the amount of calculation, which is only 1/10 of the previous amount.
  • the above examples also illustrate the accuracy of the results.

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

一种预测药-靶结合强度的绝对结合自由能微扰方法，基于约束概率分布函数针对添加约束过程的概率分布的拟合，对其余热力学循环步骤的窗口统计得到的概率分布，利用多个高斯函数的线性组合对其余窗口的概率分布进行拟合，根据拟合后的概率分布重新生成势能差，最终统计得出药物-靶标绝对结合自由能。通过对自由能微扰中间状态的动力学模拟轨迹所对应的概率分布进行拟合，使得采样的概率分布更加合理，大幅提升了能量计算结果的收敛性及准确性。

Description

预测药-靶结合强度的绝对结合自由能微扰方法 技术领域
本发明涉及药物研发技术领域，更具体的，涉及一种预测药-靶结合强度的绝对结合自由能微扰方法。
背景技术
药物产生药效的原因几乎都与药物和靶标大分子的相互结合密切相关。因此，药物-靶标结合强度的自由能预测技术在药物设计中扮演着极其重要的角色(Q.Rev.Biophys.2012,45,1-25.；Curr.Opin.Struct.Biol.2011,21,150-160.等)。若有一种方法能够准确地预测结合自由能，那么新药研发的周期与成本将大大降低。
尽管目前存在一些较为快速的结合自由能预测方法，例如对接打分函数，MM/PBSA(分子力学/泊松玻尔兹曼表面，药物设计中常用的一种能量计算方法)、MM/GBSA(分子力学/广义伯恩表面)等方法。但由于这些方法多引入了很多的经验参数，或者理论基础存在一定的缺陷，导致它们的准确性不高。自由能微扰(FEP)作为一种基于统计力学的自由能计算方法，理论上具有更高的准确度。如著名药物设计软件公司薛定谔Schrodinger在2015年发明了一种基于自由能微扰的相对结合自由能预测方法(J.Am.Chem.Soc.2015,137,2695-2703.)，该方法能预测对现有药物分子只有微小结构修饰的同骨架化合物与靶标的相对结合自由能，且往往这类预测方法需要参考已知活性的小分子结构。新药研发中遇到的大多数问题情况，如虚拟筛选，骨架跃迁和全新药物设计时，需要进行绝对结合自由能的预测。
发明内容
本发明为克服上述现有技术存在无法进行绝对结合自由能预测的技术问题，提供一种预测药-靶结合强度的绝对结合自由能微扰方法。
为解决上述技术问题，本发明的技术方案如下：
预测药-靶结合强度的绝对结合自由能微扰方法，包括以下步骤：
S1：对受体蛋白-配体体系，即Rec-Lig体系进行常规的分子动力学模拟，并基于模拟获得的轨迹得到添加约束所需的距离、夹角、二面角数据；
S2：构建热力学循环，构建
Figure PCTCN2020071663-appb-000001
和
Figure PCTCN2020071663-appb-000002
根据添加约束所需的距离、夹角、二面角数据得到，并对
Figure PCTCN2020071663-appb-000003
的自由能进行解析计算；
其中：
Figure PCTCN2020071663-appb-000004
表示将溶液中药物小分子全部微扰至湮灭状态的能量；
Figure PCTCN2020071663-appb-000005
表示将药物-靶标复合物中药物小分子全部微扰至湮灭状态的能量；
Figure PCTCN2020071663-appb-000006
表示对药物-靶标复合物中药物小分子添加简谐势约束的能量变化；
S3：利用分子动力学模拟程序，对
Figure PCTCN2020071663-appb-000007
步骤中的每一个窗口进行分子动力学模拟采样，得到
Figure PCTCN2020071663-appb-000008
Figure PCTCN2020071663-appb-000009
步骤中每一个窗口的分子动力学轨迹；
S4：根据
Figure PCTCN2020071663-appb-000010
步骤中的每一个窗口的分子动力学轨迹，通过相邻窗口的力场文件进行计算，得到每个相邻窗口之间的势能差ΔU；
S5：对于
Figure PCTCN2020071663-appb-000011
这两步，根据对应窗口的ΔU的概率分布，，使用高斯函数进行拟合，并计算得到
Figure PCTCN2020071663-appb-000012
步骤的自由能；
S6：对于
Figure PCTCN2020071663-appb-000013
这一步每个相邻窗口之间的势能差ΔU进行概率分布统计，利用约束概率分布函数进行拟合，得到平滑的概率分布，并根据拟合结果检查收敛性；
S7：根据拟合后的概率分布，重新生成
Figure PCTCN2020071663-appb-000014
这一步中每一窗口的势能差，通过计算得到
Figure PCTCN2020071663-appb-000015
的能量，并根据之前步骤得到的
Figure PCTCN2020071663-appb-000016
Figure PCTCN2020071663-appb-000017
以及
Figure PCTCN2020071663-appb-000018
计算得到药-靶结合自由能ΔA binding。
其中，所述步骤S1具体包括以下步骤：
S11：准备受体蛋白Rec-配体Lig复合物坐标文件；
S12：准备小分子配体的电荷，并对配体指定力场；
S13：将受体-配体复合物中的受体指定力场，并将该复合物在水箱中进行浸泡，准备好复合物进行分子模拟的初始体系，记为Rec-Lig体系，并到处力场文件作为Rec-Lig体系初始状态参数；
S14：将单独小分子配体浸泡在水箱中，准备好单独的小分子的分子模拟初始体系，记为Lig体系，并到处力场文件作为Lig体系初始状态参数；
S15：根据Rec-Lig体系初始状态的力场参数，进行短时间的分子动力学模拟，得到运动轨迹；
S16：基于模拟得到的轨迹，统计与小分子配体相对位置变化较小的受体蛋白的氨基酸残基，取该残基中与小分子配体相对距离最稳定的三个原子作为添加约束所需的属于受体蛋白的三个原子，选取小分子配体中最稳定的三个原子作为添加约束所需的属于配体小分子的三个原子，得到添加约束所需的距离、夹角、二面角数据。
其中，所述步骤S2具体包括以下步骤：
S21：根据准备好的Rec-Lig体系的力场参数，根据步骤S16得到的6个原子，添加一个长度、两个夹角、三个二面角约束，得到计算
Figure PCTCN2020071663-appb-000019
所需窗口的分子动力学模拟参数；
S22：进一步将Rec-Lig体系中的小分子的电荷以及范德华参数依次分别减小至0，得到计算
Figure PCTCN2020071663-appb-000020
所需窗口的分子动力学模拟参数；
S23：对于Rec-Lig体系中小分子的电荷以及范德华参数依次分别减小至0，得到计算
Figure PCTCN2020071663-appb-000021
所需窗口的分子动力学模拟参数；
S24：对
Figure PCTCN2020071663-appb-000022
的自由能进行解析计算。
其中，在所述步骤S24中，所述
Figure PCTCN2020071663-appb-000023
自由能的计算公式具体为：
Figure PCTCN2020071663-appb-000024
其中，k B为玻尔兹曼常数，V为系统的体积，r aA为一个两原子间的距离用于距离简谐势约束，θ a和θ A为两个三原子间的夹角用于夹角简谐势约束，诸多k值为用于约束一个距离r aA、两个夹角θ a，θ A、三个二面角
Figure PCTCN2020071663-appb-000025
的简谐势力常数。
其中，在所述步骤S3或S4中，采用包括AMBER、Gromacs、Namd在内的任意分子动力学模拟软件进行计算。
其中，所述步骤S5具体包括以下步骤：
S51：针对得到每个相邻状态之间的势能差进行概率分布统计，得到概率分布P(ΔU) data；
S52：将概率分布P(ΔU) data利用多个高斯函数的线性组合进行拟合，得到平滑的概率分布，具体计算公式为：
Figure PCTCN2020071663-appb-000026
其中，P(ΔU) fit为经高斯函数拟合后的概率分布，c i为第i个高斯函数的参数，μ i为第i个高斯函数的平均值，σ i为第i个高斯函数的标准差；
S53：根据得到的P(ΔU) fit，采用Benneet Acceptance Ratio计算方法，即BAR计算方法进行
Figure PCTCN2020071663-appb-000027
自由能的计算，具体计算公式为：
Figure PCTCN2020071663-appb-000028
Figure PCTCN2020071663-appb-000029
其中，ΔA i，i+1为两相邻窗口i与i+1的结合自由能差，ΔU i，i+1为基于拟合后的概率分布新生成的势能差，β＝1/k BT，k B为玻尔兹曼常数，n i与n i+1分别为两相邻状态i与i+1采样所包含的样本量，C为方程组之间的隐含变量；
S54：将
Figure PCTCN2020071663-appb-000030
和
Figure PCTCN2020071663-appb-000031
对应的计算窗口得到的ΔA i，i+1进行加和，得到
Figure PCTCN2020071663-appb-000032
和
Figure PCTCN2020071663-appb-000033
的数值。
其中，所述步骤S6具体包括以下步骤：
S61：对计算
Figure PCTCN2020071663-appb-000034
的窗口中每个窗口之间的势能差进行概率分布统计，得到概率分布P(ΔU) data；
S62：基于约束概率分布函数对概率分布P(ΔU) data进行拟合，具体表达式为：
Figure PCTCN2020071663-appb-000035
其中，P(ΔU) fit为经过约束概率分布函数拟合后的平滑的概率分布，其中a i，b i，c i，d i，h i，μ i，σ i为第i个“约束概率分布函数”的参数，n为选取拟合函数的数目，n 1和n 2为两个正整数的修正因子。
其中，所述步骤S7具体包括以下步骤：
S71：根据概率分布P(ΔU) fit，重新生成每一状态的势能差ΔU′；
S72：基于新的势能差ΔU′，利用步骤S5的计算方法重新计算得到
Figure PCTCN2020071663-appb-000036
的计算数值；
S73：根据步骤S24、步骤S5和步骤S72得到的
Figure PCTCN2020071663-appb-000037
Figure PCTCN2020071663-appb-000038
和
Figure PCTCN2020071663-appb-000039
计算出药-靶结合自由能ΔA binding。
其中，在所述步骤S73中，所述药-靶结合自由能ΔA binding的计算过程具体为：
Figure PCTCN2020071663-appb-000040
与现有技术相比，本发明技术方案的有益效果是：
本发明提供的一种预测药-靶结合强度的绝对结合自由能微扰方法，基于约束概率分布函数针对添加约束过程的概率分布的拟合，对其余热力学循环步骤的窗口统计得到的概率分布，利用多个高斯函数的线性组合对其余窗口的概率分布进行拟合，根据拟合后的概率分布重新生成势能差，最终统计得出药物-靶标绝对结合自由能；通过对自由能微扰中间状态的动力学模拟轨迹所对应的概率分布进行拟合，使得采样的概率分布更加合理，大幅提升了能量计算结果的收敛性及准确性。
附图说明
图1为在同一靶标不同小分子之间质量差距比较小时，不添加约束所涉及的热力学循环示意图；
图2为采样所得ΔU被转化的概率分布P(ΔU)示意图；
图3为采样所得概率分布P(ΔU)被一系列高斯函数拟合的结果示意图；
图4为经高斯函数拟合采样概率分布前后方法准确度的差异示意图；
图5为计算能量过程中，exp(-βΔU)P 0(ΔU)与P 0(ΔU)的关系示意图；
图6为本发明中高斯函数增强的绝对结合自由能计算方法与薛定谔公式的相对结合自由能计算方法的准确度对比示意图；
图7为本方法流程示意图；
图8为计算绝对结合自由能所涉及的热力学循环图；
图9为添加约束过程的状态采样所得概率被约束概率分布函数拟合的结果示意图。
具体实施方式
附图仅用于示例性说明，不能理解为对本专利的限制；
为了更好说明本实施例，附图某些部件会有省略、放大或缩小，并不代表实际产品的尺寸；
对于本领域技术人员来说，附图中某些公知结构及其说明可能省略是可以理解的。
下面结合附图和实施例对本发明的技术方案做进一步的说明。
实施例1
首先利用薛定谔公司测试相对结合自由能的靶标与配体分子(J.Am.Chem.Soc.2015,137,2695-2703.)对本专利报道的方法进行测试。
由于小分子间差异较小，因此利用不考虑约束力的如图1所示的较为简单的热力学循环，其中Rec代表Receptor，即药物-标靶受体；Lig代表Ligand，即药物小分子配体；Dumm代表Dummy，即虚原子，所有的虚原子与周围的环境没有任何的相互作用。通过计算得到Lig体系及Rec-Lig体系的结合自由能ΔA binding，依据如图1所示的热力学循环，在分别求得
Figure PCTCN2020071663-appb-000041
和
Figure PCTCN2020071663-appb-000042
的大小后，即可得到结合自由能ΔA binding。
进一步的，针对CDK2药物靶标的16个药物小分子为例，使用不考虑约束力的如图1所示的较为简单的热力学循环的具体步骤如下：
S1：构建药物小分子配体Lig体系以及药物-靶标的Rec-Lig体系，为中间微扰状态参数准备，具体为：
S11：准备受体-配体复合物的结构，利用分子对接的方法将16个药物小分子配体对接进大分子CDK2受体的结合口袋；
S12：将小分子配体通过Gaussian软件在HF/6-31G*水平下计算ESP电荷，使用Antechamber拟合出RESP电荷，并指定GAFF力场参数；
S13：将受体-配体复合物中的受体指定AMBER03力场对于整个复合物添加
Figure PCTCN2020071663-appb-000043
的去角八面体的水箱，添加抗衡离子使整个体系保持电中性，将最终产生的体系的坐标与初始的力场参数均保存为力场文件，作为Rec-Lig体系自由能微扰的初始状态参数；
S14：将单独的小分子配体添加
Figure PCTCN2020071663-appb-000044
的去角八面体的水箱，将最终产生的体系的坐标与初始的力场参数均保存为力场文件，作为Lig体系进行自由能微扰的初始状态的参数；
S15：将Rec-Lig体系以及Lig体系中的小分子电荷进行湮灭。对于这两个体系，分别通过5步，将小分子自中所有原子的电荷都减小至0。即每一步将小分子中每个原子的电荷修改为原来的0.8、0.6、0.4、0.2、0，将每一步修改过电荷之后的参数文件分别保存为独立的力场文件，作为中间微扰状态的参数；
S16：在小分子电荷进行湮灭后，小分子与外界还保留着范德华相互作用。 因此需要对Rec-Lig体系以及Lig体系中小分子的范德华相互作用进行湮灭。分子间的范德华相互作用以Lennard-Jones(LJ)势能计算，计算公式为:
Figure PCTCN2020071663-appb-000045
其中，ε ij为原子i和原子j之间的LJ势能参数，r ij为原子i和原子j之间距离；在本实施例中，通过减小ε ij来减弱原子间的范德华相互作用；由于随着范德华相互作用的湮灭，对这些非真实状态的采样变得越发困难，因此当范德华作用减小到比较小时，需让不同状态之间更加接近以提高采样的收敛性。因此每个原子的ε以(1-0.2i) 6的形式进行减少。按照该方式进行了5步对于范德华相互作用的修改之后，小分子与外界的所有作用全部消失，达到了对小分子湮灭的效果。将所有修改过ε的参数文件保存为独立的力场文件，作为中间微扰状态的参数；
S2：利用分子动力学模拟程序AMBER，对Lig体系及Rec-Lig体系分别用原始的力场参数文件以及步骤S1生活曾的中间微扰状态的力场文件进行分子动力学模拟。在模拟时不可使用SHAKE函数对体系进行限制，因为SHAKE函数的特性会使中间状态不稳定，因此模拟的步长选择为1fs。此外，需要对体系周期性边界条件，对长程静电互相作用作用采用PME的方法进行计算，对短程的范德华相互作用选取8A的截断值；对于每个体系的每个状态都在300K的条件下进行4ns的分子动力学模拟。在最后2ns的模拟过程中，每隔100fs提取一帧的构象，每个状态的轨迹包含了20000帧构象；
S3：将得到的Lig体系及Rec-Lig体系每个状态的分子动力学轨迹，通过相邻状态的力场文件进行单点能计算，得到每个相邻状态之间的势能差ΔU，具体为：
基于第i状态的分子动力学轨迹为基础，求解状态i+1状态与状态i的势能差的系统平均，由此求解出来的为向前微扰势能差；求解状态i-1状态与状态i的势能差的系统平均，由此求解出来的为向后微扰势能差；将向前微扰势能差和向后微扰势能差的结果结合起来，使用BAR方法进行计算，有效提高最终能量计算的准确性。
S4：针对得到每个相邻状态之间的势能差进行概率分布统计，利用多个高斯函数的线性组合进行拟合，得到平滑的概率分布，具体为：
对于所有的势能差ΔU的概率分布进行你和，经过步骤S3后，由于每个状态 的轨迹都包含20000个构象。因此每个状态的轨迹经计算后得到20000个ΔU的数据点；如图2所示，具体为：
S41：先将所得势能差ΔU转换成对应的概率分布P(ΔU) data；
S42：利用多个高斯函数的线性组合对概率分布P(ΔU) data进行拟合，具体计算公式为：
Figure PCTCN2020071663-appb-000046
其中，P(ΔU) fit为经高斯函数拟合后的概率分布，c i为第i个高斯函数的参数，μ i为第i个高斯函数的平均值，σ i为第i个高斯函数的标准差；共包含15个参数c i，μ i，σ i，其中i＝1，2，3，4，5，由此定义费用函数的形式为：
Figure PCTCN2020071663-appb-000047
其中，P(ΔU) data为动力学采样得到的概率分布，P(ΔU) fit为经高斯函数拟合后的概率分布，c i为第i个高斯函数的参数，μ i为第i个高斯函数的平均值，σ i为第i个高斯函数的标准差；利用随机梯度下降法对费用函数进行优化，具体为采用SGD法优化这15个参数c i，μ i，σ i，其中i＝1，2，3，4，5，具体为：
Figure PCTCN2020071663-appb-000048
其中，(c i，μ i，σ i) t为经t次迭代后的参数值；通过不断迭代直至第t次迭代所得(c i，μ i，σ i) t即为拟合后的函数的最有参数。经拟合后，结果的概率分布如图3所示。
S5：根据拟合后的概率分布，重新生成每一状态的势能差，通过计算得到药物与靶标的结合自由能，具体为：
S51：利用经过高斯函数拟合后的概率分布，重新生成500000至1000000个新的势能差ΔU′；
S52：并基于新生成的势能差ΔU′求解方程组，得到相邻状态之间的结合自由能之差ΔA i，i+1，具体方程为：
Figure PCTCN2020071663-appb-000049
Figure PCTCN2020071663-appb-000050
其中，ΔA i，i+1为两相邻状态i与i+1的结合自由能差，ΔU i，i+1为基于拟合后的概率分布新生成的势能差，β＝1/k BT，k B为玻尔兹曼常数，n i与n i+1分别为两相邻状态i与i+1采样所包含的样本量，C为方程组之间的隐含变量；
S53：对Rec-Lig体系所有状态势能差进行加和运算，得到能量
Figure PCTCN2020071663-appb-000051
对Lig体系所有状态能量差进行加和运算，得到能量
Figure PCTCN2020071663-appb-000052
计算得到药物与靶标的结合自由能ΔA binding，具体计算公式为：
Figure PCTCN2020071663-appb-000053
其中，
Figure PCTCN2020071663-appb-000054
表示将溶液中药物小分子全部微扰至湮灭状态的能量；
Figure PCTCN2020071663-appb-000055
表示将药物-靶标复合物中药物小分子全部扰乱至湮灭砖头的能量；其中，湮灭状态为药物小分子的电荷以及范德华力L-J势能为0，周围院子无法感知到药物小分子的存在。
在具体实施过程中，根据步骤S1-S5，对于CDK2体系的16个小分子进行了计算，这16个小分子的结构已知，IC 50可以由实验测定，结合自由能的实验值的计算公式为：
ΔG＝-RT ln IC 50；
得到的结果如图4所示，对比可得到，本发明所述的方法计算得到的绝对结合自由度与实验值之间具有非常高的相关性，因而该方法所得的结果精确度高。
在具体实施过程中，所述的预测配体受体间绝对结合自由能的高斯函数增强的自由能微扰方法具有较高的准确性，现有的自由能微扰方法仅能用于计算相对结合自由能，本发明能够计算绝对结合自由能，因而本发明所述的方法能够更广泛地应用于理性药物设计的工作中。
在具体实施过程中，拟合采样概率分布P(ΔU) data能使自由能微扰准确性及收敛性提升的方法为：
由自由能微扰的基本公式：
ΔA＝-βln∫exp(-βΔU)P 0(ΔU)dΔU；
其中，β为常数，ΔA表示两状态自由能差；采样得到的P(ΔU)需要与exp(-βΔU)相乘，结果如图5能差ΔA；在P(ΔU)与exp(-βΔU)进行相乘时，P(ΔU)处于较小ΔU的部分被与exp(-βΔU)的乘积放大，导致该部分的采样准确性对于能量计算的计算有非常大的影响；而该部分由于数据点出现的概率小，导致采样 往往不充分，因此在面临较大的微扰时会导致收敛性及准确性不足。对于采样得到的P(ΔU)进行拟合时，考虑到P(ΔU)类高斯分布的性质，用一系列高斯函数的加和进行拟合，该拟合函数充分考虑到P(ΔU)的整体性质，P(ΔU)处于较小ΔU的部分也通过外推法得到了优化，优化后的能量计算的准确性与收敛性随之极大的提升。所示，对该乘积进行积分，得到的曲线面积正比于自由
在具体实施过程中，通过本发明所述方法利用自由能微扰进行绝对自由能计算的准确度大幅提升。以针对CDK2药物靶标的16个药物小分子为例，以现有的自由能微扰方法不对概率分布P(ΔU)进行拟合计算，得到的实验值与计算值的线性回归结果如图4所示，实验与计算值之间的回归系数为0.40；而使用本发明所述方法对概率分布进行拟合计算后，得到的实验与计算值之间的回归系数为0.88，对比可知该发明专利中的自由能微扰方法比起现有的自由能微扰方法在进行药物与靶标的绝对结合自由能计算时准确度大幅提升。
在具体实施过程中，使用的高斯函数增强的自由能微扰方法进行绝对自由能计算，其结果与现有的结合自由能计算方法进行对比，对于CDK2、TYK2、JNK1、Thrombin这几个药物靶标以及相应的63个药物小分子，其对比结果如图6所示，由图可知使用高斯函数增强的自由能微扰方法进行绝对自由能计算的结果准确度与相对结合自由能计算方法相当，但绝对自由能计算不需要已知的参考药物分子，其应用范围远大于现有的相对结合自由能计算方法。
实施例2
进一步的，利用如图7、图8所示的完整版热力学循环，对本专利报道的绝对结合自由能计算方法进行测试，步骤如下：
S1：对受体蛋白-配体体系，即Rec-Lig体系进行常规的分子动力学模拟，并基于模拟获得的轨迹得到添加约束所需的距离、夹角、二面角数据；
S2：构建热力学循环，构建
Figure PCTCN2020071663-appb-000056
和
Figure PCTCN2020071663-appb-000057
根据添加约束所需的距离、夹角、二面角数据得到，并对
Figure PCTCN2020071663-appb-000058
的自由能进行解析计算；
其中：
Figure PCTCN2020071663-appb-000059
表示将溶液中药物小分子全部微扰至湮灭状态的能量；
Figure PCTCN2020071663-appb-000060
表示将药物-靶标复合物中药物小分子全部微扰至湮灭状态的能量；
Figure PCTCN2020071663-appb-000061
表示对药物-靶标复合物中药物小分子添加简谐势约束的能量变化；
S3：利用分子动力学模拟程序，对
Figure PCTCN2020071663-appb-000062
步骤 中的每一个窗口进行分子动力学模拟采样，得到
Figure PCTCN2020071663-appb-000063
Figure PCTCN2020071663-appb-000064
步骤中每一个窗口的分子动力学轨迹；
S4：根据
Figure PCTCN2020071663-appb-000065
步骤中的每一个窗口的分子动力学轨迹，通过相邻窗口的力场文件进行计算，得到每个相邻窗口之间的势能差ΔU；
S5：对于
Figure PCTCN2020071663-appb-000066
这两步，根据对应窗口的ΔU，分别计算得到
Figure PCTCN2020071663-appb-000067
步骤的自由能；
S6：对于
Figure PCTCN2020071663-appb-000068
这一步每个相邻窗口之间的势能差ΔU进行概率分布统计，利用约束概率分布函数进行拟合，得到平滑的概率分布，并根据拟合结果检查收敛性；
S7：根据拟合后的概率分布，重新生成
Figure PCTCN2020071663-appb-000069
这一步中每一窗口的势能差，通过计算得到
Figure PCTCN2020071663-appb-000070
的能量，并根据之前步骤得到的
Figure PCTCN2020071663-appb-000071
以及
Figure PCTCN2020071663-appb-000072
计算得到药-靶结合自由能ΔA binding。
在具体实施过程中，本发明提供的一种基于约束概率分布函数优化的自由能微扰方法，实现全自动得到添加约束所需的复合物拓扑文件，约束概率分布函数用于拟合实验中的概率分布，使自由能微扰(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 PCTCN2020071663-appb-000073
与1步方案差距，获得的测试结果如表1所示。
表1.选择不同数目的窗口计算得出的
Figure PCTCN2020071663-appb-000074
Figure PCTCN2020071663-appb-000075
从上表可以看出，在都以本发明提供的自动添加约束程序进行简谐势约束的前提下，设置不同的约束步骤扰动窗口数，经过约束概率分布函数拟合其概率分布后，计算得到
Figure PCTCN2020071663-appb-000076
相差较小，差值仅为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 PCTCN2020071663-appb-000077
所需窗口的分子动力学模拟参数；
S22：进一步将Rec-Lig体系中的小分子的电荷以及范德华参数依次分别减小至0，得到计算
Figure PCTCN2020071663-appb-000078
所需窗口的分子动力学模拟参数；
S23：对于Rec-Lig体系中小分子的电荷以及范德华参数依次分别减小至0，得到计算
Figure PCTCN2020071663-appb-000079
所需窗口的分子动力学模拟参数；
S24：对
Figure PCTCN2020071663-appb-000080
的自由能进行解析计算。
其中，在所述步骤S24中，所述
Figure PCTCN2020071663-appb-000081
自由能的计算公式具体为：
Figure PCTCN2020071663-appb-000082
其中，k B为玻尔兹曼常数，V为系统的体积，r aA为一个两原子间的距离用于距离简谐势约束，θ a和θ A为两个三原子间的夹角用于夹角简谐势约束，诸多k值为用于约束一个距离r aA、两个夹角θ a，θ A、三个二面角
Figure PCTCN2020071663-appb-000083
的简谐势力常数。
更具体的，所述步骤S5具体为：
根据对应窗口的ΔU，采用Benneet Acceptance Ratio计算方法，即BAR计算方法进行
Figure PCTCN2020071663-appb-000084
自由能的计算，具体计算公式为：
Figure PCTCN2020071663-appb-000085
Figure PCTCN2020071663-appb-000086
其中，ΔA i，i+1为两相邻窗口i与i+1的结合自由能差，ΔU i，i+1为基于拟合后的概率分布新生成的势能差，β＝1/k BT，k B为玻尔兹曼常数，n i与n i+1分别为两相邻状态i与i+1采样所包含的样本量，C为方程组之间的隐含变量；
将
Figure PCTCN2020071663-appb-000087
和
Figure PCTCN2020071663-appb-000088
对应的计算窗口得到的ΔA i，i+1进行加和，得到
Figure PCTCN2020071663-appb-000089
和
Figure PCTCN2020071663-appb-000090
的数值。
更具体的，所述步骤S6具体包括以下步骤：
S61：对计算
Figure PCTCN2020071663-appb-000091
的窗口中每个窗口之间的势能差进行概率分布统计，得到概率分布P(ΔU) data；
S62：基于约束概率分布函数对概率分布P(ΔU) data进行拟合，具体表达式为：
Figure PCTCN2020071663-appb-000092
其中，P(ΔU) fit为经过约束概率分布函数拟合后的平滑的概率分布，其中a i，b i，c i，d i，h i，μ i，σ i为第i个“约束概率分布函数”的参数，n为选取拟合函数的数目，n 1和n 2为两个正整数的修正因子，在本申请给出的计算案例中，取n，n 1和n 2分别为1，10和4。为说明拟合效果，图9为对添加约束过程的状态采样所得概率分布P(ΔU)被一系列约束概率分布函数拟合的结果示意图。图中灰点为采样所得概率分布点，黑线为最终拟合所得的分布曲线，灰实线和灰虚线分别为公式中的前后两项。
更具体的，所述步骤S7具体包括以下步骤：
S71：根据概率分布P(ΔU) fit，重新生成每一状态的势能差ΔU′；
S72：基于新的势能差ΔU′，利用步骤S5的计算方法重新计算得到
Figure PCTCN2020071663-appb-000093
的计算数值；
S73：根据得到的
Figure PCTCN2020071663-appb-000094
和
Figure PCTCN2020071663-appb-000095
计算出药-靶结合自由能ΔA binding。
其中，在所述步骤S73中，所述药-靶结合自由能ΔA binding的计算过程具体为：
Figure PCTCN2020071663-appb-000096
实施例3
实施例1的结果说明了高斯函数拟合可以增加
Figure PCTCN2020071663-appb-000097
和
Figure PCTCN2020071663-appb-000098
计算的收敛性，实施例2的结果说明自动添加约束的方法结合约束概率分布函数进行拟合，可以使得
Figure PCTCN2020071663-appb-000099
项能量结算结果达到非常好的收敛性。基于实施例1和实施例2的基础上，以两个针对BRD4蛋白靶标的药物小分子为例。采用自动添加约束并在约束过程中只设置一个扰动窗口的约束方案，在处理约束过程中的概率分布数据时，采用本发明中的约束概率分布函数对其概率分布散点图进行拟合。并基于图8提供的完整版热力学循环，计算绝对结合自由能ΔA binding。具体步骤如下：
S1：下载BRD4蛋白晶体结构3MXF和4HBV，准备BRD4的力场参数，选定Rec-Lig体系以及Lig体系，对Rec-Lig体系进行一小段常规的分子动力学模拟，并基于模拟获得的轨迹得到添加约束所需距离、夹角、二面角数据；
S2：根据图8所示构建热力学循环，构建
Figure PCTCN2020071663-appb-000100
步骤所需的分子动力学模拟参数，其中对
Figure PCTCN2020071663-appb-000101
这一步添加约束使用S1步骤得到的距离、夹角、二面角数据，并对
Figure PCTCN2020071663-appb-000102
这一步的自由能进行解析解计算；
S3：利用分子动力学模拟程序，对
Figure PCTCN2020071663-appb-000103
步骤中的每一个窗口进行分子动力学模拟采样，得到
Figure PCTCN2020071663-appb-000104
Figure PCTCN2020071663-appb-000105
步骤中的每一个窗口的分子动力学轨迹；
S4：根据
Figure PCTCN2020071663-appb-000106
步骤中的每一个窗口的分子动力学轨迹，通过相邻窗口的力场文件进行计算，得到每个相邻窗口之间的势能差ΔU；
S5：对于
Figure PCTCN2020071663-appb-000107
这两步，根据对应窗口的ΔU，分别计算得到
Figure PCTCN2020071663-appb-000108
步骤的自由能。
S6：对于
Figure PCTCN2020071663-appb-000109
这一步每个相邻窗口之间的势能差ΔU进行概率分布统计，利用本专利给出的“约束概率分布函数”进行拟合，得到平滑的概率分布，并根据拟合结果检查收敛性；
S7：根据拟合后的概率分布，重新生成
Figure PCTCN2020071663-appb-000110
这一步中每一窗口的势能差，通过计算得到
Figure PCTCN2020071663-appb-000111
的能量，并根据之前步骤得到的
Figure PCTCN2020071663-appb-000112
Figure PCTCN2020071663-appb-000113
以及
Figure PCTCN2020071663-appb-000114
计算得到药-靶结合自由能ΔA binding。
在具体实施过程中，根据以上步骤，考察针对BRD4蛋白靶标与两个两个药物小分子的计算结合自由能ΔA binding与实验值的差别，获得的结果如表2所示：
表2.药-靶结合自由能计算结果与实验结果对比 a
Figure PCTCN2020071663-appb-000115
a计算过程中使用自动添加约束方法
bΔA cal为计算所得药靶结合自由能
cΔA exp为实验所得药靶结合自由能
从计算与实验结果的对比可知，计算结果与实验结果非常接近，差距在2kcal/mol之内。
在最新文献中，计算
Figure PCTCN2020071663-appb-000116
至少需10个扰动窗口，且需要手动添加约束(Chem.Sci.,2016,7,207-218；J.Am.Chem.Soc.2017,139,946-957)。基于本发明的技术方案，可实现全自动得到添加简谐势约束所需的复合物拓扑文件，约 束概率分布函数可用于拟合实验所得的ΔU与其概率的散点图，使得FEP计算得到的约束窗口内的自由能变化值更精准。FEP方法是目前准确度较高的方法，但该方法需要非常大的计算量，而使用本发明的技术方案，可以在保证计算精确度的前提下，极大减少约束步骤的窗口数，使得计算窗口数目从10降至1，这同时也极大减少了计算量，计算量仅为之前的1/10。以上实施例也说明了结果的准确性。
显然，本发明的上述实施例仅仅是为清楚地说明本发明所作的举例，而并非是对本发明的实施方式的限定。对于所属领域的普通技术人员来说，在上述说明的基础上还可以做出其它不同形式的变化或变动。这里无需也无法对所有的实施方式予以穷举。凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明权利要求的保护范围之内。

Claims (9)

  1. 预测药-靶结合强度的绝对结合自由能微扰方法，其特征在于，包括以下步骤：
    S1：对受体蛋白-配体体系，即Rec-Lig体系进行常规的分子动力学模拟，并基于模拟获得的轨迹得到添加约束所需的距离、夹角、二面角数据；
    S2：构建热力学循环，构建
    Figure PCTCN2020071663-appb-100001
    和
    Figure PCTCN2020071663-appb-100002
    根据添加约束所需的距离、夹角、二面角数据得到，并对
    Figure PCTCN2020071663-appb-100003
    的自由能进行解析计算；
    其中：
    Figure PCTCN2020071663-appb-100004
    表示将溶液中药物小分子全部微扰至湮灭状态的能量；
    Figure PCTCN2020071663-appb-100005
    表示将药物-靶标复合物中药物小分子全部微扰至湮灭状态的能量；
    Figure PCTCN2020071663-appb-100006
    表示对药物-靶标复合物中药物小分子添加简谐势约束的能量变化；
    S3：利用分子动力学模拟程序，对
    Figure PCTCN2020071663-appb-100007
    步骤中的每一个窗口进行分子动力学模拟采样，得到
    Figure PCTCN2020071663-appb-100008
    Figure PCTCN2020071663-appb-100009
    步骤中每一个窗口的分子动力学轨迹；
    S4：根据
    Figure PCTCN2020071663-appb-100010
    步骤中的每一个窗口的分子动力学轨迹，通过相邻窗口的力场文件进行计算，得到每个相邻窗口之间的势能差ΔU；
    S5：对于
    Figure PCTCN2020071663-appb-100011
    这两步，根据对应窗口的ΔU的概率分布，，使用高斯函数进行拟合，并计算得到
    Figure PCTCN2020071663-appb-100012
    步骤的自由能；
    S6：对于
    Figure PCTCN2020071663-appb-100013
    这一步每个相邻窗口之间的势能差ΔU进行概率分布统计，利用约束概率分布函数进行拟合，得到平滑的概率分布，并根据拟合结果检查收敛性；
    S7：根据拟合后的概率分布，重新生成
    Figure PCTCN2020071663-appb-100014
    这一步中每一窗口的势能差，通过计算得到
    Figure PCTCN2020071663-appb-100015
    的能量，并根据之前步骤得到的
    Figure PCTCN2020071663-appb-100016
    Figure PCTCN2020071663-appb-100017
    以及
    Figure PCTCN2020071663-appb-100018
    计算得到药-靶结合自由能ΔA binding。
  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 PCTCN2020071663-appb-100019
    所需窗口的分子动力学模拟参数；
    S22：进一步将Rec-Lig体系中的小分子的电荷以及范德华参数依次分别减小至0，得到计算
    Figure PCTCN2020071663-appb-100020
    所需窗口的分子动力学模拟参数；
    S23：对于Rec-Lig体系中小分子的电荷以及范德华参数依次分别减小至0，得到计算
    Figure PCTCN2020071663-appb-100021
    所需窗口的分子动力学模拟参数；
    S24：对
    Figure PCTCN2020071663-appb-100022
    的自由能进行解析计算。
  4. 根据权利要求3所述的预测药-靶结合强度的绝对结合自由能微扰方法，其特征在于，在所述步骤S24中，所述
    Figure PCTCN2020071663-appb-100023
    自由能的计算公式具体为：
    Figure PCTCN2020071663-appb-100024
    其中，k B为玻尔兹曼常数，V为系统的体积，r aA为一个两原子间的距离用于距离简谐势约束，θ a和θ A为两个三原子间的夹角用于夹角简谐势约束，诸多k值为用于约束一个距离r aA、两个夹角θ a，θ A、三个二面角
    Figure PCTCN2020071663-appb-100025
    的简谐势力常数。
  5. 根据权利要求4所述的预测药-靶结合强度的绝对结合自由能微扰方法，其特征在于，在所述步骤S3或S4中，采用包括AMBER、Gromacs、Namd在内的任意分子动力学模拟软件进行计算。
  6. 根据权利要求4所述的预测药-靶结合强度的绝对结合自由能微扰方法，其特征在于，所述步骤S5具体包括以下步骤：
    S51：针对得到每个相邻状态之间的势能差进行概率分布统计，得到概率分布P(ΔU) data；
    S52：将概率分布P(ΔU) data利用多个高斯函数的线性组合进行拟合，得到平滑的概率分布，具体计算公式为：
    Figure PCTCN2020071663-appb-100026
    其中，P(ΔU) fit为经高斯函数拟合后的概率分布，c i为第i个高斯函数的参数，μ i为第i个高斯函数的平均值，σ i为第i个高斯函数的标准差；
    S53：根据得到的P(ΔU) fit，采用Benneet Acceptance Ratio计算方法，即BAR计算方法进行
    Figure PCTCN2020071663-appb-100027
    自由能的计算，具体计算公式为：
    Figure PCTCN2020071663-appb-100028
    Figure PCTCN2020071663-appb-100029
    其中，ΔA i，i+1为两相邻窗口i与i+1的结合自由能差，ΔU i，i+1为基于拟合后的概率分布新生成的势能差，β＝1/k BT，k B为玻尔兹曼常数，n i与n i+1分别为两相邻状态i与i+1采样所包含的样本量，C为方程组之间的隐含变量；
    S54：将
    Figure PCTCN2020071663-appb-100030
    和
    Figure PCTCN2020071663-appb-100031
    对应的计算窗口得到的ΔA i，i+1进行加和，得到
    Figure PCTCN2020071663-appb-100032
    和
    Figure PCTCN2020071663-appb-100033
    的数值。
  7. 根据权利要求6所述的预测药-靶结合强度的绝对结合自由能微扰方法，其特征在于，所述步骤S6具体包括以下步骤：
    S61：对计算
    Figure PCTCN2020071663-appb-100034
    的窗口中每个窗口之间的势能差进行概率分布统计，得到概率分布P(ΔU) data；
    S62：基于约束概率分布函数对概率分布P(ΔU) data进行拟合，具体表达式为：
    Figure PCTCN2020071663-appb-100035
    其中，P(ΔU) fit为经过约束概率分布函数拟合后的平滑的概率分布，其中a i，b i，c i，d i，h i，μ i，σ i为第i个“约束概率分布函数”的参数，n为选取拟合函数的数目，n1和n2为两个正整数的修正因子。
  8. 根据权利要求7所述的预测药-靶结合强度的绝对结合自由能微扰方法，其特征在于，所述步骤S7具体包括以下步骤：
    S71：根据概率分布P(ΔU) fit，重新生成每一状态的势能差ΔU′；
    S72：基于新的势能差ΔU′，利用步骤S5的计算方法重新计算得到
    Figure PCTCN2020071663-appb-100036
    的计算数值；
    S73：根据步骤S24、步骤S5和步骤S72得到的
    Figure PCTCN2020071663-appb-100037
    Figure PCTCN2020071663-appb-100038
    和
    Figure PCTCN2020071663-appb-100039
    计算出药-靶结合自由能ΔA binding。
  9. 根据权利要求8所述的预测药-靶结合强度的绝对结合自由能微扰方法，其特征在于，在所述步骤S73中，所述药-靶结合自由能ΔA binding的计算过程具体为：
    Figure PCTCN2020071663-appb-100040
PCT/CN2020/071663 2019-01-17 2020-01-13 预测药-靶结合强度的绝对结合自由能微扰方法 Ceased WO2020147668A1 (zh)

Applications Claiming Priority (4)

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

Publications (1)

Publication Number Publication Date
WO2020147668A1 true WO2020147668A1 (zh) 2020-07-23

Family

ID=71613710

Family Applications (1)

Application Number Title Priority Date Filing Date
PCT/CN2020/071663 Ceased WO2020147668A1 (zh) 2019-01-17 2020-01-13 预测药-靶结合强度的绝对结合自由能微扰方法

Country Status (1)

Country Link
WO (1) WO2020147668A1 (zh)

Citations (3)

* 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 中国科学院上海药物研究所 一种模拟配体分子与靶标受体反应并计算预测该反应的热力学与动力学参数的方法和系统
CN109859806A (zh) * 2019-01-17 2019-06-07 中山大学 一种预测药物-靶标结合强度的绝对自由能微扰方法

Patent Citations (3)

* 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 中国科学院上海药物研究所 一种模拟配体分子与靶标受体反应并计算预测该反应的热力学与动力学参数的方法和系统
CN109859806A (zh) * 2019-01-17 2019-06-07 中山大学 一种预测药物-靶标结合强度的绝对自由能微扰方法

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
ZHE LI , YIYOU HUANG, YINUO WU, JINGYI CHEN, DEYAN WU, CHANG-GUO ZHAN, AND HAI-BIN LUO: "Absolute Binding Free Energy Calculation and Design of a Subnanomolar Inhibitor of Phosphodiesterase-10", JOURNAL OF MEDICINAL CHEMISTRY, 28 January 2019 (2019-01-28), pages 2099 - 2111, XP055719585, ISSN: 0022-2623, DOI: 10.1021/acs.jmedchem.8b01763 *

Similar Documents

Publication Publication Date Title
CN109859806B (zh) 一种预测药物-靶标结合强度的绝对自由能微扰方法
CN111161810B (zh) 一种基于约束概率分布函数优化的自由能微扰方法
Huang et al. An accurate free energy estimator: based on MM/PBSA combined with interaction entropy for protein–ligand binding affinity
Puranen et al. Accurate conformation‐dependent molecular electrostatic potentials for high‐throughput in silico drug discovery
Card et al. A family of phosphodiesterase inhibitors discovered by cocrystallography and scaffold-based drug design
Lind et al. Support vector machines for the estimation of aqueous solubility
Fennell et al. Modeling aqueous solvation with semi-explicit assembly
Nam Acceleration of ab initio QM/MM calculations under periodic boundary conditions by multiscale and multiple time step approaches
CN110047559B (zh) 蛋白质与药物结合自由能的计算方法、系统、设备及介质
Huang et al. CHADS2, CHA2DS2-VASc and R2CHADS2 scores predict mortality in patients with coronary artery disease
Procacci et al. SAMPL6 host–guest blind predictions using a non equilibrium alchemical approach
WO2022006771A1 (zh) 分子力场多目标拟合算法库系统及工作流程方法
Wu et al. Quantitative and systems pharmacology 2. In silico polypharmacology of G protein-coupled receptor ligands via network-based approaches
Peng et al. Improving the accuracy of predicting protein–ligand binding-free energy with semiempirical quantum chemistry charge
CN110426512A (zh) 区分过氧化酶激活增殖受体γ完全激动剂、部分激动剂、拮抗剂活性的方法
Csala et al. Multiset sparse partial least squares path modeling for high dimensional omics data analysis
CN111863141A (zh) 分子力场多目标拟合算法库系统及工作流程方法
CN115881245A (zh) 虚拟筛选共价抑制剂的方法、装置、设备和计算机可读存储介质
Hu et al. Hamiltonian diversity: effectively measuring molecular diversity by shortest hamiltonian circuits
Su et al. Robust protein–ligand interaction modeling through integrating physical laws and geometric knowledge for absolute binding free energy calculation
WO2020147668A1 (zh) 预测药-靶结合强度的绝对结合自由能微扰方法
Liu et al. Divide-and-Conquer ABFE: Improving Free Energy Calculations by Enhancing Water Sampling
EP3614388A1 (en) Method and device for calculating binding free energy, and program
Chatterjee Accelerating rare events and building kinetic Monte Carlo models using temperature programmed molecular dynamics
Fukunishi et al. Protein− Ligand Binding Free Energy Calculation by the Smooth Reaction Path Generation (SRPG) Method

Legal Events

Date Code Title Description
121 Ep: the epo has been informed by wipo that ep was designated in this application

Ref document number: 20741901

Country of ref document: EP

Kind code of ref document: A1

NENP Non-entry into the national phase

Ref country code: DE

32PN Ep: public notification in the ep bulletin as address of the adressee cannot be established

Free format text: NOTING OF LOSS OF RIGHTS PURSUANT TO RULE 112(1) EPC (EPO FORM 1205N DATED 03/12/2021)

122 Ep: pct application non-entry in european phase

Ref document number: 20741901

Country of ref document: EP

Kind code of ref document: A1