CN110047559B - Calculation method, system, device and medium for binding free energy of protein and drug - Google Patents

Abstract

The present disclosure discloses a calculation method, system, equipment and medium for the binding free energy of proteins and drugs, and constructs a function of the binding free energy of proteins and drugs, and the function of the binding free energy of proteins and drugs is the binding free energy of protein and drug complexes. The functional relationship between energy and electrostatic term, polar solvation energy term, van der Waals term, non-polar solvation energy and interaction entropy; collect the energy terms of each protein-drug complex in the training set, use each protein in the training set Each energy term of the complex with the drug and the experimental value of each protein and drug, perform multivariate linear fitting on the function of the binding free energy of the protein and the drug to obtain the function of the binding free energy of the protein and the drug; The energy terms of the protein-drug complex are input into the trained function of the binding free energy of the protein and the drug, and the binding free energy of the protein and the drug is output.

Description

Calculation method, system, device and medium for binding free energy of protein and drug

technical field

The present disclosure relates to a method, system, device and medium for calculating the binding free energy of proteins and drugs.

Background technique

The statements in this section merely mention background related to the present disclosure and do not necessarily constitute prior art.

The interactions between biomolecules, such as receptor-ligand, antigen-antibody, DNA-protein, sugar-lectin, RNA-ribosome, etc., largely determine the physiological functions of organisms, such as self-replication, metabolism and information processing, etc. Therefore, studying the recognition principle between biomolecules plays a crucial role in biology. Among them, the combination of protein and drug is the core of drug design and drug interaction. At the same time, it has direct medical research significance and application value in the experimental field and computing field in biomedicine. Especially in virtual screening, by studying the interaction between the target protein and the drug, the existing methods are used to evaluate the binding ability of the drug, and finally the drug with stronger binding ability is selected.

As far as the inventors know, in general, we use the binding free energy to describe the binding affinity between the two. At present, the methods that can accurately calculate the binding free energy are free energy perturbation (FEP) and thermodynamic integration (TI). There are many non-physical intermediate states of , which makes it difficult for the calculated results to reach a convergent state. The linear interaction energy (LIE) method is also a classical method, which uses the interaction energy and tunable parameters to estimate the binding free energy, but this method is only suitable for complexes with similar interaction mechanisms, so its The universality is relatively poor, and the computational cost is very high.

Another more common method is Molecular Mechanics/Poisson-Boltzmann Surface Area (MM/PBSA). Compared with the method, this method only needs to traverse the phase space of the two states of target and drug binding and non-binding, and does not need to consider the intermediate state, so the calculation amount is greatly reduced, so it is more commonly used. In MM/PBSA, the normalized vibration (Nmode) method is usually used to approximate the entropy change. The method employs a harmonic oscillator approximation, which treats translation, rotation, and vibration as uncoupled, and calculates the contributions of translation, rotation, and vibration. Due to the huge computational cost, the anharmonic contributions are often ignored.

In addition, studies have found that in the same trajectory, the change in vibrational entropy even reaches 5 kcal/mol [Kuhn, B.; Kollman, P.A. Binding of a Diverse Set of Ligands to Avidin and Streptavidin: An Accurate Quantitative Prediction of Their Relative Affinities by aCombination of Molecular Mechanics and Continuum Solvent Models. J. Med. Chem. 2000, 43, 3786-3791.]. Therefore, many researchers choose to ignore the contribution of entropy to the binding free energy, making the end result often less convincing.

SUMMARY OF THE INVENTION

In order to solve the deficiencies of the prior art, the present disclosure provides a method, system, device and medium for calculating the binding free energy of proteins and drugs;

In a first aspect, the present disclosure provides a method for calculating the binding free energy of proteins and drugs;

Methods for calculating the free energy of protein-drug binding, including:

Construct the function of the binding free energy of the protein and the drug, and the function of the binding free energy of the protein and the drug is the binding free energy of the protein and the drug complex and the electrostatic term, the polar solvation energy term, the van der Waals term, the non-polar solvation term The functional relationship between energy and interaction entropy;

Collect the energy terms of each protein-drug complex in the training set, including: electrostatic term, polar solvation energy term, van der Waals term, non-polar solvation energy and interaction entropy; use each protein in the training set Each energy term of the complex with the drug and the experimental value of each protein and drug, perform multivariate linear fitting on the function of the binding free energy of the protein and the drug, and obtain the function of the binding free energy of the protein and the drug;

