CN114121148B - A method for calculating protein-ligand binding free energy based on cluster model - Google Patents

Abstract

Compared with other methods for calculating the protein-ligand binding free energy, the GFN2-xTB method is easier to use, only initial coordinates and element composition are needed to form an input structure, and the combined application of the GFN2-xTB method and the cluster model can effectively reduce the calculation time cost on the premise of ensuring relatively good accuracy. This new method based on cluster models and GFN2-xTB should have great potential in future biomacromolecule-related binding free energy calculations.

Description

A method for calculating protein-ligand binding free energy based on cluster model

technical field

The invention relates to a method for calculating protein-ligand binding free energy based on a cluster model, which belongs to the technical field of computational biology.

Background technique

Protein-ligand interactions play an important role in many biochemical processes and biomedical applications (eg, immune response, signal transduction, drug design). The study of the affinity (binding free energy) between proteins and ligands not only helps to understand the biological functions of proteins, but also has very important significance for drug development and drug mechanism research. For example, the major protease and spike protein of the current severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) virus causing the COVID-19 pandemic could serve as a good drug target. Accurate calculation of the binding free energy between SARS-CoV-2 protease or spike protein and the corresponding ligand is expected to screen out high-affinity ligands and inhibitors for effective treatment. However, protein-ligand binding is a complex process, involving thousands or even tens of thousands of atoms, and many factors will affect the calculation results of binding free energy. Therefore, how to accurately and efficiently calculate protein-ligand binding free energy is an extremely important topic in the field of computational biology.

In the past few decades, in order to accurately calculate the free energy of protein-ligand binding, various methods have been proposed and developed, mainly based on molecular mechanics, such as free energy perturbation (FEP) method and molecular mechanics/ Poisson Boltzmann Surface Area (MM/PBSA) and Molecular Mechanics/Generalized Boltzmann Surface Area (MM/GBSA) methods. Free energy perturbation (FEP) is a widely used alchemical free energy method. Over the past few years, thanks to several improvements in computing power, classical force field accuracy, sampling methods, and simulation settings, FEP can be applied accurately and reliably to the calculation of protein-ligand binding free energies. In detail, the FEP method relies on the atomic analysis of the energy difference between two similar compounds and transforms between compounds through alchemical pathways. This is a rigorous set of statistical mechanics methods that require the calculation of free energy changes in the alchemical process by slowly changing the potential energy. MM/PBSA and MM/GBSA are methods for calculating the free energy of macromolecular binding by combining molecular mechanics calculations and continuum solvation models, which are widely used to evaluate docking sites, determine structural stability, predict binding End-point free energy method for affinity. The method first performs MD simulations of protein-ligand complexes under the explicit solvent model; then removes all solvent molecules and charged ions from each MD snapshot, and uses implicit PBSA or GBSA solvent models to estimate the solvation energy; Finally, the solute conformational entropy change is calculated from a selected set of snapshots, and these individual energy components are summed to obtain the final binding free energy.

Recently, in addition to the method of molecular mechanics, the semi-empirical quantum mechanical method GFN-xTB of "geometry, frequency, non-covalent, extended tight combination" has also received extensive attention and application. This method is an extended semi-empirical tight-binding model proposed by Grimme et al. on the basis of Kohn–Sham density functional theory (DFT). Compared with other traditional semi-empirical quantum mechanics methods, the GFN2-XTB method has high computational efficiency and accuracy, and its parameterized elements cover most of the periodic table up to Z=86. The introduction of the GFN2-XTB method enables computational studies of structure optimization, vibrational frequency and non-covalent interactions for macromolecular systems with more than one thousand atoms. At present, the GFN2-xTB method has been applied to the structure optimization screening and interaction energy calculation of some chemical macromolecules, and has shown good performance, so it is expected to be extended and applied to the binding free energy calculation of protein-ligand systems .

Although there are different methods for researching and predicting the free energy of protein-ligand binding, these methods all have some corresponding deficiencies and defects. Firstly, the calculation of protein-ligand binding free energy using FEP and MM/PBSA (MM/GBSA) methods is usually cumbersome, and the corresponding force field parameters need to be found, and the simulation settings and data analysis are also relatively complicated; secondly, the FEP method to estimate the protein-ligand The accuracy of bulk binding energy is very high, the error is about 2-3kcal/mol, but it is extremely time-consuming and difficult to converge. In comparison, the calculation cost required by the MM/PBSA (MM/GBSA) method is much lower, but its accuracy is average, and the error is generally around 10kcal/mol. In a high-charge system, the error is even greater than 20kcal/mol .

In addition, the semi-empirical tight-binding GFN2-xTB method has shown certain potential in calculating macromolecular interactions, but for large protein-ligand systems, the calculation frequency of the GFN2-xTB method is extremely time-consuming due to the large amount of calculation. Not yet available directly.

Contents of the invention

In order to solve the above-mentioned technical problems, the present invention proposes a cluster model combined with GFN2-xTB method to calculate the binding free energy of protein-ligand. Compared with other methods for calculating protein-ligand binding free energy, GFN2-xTB method is easier to use and only needs initial The input structure can be formed by the composition of coordinates and elements, and the joint application with the cluster model can effectively reduce the time consumed for calculation while maintaining good calculation accuracy.

The first object of the present invention is to provide a method for calculating protein-ligand binding free energy based on the cluster model, comprising the following steps:

S1. Obtain the heavy atom information of the protein and ligand structures respectively;

S2. Open the heavy atom information of the obtained protein and ligand structure in the three-dimensional molecular model software, and intercept the protein residues whose cut-off distance is within the preset range, and save the intercepted protein-ligand structure;

S3. Perform hydrogenation saturation on the intercepted protein-ligand structure to obtain a cluster model after hydrogenation saturation;

S4. Input the cluster model after hydrogenation saturation into the xtb program, set a strong constraint on the heavy atoms of the truncated protein in the cluster model after hydrogenation saturation, and set the environmental model of protein-ligand interaction;

S5. Divide the cluster model with strong constraints into three structures: the overall cluster model, truncated protein residues, and ligands, and use the xtb program to perform Hessian calculations on the three structures, and read the overall cluster model from the output file The free energy of G _com , the free energy of truncated protein residues G _pro and the free energy of ligand G _lig , and the protein-ligand binding free energy ΔG was calculated from the above three free energies.

Further, the protein-ligand binding free energy ΔG is calculated according to the protein and ligand binding free energy formula: ΔG=G _com -G _pro -G _lig ;

Among them, ΔG is the binding free energy of protein-ligand,

G _com is the free energy of the cluster model as a whole,

G _pro is the free energy of the truncated protein residue,

G _lig is the free energy of the ligand.

范围。Further, the preset range of the cut-off distance is

scope.

Further, in the present invention, when using the xtb program to perform Hessian calculations, the GFN2-xTB method is used. For the specific calculation process, see C.Bannwarth.; S.Ehlert, et al.J.Chem.Theory Comput.2019, 15 , 1652-1671. and C. Bannwarth.; E. Caldeweyher, et al. Wiley Interdiscip. Rev.: Comput. Mol. Sci. 2021, 11, e1493.

Further, the three-dimensional molecular modeling software is one or both of PyMOL, VMD or GaussView.

Among them, VMD and PyMOL can be used to intercept proteins, and GaussView and PyMOL can hydrogenate the structure, because PyMOL can complete these two steps alone, and the generated PDB structure includes the charge information of protein residues, preferably using PyMOL software for processing.

Further, the force constant of the strong constraint is chosen to be 0.5-1.0 Hartree/Bohr ² .

Further, the environment model is an implicit water solvent model, a vacuum model, a methanol solvent model, a toluene solvent model or an ethane model.