Input the electrostatic term, polar solvation energy term, van der Waals term, non-polar solvation energy and interaction entropy of the protein-drug complex to be calculated into the fitted function of the binding free energy of protein and drug, Export the binding free energy of protein and drug.

In a second aspect, the present disclosure also provides a protein-drug binding free energy calculation system;

Protein and drug binding free energy calculation system, including:

The function building block of the binding free energy of protein and drug is configured to construct the function of binding free energy of protein and drug, and the function of binding free energy of protein and drug is the binding free energy of protein and drug complex and electrostatic term, polarity The functional relationship of solvation energy, van der Waals term, non-polar solvation energy and interaction entropy;

The multivariate linear fitting module is configured to collect the energy terms of each protein-drug complex in the training set. The energy terms include: electrostatic terms, polar solvation energy terms, van der Waals terms, non-polar solvation energy and mutual Action entropy; use the energy terms of each protein-drug complex in the training set and the experimental values of each protein and drug to perform multivariate linear fitting on the function of the binding free energy of protein and drug to obtain the binding of protein and drug. function of free energy;

The binding free energy calculation module is configured to input the electrostatic term, polar solvation energy term, van der Waals term, non-polar solvation energy and interaction entropy of the protein-drug complex to be calculated into the trained protein-drug As a function of the binding free energy of , output the binding free energy of the protein and the drug.

In a third aspect, the present disclosure also provides an electronic device, including a memory, a processor, and computer instructions stored in the memory and executed on the processor, and when the computer instructions are executed by the processor, the first aspect is completed. steps of the method.

In a fourth aspect, the present disclosure further provides a computer-readable storage medium for storing computer instructions, which, when executed by a processor, complete the steps of the method in the first aspect.

Compared with the prior art, the beneficial effects of the present disclosure are:

The present disclosure integrates the respective advantages and disadvantages of PBSA and IE, calculates various energy terms of a large number of protein and drug complexes with experimental values through PBSA and IE, and performs multivariate linear fitting between the calculated energy terms and experimental values to obtain a set of linear Fitting function, through which the binding free energy of the protein-drug complex is calculated accurately and simply.

Since the energy terms of the fitting function are calculated by the mainstream method, and are also fitted by a large number of complexes, from the final calculation results, the calculated binding free energy is numerically very close to the experimental value, The overall correlation is also much better than conventional methods.

At the same time, from the test results on the test set, the rules obtained are consistent with the rules obtained on the training set. Therefore, it has excellent accuracy and universality. Furthermore, since the main advantage of the IE method (ie, fast computation speed) is fully utilized in the present disclosure, the present disclosure also has the advantage of low computational cost.

Based on the above-mentioned advantages, the present disclosure will shine in describing the affinity of protein-drug-like systems. Especially in the virtual screening of drugs, it will provide a computationally reliable, widely applicable and low-cost new avenue for this.

Description of drawings

The accompanying drawings that form a part of the present application are used to provide further understanding of the present application, and the schematic embodiments and descriptions of the present application are used to explain the present application and do not constitute improper limitations on the present application.

FIG. 1 is a method flow diagram of one or more embodiments.

Detailed ways

It should be noted that the following detailed description is exemplary and intended to provide further explanation of the application. Unless otherwise defined, all technical and scientific terms used in this disclosure have the same meaning as commonly understood by one of ordinary skill in the art to which this application belongs.

It should be noted that the terminology used herein is for the purpose of describing specific embodiments only, and is not intended to limit the exemplary embodiments according to the present application. As used herein, unless the context clearly dictates otherwise, the singular is intended to include the plural as well, furthermore, it is to be understood that when the terms "comprising" and/or "including" are used in this specification, it indicates that There are features, steps, operations, devices, components and/or combinations thereof.

As shown in FIG. 1 , the first embodiment, an embodiment of the present disclosure provides a method for calculating the binding free energy of proteins and drugs;

Methods for calculating the free energy of protein-drug binding, including:

Construct the function of the binding free energy of the protein and the drug, and the function of the binding free energy of the protein and the drug is the binding free energy of the protein and the drug complex and the electrostatic term, the polar solvation energy term, the van der Waals term, the non-polar solvation term The functional relationship between energy and interaction entropy;

Collect the energy terms of each protein-drug complex in the training set, including: electrostatic term, polar solvation energy term, van der Waals term, non-polar solvation energy and interaction entropy; use each protein in the training set Each energy term of the complex with the drug and the experimental value of each protein and drug, perform multivariate linear fitting on the function of the binding free energy of the protein and the drug, and obtain the function of the binding free energy of the protein and the drug;