Further, in step S1, the heavy atom information of the protein and ligand structure is obtained from the PDB-BIND website (http://www.pdbbind-cn.org/) and the protein database website (http://www.rcsb.org/ pdb/) get. For unknown protein-ligands, it is necessary to obtain possible binding sites by molecular docking method, and then calculate according to this method to obtain binding free energy.

Further, in step S3, the hydrogenation saturation is performed using three-dimensional molecular modeling software.

The beneficial effects of the present invention are:

Compared with other methods for calculating protein-ligand binding free energy, the GFN2-xTB method of the present invention is easier to use, only needs initial coordinates and element composition to form the input structure, and the joint application with the cluster model ensures a relatively good Under the premise of high accuracy, the calculation time cost is more effectively reduced. This new approach based on the cluster model and GFN2-xTB should have great potential for future calculations of binding free energies related to biomacromolecules.

Description of drawings:

Figure 1 is a calculation flow chart;

Figure 2 is the cluster model (PDB: 1A28) under different cut-off distances;

Figure 3 is the dependence of binding free energy and cut-off distance calculated by GFN2-xTB method;

Figure 4 is the MAE value and calculation time of the protein ligand binding free energy calculated using the GFN2-xTB method binding cluster model;

Figure 5 shows the corresponding energy items in the output file.

detailed description

The present invention will be further described below in conjunction with specific examples, so that those skilled in the art can better understand the present invention and implement it, but the given examples are not intended to limit the present invention.

In order to accurately and effectively use GFN2-xTB to calculate the binding free energy of protein-ligand, we adopted the cluster model. Figure 1 shows the specific calculation process. In detail, first select the protein-ligand complex to be calculated from the PDB-BIND website (http://www.pdbbind-cn.org/), and then select the protein-ligand complex from the protein database website (http://www.rcsb.org/pdb/) to download the corresponding PDB structure. Secondly, we extract the heavy atoms of the protein and ligand structure from the PDB file, and open it in the PyMOL package to intercept the cluster model, namely The protein is truncated at a different cutoff distance than the ligand; meanwhile, addition of hydrogen atoms in the PyMOL software is required to saturate the selected residues. It is worth noting that here we treat the residue as a truncated object, i.e. all atoms in the selected residue are preserved, which is beneficial to determine the total charge of the truncated protein. This will generate the input structure required by the xtb (version 6.3.3) program, the xtb package is freely available at https://github.com/grimme-lab/xtb.

After obtaining the desired input structure, we used the GFN2-xTB method of the xtb program to perform a preliminary optimization of the initial structure. In order to calculate the binding free energy of the protein-ligand interaction in the aqueous phase, we further divided the optimized structure into For a complex and two monomers, the Hessian calculations of GFN2-xTB were performed respectively. In the above two steps, we calculated using an implicit water-solvent model, that is, a generalized Born(GB) model (GBSA( _H2O )) with a surface area (SA) contribution. Furthermore, since we are using a cluster model, we placed a strong constraint on the heavy atoms of the truncated protein (the force constant was chosen to be 1.0 Hartree/Bohr ² ).

In the embodiment of the present invention, the method for calculating the binding free energy of proteins and ligands determined based on the existing coordinates in the PDB is to verify the accuracy of the method. For unknown protein-ligands, it is necessary to obtain possible binding sites by molecular docking method, and then calculate according to this method to obtain binding free energy.

Example 1: Optimal Cluster Model Scheme

时，蛋白中只有一个残基，这明显不够。随着截距的增大，残基/原子的数量越来越多，配体可以被残基包围。截断为

时，原子总数接近1500个，接近GFN2-xTB方法的上限。为确定截断值的准确可靠，我们选定了八个蛋白配体体系并使用GFN2-xTB方法计算了不同截断距离下的结合自由能，并评估计算预测的结合自由能与截断值的依赖关系。如图3所示，发现计算得到的结合自由能随着截断距离的增加而趋于平稳。当截止距离大于

时，计算出的结合自由能变化很小。在这可以判断截断距离的预设范围为

In this example, the 1A28 protein is taken as an example to show its cluster models at different cutoff distances, as shown in FIG. 2 . When the cutoff value is

, there is only one residue in the protein, which is obviously not enough. As the intercept increases, the number of residues/atoms increases and the ligand can be surrounded by residues. truncated to

, the total number of atoms is close to 1500, which is close to the upper limit of the GFN2-xTB method. In order to determine the accuracy and reliability of the cutoff value, we selected eight protein ligand systems and used the GFN2-xTB method to calculate the binding free energy at different cutoff distances, and evaluated the dependence of the predicted binding free energy on the cutoff value. As shown in Fig. 3, it was found that the calculated binding free energies leveled off as the cut-off distance increased. When the cut-off distance is greater than

, the calculated binding free energy changes very little. Here you can judge the preset range of the cut-off distance as

Example 2:

In this example, cluster models with different cutoff values combined with the GFN2-xTB method were used to calculate the binding free energies of 25 groups of protein-ligand complexes according to the above process, and further determine the cutoff range of the cluster model.

或

时，使用GFN2-xTB方法计算蛋白质-配体结合自由能的MAE为10.4kcal/mol和6.1kcal/mol，这两种方法的误差与前面所提及的MM/PBSA(MM/GBSA)误差相当，而且其还存在未将与配体相互作用的残基截取的问题，因此该截断值不适用于簇模型。而当截断值为

时，使用GFN2-xTB方法计算蛋白质-配体结合自由能的MAE值大多数都小于10.0kcal/mol，平均MAE值约为5.0kcal/mol，接近于FEP方法，具有相对较好的准确性。同时在图4b中给出了计算蛋白质-配体簇模型的所耗时长(消耗时长包括初始结构优化以及单体和复合物的频率计算，在Intel XeonPlatinum9242@2.3GHz集群上使用了12核心)，我们发现当截断值为

时，平均计算时间约为26h，时间成本较高，而截断值为

时，计算成本适中。因此，在簇模型结合GFN2-xTB方法计算蛋白质-配体结合自由能中，截断值为

都是较好的选择。此外考虑到使用截断值为

的簇模型结合GFN2-xTB方法计算时，大部分蛋白质-配体的簇模型的计算时间低于两小时，平均计算时间为4574s，可有效的降低了计算的时间成本，同时在个别情况下，能够更完全可靠的将与配体相互作用的蛋白质残基截取。因而我们认为截断值为

是一个更安全有效的选择。As shown in Figure 4a, we find that when the cutoff distance is

or

When using the GFN2-xTB method to calculate the protein-ligand binding free energy MAE is 10.4kcal/mol and 6.1kcal/mol, the errors of these two methods are comparable to the MM/PBSA (MM/GBSA) error mentioned above , and it also has the problem of not truncating residues that interact with the ligand, so this cutoff value is not suitable for cluster models. And when the cutoff value is

When using GFN2-xTB method to calculate the MAE value of protein-ligand binding free energy, most of them are less than 10.0kcal/mol, and the average MAE value is about 5.0kcal/mol, which is close to the FEP method and has relatively good accuracy. At the same time, the time taken to calculate the protein-ligand cluster model is given in Figure 4b (the time spent includes initial structure optimization and frequency calculation of monomers and complexes, using 12 cores on the Intel Xeon Platinum9242@2.3GHz cluster), We found that when the cut-off value

When , the average calculation time is about 26h, the time cost is high, and the cut-off value is

, the computational cost is moderate. Therefore, in the cluster model combined with GFN2-xTB method to calculate the protein-ligand binding free energy, the cutoff value is

are better choices. Also consider using a cutoff value of

When the cluster model combined with the GFN2-xTB method is calculated, the calculation time of most protein-ligand cluster models is less than two hours, and the average calculation time is 4574s, which can effectively reduce the calculation time cost. At the same time, in some cases, It can more completely and reliably intercept the protein residues that interact with the ligand. Therefore, we consider the cut-off value to be

is a safer and more effective option.

Example 3:

Taking the protein-ligand complex 1AU0 as an example, the binding free energy is calculated according to the method of the present invention, and compared with the experimental binding free energy calculated from the information published on the website.

(1) First select the protein-ligand complex 1AU0 from the PDB-BIND website. The basic information of the complex has been given on the page, including the protein name, corresponding ligand (SDK) and pK _d value (7.66) Wait. We can obtain the corresponding experimental binding free energy value through pK _d =-logK _d and ΔG = RTlnK _d .

(2) After obtaining the relevant information of protein 1AU0, we search and download the corresponding PDB structure from the protein database website, and need to extract the heavy atoms of the protein and ligand structure from the PDB file. If you don’t do any processing and directly use PyMOL to hydrogenate and saturate, there will be 3329 atoms. Using GFN2-xTB to directly optimize and calculate Hessian, the time cost is too high. Therefore, we used the cluster model combined with GFN2-xTB to calculate the binding free energy of protein and ligand.

范围内的蛋白质残基并保存截取好的蛋白质-配体结构。具体操作：在命令行PyMOL＞中输入“select ligand,resn SDK”和“select 6A,byres ligand expand6”。然后再保存处理好的结构：file＞Export Molecule＞selection＞6A＞PDB Options＞Save。(3) Open the PDB structure containing only protein and ligand heavy atoms in PyMOL, take the 1AU0 ligand (SDK) as the cut-off point, and intercept its

range of protein residues and preserve the truncated protein-ligand structure. Specific operation: Enter "select ligand, resn SDK" and "select 6A, byres ligand expand6" in the command line PyMOL>. Then save the processed structure: file>Export Molecule>selection>6A>PDB Options>Save.

(4) At the same time, it needs to be saturated with hydrogen and stored in PyMOL. At this time, the protein-ligand complex cluster model has only 489 atoms. This will also generate the required input structure for the xtb program, the cluster model after hydrogenation saturation.

(5) After obtaining the desired input structure, we first constrainedly optimize it using the GFN2-xTB method of the xtb program.

(6) In order to calculate the binding free energy of protein-ligand interaction in aqueous phase, we further divided the optimized structure into a complex (cluster model overall) and two monomers (truncated protein residues and ligand body), and the Hessian calculation of GFN2-xTB was performed on them respectively.

In the fifth and sixth steps, we calculate using the implicit water solvent model, namely GBSA(H ₂ O). Since we are using a cluster model, we additionally set an input file md.inp in which a strong constraint can be placed on the heavy atoms of the truncated protein (the force constant is chosen to be 1.0 Hartree/Bohr ² ). It should be noted that the Hessian calculation for the ligand does not need to read this input file. After the calculation is completed, add and subtract the corresponding energy components in the output file (as shown in Figure 5) according to the protein and ligand binding free energy formula, that is, read and calculate the corresponding energy of the complex as G _com (Figure 5a) , the corresponding energies of the two monomers (protein residue and ligand) are recorded as G _pro (Fig. 5b) and G _lig (Fig. 5c), and the corresponding ΔG is calculated as G _com -G _pro -G _lig . Thus, the final binding free energy is obtained.

The free energy formula for protein and ligand binding is:

ΔG＝G _com -G _pro -G _lig (1)

Among them, G _com is the free energy of the protein-ligand complex, G _pro is the free energy of the protein, and G _lig is the free energy of the ligand.

The corresponding calculation results of protein 1AU0 are given in Table 1. The table not only gives the binding free energy and time-consuming time of the protein-ligand cluster model calculated by the GFN2-xTB method, but also gives other energy contribution items of the free energy. It can be seen that the binding cluster model using the GFN2-xTB method to calculate the protein-ligand binding free energy can effectively reduce the time cost and maintain relatively good accuracy.

Table 1 Binding free energy results of protein 1AU0 calculated by GFN2-XTB method

PDB G_com/Eh G_pro/Eh G_lig/Eh ΔG (kcal/mol) Experiment ΔG(kcal/mol) WALL Time/s 1AU0 -776.618 -649.127 -127.480 -6.90 -10.45 4955

Note: 1Eh=627.51kcal/mol

The above-mentioned embodiments are only preferred embodiments for fully illustrating the present invention, and the protection scope of the present invention is not limited thereto. Equivalent substitutions or transformations made by those skilled in the art on the basis of the present invention are all within the protection scope of the present invention. The protection scope of the present invention shall be determined by the claims.

Claims (6)

1. A method for calculating protein-ligand binding free energy based on a cluster model, comprising the steps of:

s1, respectively obtaining heavy atom information of protein and ligand structures;

s2, opening the obtained heavy atom information of the protein and ligand structures in three-dimensional molecular model software, intercepting protein residues with the intercepting distance within a preset range, and storing the intercepted protein and the ligand together into a protein-ligand structure;

s3, carrying out hydrogenation saturation on the intercepted protein-ligand structure to obtain a cluster model after hydrogenation saturation;

s4, inputting the cluster model after hydrogenation saturation into a GFN2-xTB program, setting a strong constraint on heavy atoms of truncated proteins in the cluster model after hydrogenation saturation, and setting an environment model of protein-ligand interaction; wherein the strongly constrained force constant is selected to be 0.5-1.0 Hartree/Bohr ² (ii) a The environment model is an implicit hydrosolvent model, a vacuum model, a methanol solvent model, a toluene solvent model or an ethane model;

s5, dividing the cluster model with strong constraint into three structures of a cluster model whole body, a truncated protein residue and a ligand, respectively performing Hessian calculation on the three structures by using a GFN2-xTB method, and reading the free energy G of the cluster model whole body from an output file _com Free energy of truncated protein residue G _pro And free energy G of the ligand _lig And calculating the binding free energy delta G of the protein-ligand according to the three free energies, wherein the cluster model is a protein-ligand structure which is subjected to hydrogenation saturation and has strong constraint as a whole.

2. The method of claim 1, wherein the free energy of protein-ligand bindingΔ G according to the protein and ligand binding free energy formula: Δ G = G _com -G _pro -G _lig Calculating;

wherein AG is the binding free energy of protein-ligand,

G _com is the free energy of the cluster model as a whole,

G _pro in order to truncate the free energy of the protein residue,

G _lig is the free energy of the ligand.

3. The method of claim 1, wherein the predetermined range of the truncation distance is

A range.

4. The method according to claim 1, wherein the three-dimensional molecular modeling software is one or two of PyMOL, VMD, or GaussView.

5. The method of claim 1, wherein in step S1, the heavy atom information of the protein and ligand structures is obtained from PDB-BIND website and protein database website.

6. The method of claim 1, wherein in the step S3, the hydrogenation saturation is performed by using three-dimensional molecular modeling software.

Images