Input the electrostatic term, polar solvation energy term, van der Waals term, non-polar solvation energy and interaction entropy of the protein-drug complex to be calculated into the fitted function of the binding free energy of protein and drug, Export the binding free energy of protein and drug.

Further, the method also includes:

Input the electrostatic term, polar solvation energy term, van der Waals term, non-polar solvation energy and interaction entropy of the protein-drug complex of the test set into the fitted function of the binding free energy of protein and drug, The estimated value of binding free energy of protein and drug is output; by comparing the relationship between the calculated estimated value and the experimental value, the rationality, generality and stability of the fitting function are evaluated.

Further, the specific steps of constructing the function of the binding free energy of the protein and the drug are:

PBSA_IE=a(ΔE _ele +ΔG _pb )+bΔE _vdw +cΔG _np +dIE+e

Among them, PBSA_IE represents the binding free energy of the protein and the drug; ΔG _pb represents the polar solvation energy term; ΔE _vdw represents the van der Waals energy term of the protein and the drug; ΔG _np represents the non-polar solvation energy term; IE represents the protein-drug complex The interaction entropy of ; a, b, c, d, and e represent the parameters to be fitted, respectively.

Further, before performing multivariate linear fitting, the steps of preprocessing the training set and the test set include:

The same preprocessing steps were used for both the training set and the test set. The protein library IDs (PDB IDs) of the proteins and drugs in the training set and the test set were all derived from PDB bind (URL: http://www.pdbbind-cn.org/ ) in the "Protein-ligand complexes: The general set minus refined set" dataset. From this dataset, we randomly selected 84 complexes from a total of 1047 complexes from 1a to 1u in alphabetical order as a training set, and randomly selected 44 complexes from a total of 860 complexes from 1v to 2h as a test set. The training and test sets do not contain complexes with more than 2 charges or containing metal ions). The crystal structure of the complex was obtained from the PDB (website: https://www.rcsb.org/) database, and the experimental values were queried from PDB bind.

(110) Drug pretreatment:

Two-step quantum processing of drugs using Gaussian 03 software:

The first step: optimize the drug through the HF/6-31G** method of the Gaussian 03 software to find the optimal structure, that is, the structure with the lowest energy;

Step 2: Perform single-point energy calculation on the Gaussian group B3LYP/cc-PVTZ to obtain the electrostatic potential of the drug; fit the charge of each atom by the constrained electrostatic potential method;

(120) Protein pretreatment:

Since the obtained protein crystal structure generally does not contain the information of hydrogen atoms, the missing hydrogen atoms of the protein are supplemented by the LEaP module in AMBER 12;

(130) Pretreatment of protein and drug complexes:

同时，根据每个体系的带电性，加入反向离子对体系进行电性平衡；The protein and drug complexes are placed in the AMBER ff12SB force field and placed in a periodic TIP3P water box. The minimum distance between the solute and the edge of the water box is

At the same time, according to the electrification of each system, counter ions are added to balance the electrical properties of the system;

Further, each system in the training set and test set uses the AMBER package for MD simulation. The specific simulation methods are as follows:

(140) Energy optimization of protein and drug complex system:

First use the steepest descent method to minimize the energy of the protein-drug complex, and then use the conjugate gradient minimization method until the energy of the system reaches convergence;

(150) composite system is warming up to normal temperature:

将每个体系的温度从0K逐步升高到300K，利用朗之万动力学调节温度，碰撞频率为1.0ps^-1，所有涉及氢原子的键都通过SHAKE算法进行约束，模拟步长为2fs；A 300ps restricted MD simulation was performed with the constraining force constant of

The temperature of each system was gradually increased from 0K to 300K, the temperature was adjusted by Langevin dynamics, the collision frequency was 1.0ps ^-1 , all bonds involving hydrogen atoms were constrained by the SHAKE algorithm, and the simulation step size was 2fs;

(160) Non-restrictive MD simulation of composite system:

Perform an unrestricted MD simulation of the protein-drug system with a duration of 2ns and steps of 2fs:

For the first 1 ns, the trajectory is recorded every 2000 steps, and the trajectory contains the information of all atomic positions in the water molecules of the complex, counter ion and water box. The purpose of this MD simulation is to move the complex to a relatively balanced conformation;

The trajectory was recorded every 5 steps for the next 1 ns. The trajectory only contained the position information of the complex, and the energy terms of the combined free energy were obtained by sampling the trajectory.

Further, the specific steps for collecting the energy terms of each protein and drug system in the training set and the test set are as follows:

(170)

为复合物的静电能，

为蛋白质的静电能，

为药物的静电能；其通过如下公式得出：in,

is the electrostatic energy of the complex,

is the electrostatic energy of the protein,

is the electrostatic energy of the drug; it is given by the following formula:

Among them, N is the total number of atoms, q _i and q _j are the charge amounts of the i-th and j-th atoms, and R _ij is the distance between the i-th atom and the j-th atom.

(180) The expression of the van der Waals energy term for proteins and drugs is:

和

分别为复合物、蛋白质和药物的范德华势，其通过如下公式得出in,

and

are the van der Waals potentials of the complex, protein, and drug, respectively, which are given by

Among them, A _ij and B _ij are the Lanna-Jones potentials between the ith atom and the j th atom, and A _ij and B _ij are derived from the AMBERff12SB force field.

(190) Polar solvation energy term:

ΔG _pb is the polar solvation energy term, obtained by solving the Poisson-Boltzmann (PB) equation from the implicit solvent model:

where ε(r) is the dielectric constant inside or outside the molecule when r falls inside or outside the molecular surface; k(r) is the Debye-Huckel constant in the solvent region; ρ(r) is the charge of the solute molecule Distribution function; Φ(r) is the potential distribution at r. In this method, this term is obtained by solving the PB equation with the PBSA program in the AMBER software package, and the internal and external dielectric constants are set to 1 and 80, respectively.

(200) Non-polar solvation energy term:

It is derived from the solvent accessible surface area SASA:

ΔG _np = γSASA+β

和0.92kcal/mol。where the empirical constants γ and β are respectively set as

and 0.92kcal/mol.

(210) The interaction entropy of protein-drug complexes:

为蛋白质与药物的相互作用能，其可以通过对MD模拟进行平均来评估：where K is the Boltzmann constant, T is the temperature,

is the protein-drug interaction energy, which can be estimated by averaging the MD simulations:

but

是蛋白质与药物相互作用能在平均能量附近波动的振幅；Among them, N is the total number of conformation frames obtained, i is the conformation of the i-th frame, β is a constant, t _i is the time of the i-th frame, and <.> represents the mathematical symbol for averaging;

is the amplitude of the protein-drug interaction energy fluctuation around the mean energy;

Based on the concept of interaction entropy (IE), the present disclosure uses it to calculate the entropy change of the complex instead of the conventional Nmode method. It is derived from a strict theoretical formula, and it has a more convincing theoretical basis than the Normal mode method derived approximately.

Second, it directly calculates the average value of the interaction entropy through the kinetically simulated conformational trajectory of the protein-drug complex, so the interaction energy of the complex in the "gas phase" of each frame of simulation can be used for the calculation of IE, so in The stability of computational entropy is also better than that of current computational entropy methods.

Furthermore, the Normal mode method, which is based on the absolute entropy value of complexes, proteins and drugs, often causes numerical errors due to too large absolute values, while the IE method perfectly avoids this shortcoming by directly calculating the entropy change of the system. It also greatly speeds up the calculation. Therefore, the method has the characteristics of strong reliability and low computational cost.

Due to the inherent defects of PBSA, its calculation results are often inaccurate and reliable. Although there is an excellent correlation between the theoretical and experimental values calculated by combining the PBSA and IE methods (ie, PBSA+IE method), there is a large absolute error in the numerical value, and the final binding free energy is often overestimated. .

The present disclosure belongs to the scoring function method based on PBSA and interaction entropy IE. By performing multivariate linear fitting on a training set consisting of a large number of randomly selected protein and drug complexes, each energy term is given different weights, thereby making up for the PBSA is inherently flawed to accurately predict the binding free energies of other protein-drug complexes. At the same time, in order to verify the universality of the method, we randomly selected multiple protein-drug complexes with no repetition from the training set as the test set to evaluate the scoring function.

Combining the energy terms of free energy to calculate

The binding free energy is expressed as the sum of the enthalpy term and entropy:

ΔG=ΔH-TΔS (1)

Since the present disclosure adopts the IE method to calculate the entropy term, the above formula can be written as:

ΔG=ΔH+IE (2)

PBSA calculates the enthalpy term

The free energy of protein and drug can be regarded as two parts: gas-phase binding free energy and solvation energy:

ΔH=ΔE _MM +ΔG _sol (3)

where ΔE _MM is the gas phase energy, and its form is as follows:

ΔE _MM = ΔE _ele + ΔE _vdw (4)

where ΔE _ele , its expression is:

和

分别为复合物、蛋白质和药物的静电能，其通过如下公式得出in,

and

are the electrostatic energies of the complex, protein and drug, respectively, which are obtained by the following formula

Among them, N is the total number of atoms, q _i and q _j are the charges of the ith and jth atoms, and R _ij is the distance between the ith and jth atoms (the same below). In formula (4), ΔE _vdw is the van der Waals energy term of protein and drug:

和

分别为复合物、蛋白质和药物的范德华势，其通过如下公式得出：in,

and

are the van der Waals potentials of the complex, protein, and drug, respectively, which are given by:

Among them, A _ij and B _ij are the Lanna-Jones potentials between the ith atom and the jth atom, which are derived from the AMBER ff12SB force field.

In formula (3), ΔG _sol is the solvation energy of the complex, which can be divided into the following two parts:

ΔG _sol = ΔG _pb + ΔG _np (9)

where ΔG _pb is the polar solvation energy term, obtained by solving the Poisson-Boltzmann (PB) equation from the implicit solvent model:

where ε(r) is the dielectric constant inside or outside the molecule when r falls inside or outside the molecular surface; k(r) is the Debye-Huckel constant in the solvent region; ρ is the charge distribution of the solute molecule; Φ( r) is the potential distribution at r. In this method, this term is obtained by solving the PB equation with the PBSA program in the AMBER software package, and the internal and external dielectric constants are set to 1 and 80, respectively.

ΔG _np in Eq. (9) is the non-polar solvation energy term, which is derived from the empirical solvent accessible surface area SASA:

ΔG _np = γSASA+β (11)

和0.92kcal/mol。所计算的构象是从最后1ns中的构象空间中每隔1000帧取出共100帧构象用于计算PBSA。In the embodiment of the present application, the empirical constants γ and β are respectively set as

and 0.92kcal/mol. The calculated conformations are taken from the conformation space in the last 1 ns every 1000 frames for a total of 100 frames for calculating PBSA.

IE calculates the entropy term

In the method of interaction entropy adopted in the implementation of this application, the gas phase part of the binding free energy is derived by the following formula:

和

分别为蛋白质与药物、蛋白质-水和药物-水的相互作用能。

是蛋白质与药物相互作用能的平均值，

是蛋白质与药物相互作用能在平均能量附近波动的振幅。因此，我们将IE定义为：where E _p , _El and E _w are the internal energies of protein, drug and water, respectively,

and

are the interaction energies of protein and drug, protein-water and drug-water, respectively.

is the average value of protein-drug interaction energy,

is the amplitude of the fluctuation of the protein-drug interaction energy around the mean energy. Therefore, we define IE as:

为蛋白质与药物的相互作用能，其可以通过对MD模拟进行平均来评估：where K is the Boltzmann constant, T is the temperature,

is the protein-drug interaction energy, which can be estimated by averaging the MD simulations:

but

N is the total number of conformation frames obtained, i is the conformation of the ith frame, β is a constant, and t _i is the time of the ith frame.

Multivariate Linear Fitting: PBSA_IE

According to the essential relationship between the energy terms of the binding free energy, the present disclosure adopts the following method to perform linear fitting:

PBSA_IE=a(ΔE _ele +ΔG _pb )+bΔE _vdw +cΔG _np +dIE+e (16)

In the second embodiment, an embodiment of the present disclosure also provides a protein-drug binding free energy calculation system;

Protein and drug binding free energy calculation system, including:

The function building block of the binding free energy of protein and drug is configured to construct the function of binding free energy of protein and drug, and the function of binding free energy of protein and drug is the binding free energy of protein and drug complex and electrostatic term, polarity The functional relationship of solvation energy, van der Waals term, non-polar solvation energy and interaction entropy;

The multivariate linear fitting module is configured to collect the energy terms of each protein-drug complex in the training set. The energy terms include: electrostatic terms, polar solvation energy terms, van der Waals terms, non-polar solvation energy and mutual Action entropy; use the energy terms of each protein-drug complex in the training set and the experimental values of each protein and drug to perform multivariate linear fitting on the function of the binding free energy of protein and drug to obtain the binding of protein and drug. function of free energy;

The binding free energy calculation module is configured to input the electrostatic term, polar solvation energy term, van der Waals term, non-polar solvation energy and interaction entropy of the protein-drug complex to be calculated into the trained protein-drug As a function of the binding free energy of , output the binding free energy of the protein and the drug.

In the third embodiment, the embodiment of the present disclosure also provides an electronic device, including a memory, a processor, and computer instructions stored in the memory and executed on the processor, and when the computer instructions are executed by the processor, the first step is completed. Steps of the method of an embodiment.

The present disclosure also provides an electronic device, including a memory, a processor, and computer instructions stored in the memory and executed on the processor, and when the computer instructions are executed by the processor, each operation in the method is completed. For brevity, It is not repeated here.

It should be understood that in the present disclosure, the processor may be a central processing unit CPU, and the processor may also be other general-purpose processors, digital signal processors DSP, application-specific integrated circuits ASIC, off-the-shelf programmable gate array FPGA or other Programming logic devices, discrete gate or transistor logic devices, discrete hardware components, etc. A general purpose processor may be a microprocessor or the processor may be any conventional processor or the like.

The memory may include read-only memory and random access memory and provide instructions and data to the processor, and a portion of the memory may also include non-volatile random access memory. For example, the memory may also store device type information.

In the fourth embodiment, an embodiment of the present disclosure further provides a computer-readable storage medium for storing computer instructions, and when the computer instructions are executed by a processor, the steps of the method described in the first embodiment are completed.

In the implementation process, each step of the above-mentioned method can be completed by a hardware integrated logic circuit in a processor or an instruction in the form of software. The steps of the method disclosed in conjunction with the present disclosure can be directly embodied as executed by a hardware processor, or executed by a combination of hardware and software modules in the processor. The software modules may be located in random access memory, flash memory, read-only memory, programmable read-only memory or electrically erasable programmable memory, registers and other storage media mature in the art. The storage medium is located in the memory, and the processor reads the information in the memory, and completes the steps of the above method in combination with its hardware. To avoid repetition, detailed description is omitted here. Those of ordinary skill in the art can realize that the units, ie algorithm steps, of each example described in conjunction with the embodiments disclosed in the present disclosure can be implemented by electronic hardware or a combination of computer software and electronic hardware. Whether these functions are performed in hardware or software depends on the specific application and design constraints of the technical solution. Skilled artisans may implement the described functionality using different methods for each particular application, but such implementations should not be considered beyond the scope of this application.

Those skilled in the art can clearly understand that, for the convenience and brevity of description, the specific working process of the above-described systems, devices and units may refer to the corresponding processes in the foregoing method embodiments, which will not be repeated here.

In the several embodiments provided in this application, it should be understood that the disclosed system, apparatus and method may be implemented in other manners. For example, the apparatus embodiments described above are only illustrative. For example, the division of the units is only a division of a logical function. In actual implementation, there may be other division methods, for example, multiple units or components may be combined Either it can be integrated into another system, or some features can be omitted, or not implemented. In addition, the shown or discussed mutual coupling or direct coupling or communication connection may be through some interfaces, indirect coupling or communication connection of devices or units, which may be in electrical, mechanical or other forms.

The functions, if implemented in the form of software functional units and sold or used as separate products, may be stored in a computer-readable storage medium. Based on this understanding, the technical solution of the present application can be embodied in the form of a software product in essence, or the part that contributes to the prior art or the part of the technical solution. The computer software product is stored in a storage medium, including Several instructions are used to cause a computer device (which may be a personal computer, a server, or a network device, etc.) to execute all or part of the steps of the methods described in the various embodiments of the present application. The aforementioned storage medium includes: U disk, mobile hard disk, Read-Only Memory (ROM, Read-Only Memory), Random Access Memory (RAM, Random Access Memory), magnetic disk or optical disk and other media that can store program codes .

In conclusion, the PBSA_IE method can extensively, stably and accurately calculate the binding free energy of protein-drug complexes. At the same time, the advantages and disadvantages of PBSA and the method are also integrated, and they also have good advantages in computational cost.

The above descriptions are only preferred embodiments of the present application, and are not intended to limit the present application. For those skilled in the art, the present application may have various modifications and changes. Any modification, equivalent replacement, improvement, etc. made within the spirit and principle of this application shall be included within the protection scope of this application.

Claims (6)

1. A method for calculating free energy of protein and drug combination, which is characterized by comprising the following steps:

constructing a function of the free energy of protein and drug binding as a function of the free energy of protein and drug complex as a function of electrostatic, polar solvation energy, van der waals, nonpolar solvation energy, and entropy of interaction;

acquiring energy items of each protein and drug compound in the training set, wherein the energy items comprise: electrostatic terms, polar solvation energy terms, van der waals terms, nonpolar solvation energy, and entropy of interaction; performing multi-element linear fitting on the function of the free energy of the combination of the protein and the medicine by using each energy item of each protein and medicine compound in a training set and the experimental value of each protein and medicine to obtain the function of the free energy of the combination of the protein and the medicine;

inputting the electrostatic term, the polar solvation energy term, the van der waals term, the nonpolar solvation energy and the interaction entropy of the protein and drug compound to be calculated into the fitted function of the free energy of the protein and the drug, and outputting the free energy of the protein and the drug;

the specific steps for constructing the function of the free energy of the combination of the protein and the drug are as follows:

PBSA_IE＝a(ΔE_ele+ΔG_pb)+bΔE_vdw+cΔG_np+dIE+e

wherein PBSA _ IE represents the free energy of the protein combined with the drug; Δ G_pbRepresents a polar solvation energy term; delta E_vdwRepresents the van der waals energy term of protein and drug; Δ G_npRepresents a non-polar solvation energy term; IE represents the entropy of interaction of the protein with the drug complex; a. b, c, d and e respectively represent parameters to be fitted;

each system in the training set adopts an AMBER program package to perform MD simulation, and the specific simulation mode is as follows:

(140) optimizing the energy of a protein and drug composite system: minimizing the energy of the protein-drug compound by a steepest descent method, and then minimizing by a conjugate gradient until the energy of the system reaches convergence;

(150) heating the composite system to normal temperature: a restrictive MD simulation was performed at 300ps with a constraint constant of

The temperature of each system was gradually raised from 0K to 300K, the temperature was adjusted using Langew dynamics, and the collision frequency was 1.0ps^-1All bonds related to hydrogen atoms are constrained by a SHAKE algorithm, and the simulation step length is 2 fs;

(160) composite system non-limiting MD simulation: non-limiting MD simulations with a 2ns step size of 2fs were performed on protein-drug systems: the first 1ns recorded a trace every 2000 steps, which contained information on the position of all atoms in the water molecules of the complex, the counterion and the water box, the purpose of this MD simulation being to move the complex to a relatively equilibrium conformation; recording the track every 5 steps in the last 1ns, wherein the track only contains the position information of the compound, and each energy item combining free energy is obtained by sampling the track;

the polar solvation energy term is obtained by solving a Poisson-Boltzmann (PB) equation through an implicit solvent model:

wherein ε (r) is the dielectric constant of the molecule when r falls inside or outside the surface of the molecule; k (r) represents the Debye-Huckel constant in the solvent region; ρ (r) represents the charge distribution function of the solute molecules; Φ (r) is the potential distribution at r;

the non-polar solvation energy term: it is given by the solvent accessible surface area SASA:

ΔG_np＝γSASA+β

wherein the empirical constants γ and β are set to

And 0.92 kcal/mol;

entropy of interaction of the protein with the drug complex:

wherein K is Boltzmann constant, T is temperature,

as the interaction energy of the protein with the drug, it can be evaluated by averaging MD simulations:

then

Where N is the total number of acquired constellations, i is the ith frame constellation, β is a constant, t_iIs the time of the i-th frame,<.>a mathematical symbol representing an averaging;

is the amplitude of the fluctuation of the protein interaction energy with the drug around the mean energy.

2. The method of claim 1, wherein the step of preprocessing the training set prior to performing the multivariate linear fit comprises:

the same pretreatment steps are adopted in the training set and the testing set, and PDB IDs of the protein and the medicine in the training set and the testing set are derived from a PDB bind database;

(110) pretreatment of the medicine:

the medicine is subjected to two-step quantum processing by utilizing Gaussian 03 software:

the first step is as follows: optimizing the medicine by using an HF/6-31G method of Gaussian 03 software to find an optimal structure, namely a structure with the lowest energy;

the second step is that: performing single-point energy calculation on a Gaussian group B3LYP/cc-PVTZ to obtain the electrostatic potential of the medicine; fitting the charge of each atom by a constrained electrostatic potential method;

(120) protein pretreatment:

supplementing the missing hydrogen atoms of the protein through a LEAP module in AMBER 12;

(130) pretreatment of a composite system:

the protein and drug complex was placed in an AMBERff 12SB force field and placed in a periodic TIP3P water box with solute at a minimum distance from the edges of the box

Meanwhile, according to the charging property of each system, reverse ions are added to carry out electric balance on the system.

3. The method of claim 1, wherein the step of collecting the energy terms for each protein and drug system in the training set and test set comprises:

(170)

wherein,

is the electrostatic energy of the composite and is,

is the electrostatic energy of the protein and is,

is the electrostatic energy of the drug; it is given by the following formula:

wherein N is the total number of atoms, q_i、q_jIs the charge of the ith and jth atom, R_ijIs the distance between the ith atom and the jth atom;

(180) the expression of the van der waals energy term for proteins and drugs is:

wherein,

and

van der Waals potentials of the complex, protein and drug, respectively, are given by the following formula

Wherein A is_ijIs the Lanna-Jones potential between the ith atom and the jth atom, B_ijIs a Lanna-Jones potential between the ith atom and the jth atom, A_ij、B_ijDerived from the AMBER ff12SB force field.

4. A protein and drug binding free energy calculation system, comprising:

a function construction module of the free energy of combination of the protein and the drug, which is configured to construct a function of the free energy of combination of the protein and the drug, wherein the function of the free energy of combination of the protein and the drug is a function relation of the free energy of combination of the protein and the drug and electrostatic terms, polar solvation energy, van der waals terms, nonpolar solvation energy and interaction entropy;

a multivariate linear fit module configured to acquire energy terms for each protein and drug complex in the training set, the energy terms including: electrostatic terms, polar solvation energy terms, van der waals terms, nonpolar solvation energy, and entropy of interaction; performing multi-element linear fitting on the function of the free energy of the combination of the protein and the medicine by using each energy item of each protein and medicine compound in a training set and the experimental value of each protein and medicine to obtain the function of the free energy of the combination of the protein and the medicine;

a binding free energy calculation module configured to input an electrostatic term, a polar solvation energy term, a van der waals term, a non-polar solvation energy and an interaction entropy of the protein-drug complex to be tested into a function of trained binding free energies of the protein and the drug, and output the binding free energies of the protein and the drug;

the specific steps for constructing the function of the free energy of the combination of the protein and the drug are as follows:

PBSA_IE＝a(ΔE_ele+ΔG_pb)+bΔE_vdw+cΔG_np+dIE+e

wherein PBSA _ IE represents the free energy of the protein combined with the drug; Δ G_pbRepresents a polar solvation energy term; delta E_vdwRepresents the van der waals energy term of protein and drug; Δ G_npRepresents a non-polar solvation energy term; IE represents the entropy of interaction of the protein with the drug complex; a. b, c, d and e respectively represent parameters to be fitted;

each system in the training set adopts an AMBER program package to perform MD simulation, and the specific simulation mode is as follows:

(140) optimizing the energy of a protein and drug composite system: minimizing the energy of the protein-drug compound by a steepest descent method, and then minimizing by a conjugate gradient until the energy of the system reaches convergence;

(150) heating the composite system to normal temperature: a restrictive MD simulation was performed at 300ps with a constraint constant of

The temperature of each system was gradually raised from 0K to 300K, the temperature was adjusted using Langew dynamics, and the collision frequency was 1.0ps^-1All bonds related to hydrogen atoms are constrained by a SHAKE algorithm, and the simulation step length is 2 fs;

(160) composite system non-limiting MD simulation: non-limiting MD simulations with a 2ns step size of 2fs were performed on protein-drug systems: the first 1ns recorded a trace every 2000 steps, which contained information on the position of all atoms in the water molecules of the complex, the counterion and the water box, the purpose of this MD simulation being to move the complex to a relatively equilibrium conformation; recording the track every 5 steps in the last 1ns, wherein the track only contains the position information of the compound, and each energy item combining free energy is obtained by sampling the track;

the polar solvation energy term is obtained by solving a Poisson-Boltzmann (PB) equation through an implicit solvent model:

wherein ε (r) is the dielectric constant of the molecule when r falls inside or outside the surface of the molecule; k (r) represents the Debye-Huckel constant in the solvent region; ρ (r) represents the charge distribution function of the solute molecules; Φ (r) is the potential distribution at r;

the non-polar solvation energy term: it is given by the solvent accessible surface area SASA:

ΔG_np＝γSASA+β

wherein the empirical constants γ and β are set to

And 0.92 kcal/mol;

entropy of interaction of the protein with the drug complex:

wherein K is Boltzmann constant, T is temperature,

as the interaction energy of the protein with the drug, it can be evaluated by averaging MD simulations:

then

Where N is the total number of acquired constellations, i is the ith frame constellation, β is a constant, t_iIs the time of the i-th frame,<.>a mathematical symbol representing an averaging;

is the amplitude of the fluctuation of the protein interaction energy with the drug around the mean energy.

5. An electronic device comprising a memory and a processor and computer instructions stored on the memory and executable on the processor, the computer instructions when executed by the processor performing the steps of the method of any of claims 1 to 3.

6. A computer-readable storage medium storing computer instructions which, when executed by a processor, perform the steps of the method of any one of claims 1 to 3.

Images