CN109859806A - A kind of prediction drug-target bond strength absolute freedom energy perturbation method - Google Patents
A kind of prediction drug-target bond strength absolute freedom energy perturbation method Download PDFInfo
- Publication number
- CN109859806A CN109859806A CN201910045003.XA CN201910045003A CN109859806A CN 109859806 A CN109859806 A CN 109859806A CN 201910045003 A CN201910045003 A CN 201910045003A CN 109859806 A CN109859806 A CN 109859806A
- Authority
- CN
- China
- Prior art keywords
- lig
- state
- drug
- probability distribution
- rec
- 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.)
- Granted
Links
- 238000000034 method Methods 0.000 title claims abstract description 49
- 239000003596 drug target Substances 0.000 title claims abstract description 22
- 238000009826 distribution Methods 0.000 claims abstract description 48
- 238000005381 potential energy Methods 0.000 claims abstract description 45
- 239000003814 drug Substances 0.000 claims abstract description 37
- 238000004364 calculation method Methods 0.000 claims abstract description 34
- 229940079593 drug Drugs 0.000 claims abstract description 33
- 238000000329 molecular dynamics simulation Methods 0.000 claims abstract description 17
- 238000005070 sampling Methods 0.000 claims abstract description 11
- 239000003446 ligand Substances 0.000 claims description 30
- 150000003384 small molecules Chemical class 0.000 claims description 29
- 238000012986 modification Methods 0.000 claims description 8
- 230000004048 modification Effects 0.000 claims description 8
- 230000001172 regenerating effect Effects 0.000 claims description 6
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Substances O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 claims description 6
- 238000012900 molecular simulation Methods 0.000 claims description 4
- 238000002791 soaking Methods 0.000 claims description 4
- 238000005411 Van der Waals force Methods 0.000 claims description 3
- 239000011449 brick Substances 0.000 claims description 3
- 238000004088 simulation Methods 0.000 abstract description 4
- 230000003993 interaction Effects 0.000 description 9
- 230000008569 process Effects 0.000 description 7
- 108010024986 Cyclin-Dependent Kinase 2 Proteins 0.000 description 5
- 102100036239 Cyclin-dependent kinase 2 Human genes 0.000 description 5
- 238000010586 diagram Methods 0.000 description 5
- 238000009510 drug design Methods 0.000 description 5
- 230000000694 effects Effects 0.000 description 3
- 239000002547 new drug Substances 0.000 description 3
- 150000001875 compounds Chemical class 0.000 description 2
- 230000006872 improvement Effects 0.000 description 2
- 238000000324 molecular mechanic Methods 0.000 description 2
- 101000950695 Homo sapiens Mitogen-activated protein kinase 8 Proteins 0.000 description 1
- 101000844245 Homo sapiens Non-receptor tyrosine-protein kinase TYK2 Proteins 0.000 description 1
- 102100037808 Mitogen-activated protein kinase 8 Human genes 0.000 description 1
- 102100032028 Non-receptor tyrosine-protein kinase TYK2 Human genes 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 238000005094 computer simulation Methods 0.000 description 1
- 230000003247 decreasing effect Effects 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 238000013461 design Methods 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 238000009509 drug development Methods 0.000 description 1
- 230000000857 drug effect Effects 0.000 description 1
- 230000009881 electrostatic interaction Effects 0.000 description 1
- UVCJGUGAGLDPAA-UHFFFAOYSA-N ensulizole Chemical compound N1C2=CC(S(=O)(=O)O)=CC=C2N=C1C1=CC=CC=C1 UVCJGUGAGLDPAA-UHFFFAOYSA-N 0.000 description 1
- 238000002474 experimental method Methods 0.000 description 1
- 238000013213 extrapolation Methods 0.000 description 1
- 238000011478 gradient descent method Methods 0.000 description 1
- 150000002500 ions Chemical class 0.000 description 1
- 238000012417 linear regression Methods 0.000 description 1
- 229920002521 macromolecule Polymers 0.000 description 1
- 238000003032 molecular docking Methods 0.000 description 1
- 230000007935 neutral effect Effects 0.000 description 1
- 230000000737 periodic effect Effects 0.000 description 1
- 229920009537 polybutylene succinate adipate Polymers 0.000 description 1
- 238000012827 research and development Methods 0.000 description 1
- 230000007704 transition Effects 0.000 description 1
- 238000003041 virtual screening Methods 0.000 description 1
Abstract
The present invention provides a kind of prediction drug-target bond strength absolute freedom energy perturbation method, comprising: building Lig system and Rec-Lig system；Molecular dynamics simulation is carried out to the parameter of each state of Lig system and Rec-Lig system, obtains the molecular dynamics track of each state；Carrying out single-point by the field of force file of adjacent states can calculate, and obtain the potential energy difference between each adjacent states；Probability distribution statistical is carried out to the potential energy difference between each adjacent states, is fitted using the linear combination of multiple Gaussian functions；According to the probability distribution after fitting, the potential energy difference of each state is regenerated, the Conjugated free energy of drug and target is obtained by calculation.Absolute freedom energy perturbation method provided by the invention, probability distribution corresponding to the dynamics simulation track of all free energy perturbation intermediate state is fitted using Gaussian function, so that the probability distribution of sampling is more reasonable, the convergence and accuracy of energy computation results are greatly improved.
Description
Technical Field
The invention relates to the technical field of drug research and development, in particular to an absolute free energy perturbation method for predicting drug-target binding strength.
Background
The reasons for the drug effect are almost closely related to the mutual combination of the drug and the target macromolecules. Therefore, free energy prediction technology of drug-target binding strength plays an extremely important role in drug design (Q.Rev.Biophys.2012, 45, 1-25.; curr.Opin.struct.biol.2011, 21, 150-. If a method is available to accurately predict binding free energy, the cycle and cost of new drug development will be greatly reduced.
Although there are some fast combined free energy prediction methods, such as docking scoring function, MM/PBSA (molecular mechanics/poisson boltzmann surface, an energy calculation method commonly used in drug design), MM/GBSA (molecular mechanics/generalized berne surface), etc. However, because these methods introduce many empirical parameters or the theoretical basis has certain defects, their accuracy is not high. Free Energy Perturbation (FEP), as a free energy calculation method based on statistical mechanics, theoretically has higher accuracy. For example, the famous medicine design software company Schrodinger invented a relative binding free energy prediction method (j.am.chem.soc.2015, 137, 2695-doping 2703.) in 2015 based on free energy perturbation, but the method can only predict the relative binding free energy of a homoframework compound (the difference is less than 10 atoms) and a target which only has a tiny structure modification on the existing medicine molecules, and the prediction method needs to refer to a small molecular structure with known activity and cannot predict the absolute binding free energy of different framework compounds and the target. The prediction of absolute binding free energy is required for most of the problem cases encountered in the development of new drugs, such as virtual screening, framework transition and full-new drug design.
Disclosure of Invention
The invention provides an absolute free energy perturbation method for predicting the binding strength of a drug-target, aiming at overcoming the technical problem that the absolute binding free energy prediction cannot be carried out in the prior art.
In order to solve the technical problems, the technical scheme of the invention is as follows:
a method of predicting absolute free energy perturbation of drug-target binding strength comprising the steps of:
s1: constructing a drug small molecular ligand Lig system and a drug-target Rec-Lig system, and preparing for intermediate perturbation state parameters;
s2: performing molecular dynamics simulation on parameters of each state of the Lig system and the Rec-Lig system by using a molecular dynamics simulation program to obtain a molecular dynamics track of each state of the Lig system and the Rec-Lig system;
s3: carrying out single-point energy calculation on the molecular dynamics track of each state of the obtained Lig system and Rec-Lig system through the force field files of adjacent states to obtain the potential energy difference delta U between the adjacent states;
s4: carrying out probability distribution statistics on potential energy difference between each two adjacent states, and fitting by utilizing linear combination of a plurality of Gaussian functions to obtain smooth probability distribution;
s5: and regenerating the potential energy difference of each state according to the fitted probability distribution, and calculating to obtain the binding free energy of the drug and the target.
Wherein, the step S1 specifically includes the following steps:
s11: preparing a structure of a receptor-ligand complex;
s12: calculating the charge of the small molecule ligand and assigning a force field to the ligand;
s13: assigning a force field to a receptor in the receptor-ligand complex, soaking the complex in a water tank, preparing an initial system for molecular simulation of the complex, recording the initial system as a Rec-Lig system, and exporting a force field file as an initial state parameter of the Rec-Lig system;
s14: soaking an individual small molecule ligand in a water tank, preparing an individual small molecule molecular simulation initial system, recording the system as a Lig system, and exporting a force field file as an initial state parameter of the Lig system;
s15: modifying a force field corresponding to the Rec-Lig system, reducing the charges of all atoms in a micromolecular ligand in a force field file to 0 through multiple modifications, and storing a new force field file in each step as a middle perturbation state parameter of the Rec-Lig system; the same modification is carried out on the force field file of the Lig system, the charges of all atoms of the small molecular ligand are reduced to 0, and the parameters of the middle perturbation state of the Lig system are stored;
s16: after the small molecule charge is 0, reducing Van der Waals parameters of all atoms in the small molecule ligand in a Rec-Lig system force field file to 0, and storing a new force field file as an intermediate perturbation state parameter of the Rec-Lig system in each step; and (3) carrying out the same treatment on the force field file of the Lig system, reducing the Van der Waals parameters of all atoms in the micromolecular ligand in the force field file to 0, and saving a new force field file as an intermediate perturbation state parameter of the Lig system in each step.
Wherein, in step S2, the molecular dynamics simulation program is one of AMBER program, Gromacs program and Namd program.
Wherein, the step S3 includes the following steps:
s31: respectively calculating the single point energy of each frame by using the force fields of the (i + 1) th state and the (i-1) th state based on the molecular dynamics track of the (i) th state;
s32: and calculating potential energy difference delta U of the i-th to i + 1-th states and the i-th to i-1-th states based on the track of each frame according to the single point energy of each frame.
Wherein, the step S4 specifically includes:
s41: carrying out probability distribution statistics on potential energy difference between each two adjacent states to obtain probability distribution P (delta U)data；
S42: will probability distribution P (Δ U)dataFitting by utilizing the linear combination of a plurality of Gaussian functions to obtain smooth probability distribution, wherein the specific calculation formula is as follows:
wherein, P (delta U)fitIs probability distribution after Gaussian function fitting, ciIs a parameter of the ith Gaussian function, muiIs the mean value, σ, of the ith Gaussian functioniIs the standard deviation of the ith gaussian function.
Wherein, the step S5 specifically includes the following steps:
s51: regenerating potential energy difference delta U' of each state according to the probability distribution obtained by fitting;
s52: based on the new potential energy difference delta U', solving an equation set to obtain the difference delta A of the combined free energy between the adjacent statesi，i+1The specific equation is as follows:
wherein ,ΔAi，i+1Is the combined free energy difference, delta U, of two adjacent states i and i +1i，i+1For the newly generated potential energy difference based on the fitted probability distribution, β ═ 1/kBT，kBIs the Boltzmann constant, niAnd ni+1Respectively sampling sample quantities contained in two adjacent states i and i +1, wherein C is an implicit variable between equation sets;
s53: adding potential energy differences of all states of the Rec-Lig system to obtain energyAnd adding energy differences of all states of the Lig system to obtain energyCalculating to obtain the binding free energy delta A of the drug and the targetbindingThe specific calculation formula is as follows:
wherein ,represents the energy required to perturb all small drug molecules in solution to an annihilation state;represents the energy that perturbs all the small drug molecules in the drug-target complex to annihilate the bricks; the annihilation state is that the charge and van der waals force L-J potential energy of the drug small molecule are 0, and the surrounding yards cannot sense the existence of the drug small molecule.
Compared with the prior art, the technical scheme of the invention has the beneficial effects that:
according to the absolute free energy perturbation method for predicting the drug-target bonding strength, probability distribution corresponding to the dynamic simulation track of all free energy perturbation intermediate states is fitted by utilizing a Gaussian function, so that the probability distribution of sampling is more reasonable, and the convergence and the accuracy of an energy calculation result are greatly improved. The improvement of convergence and accuracy enables the free energy perturbation method to process larger perturbation, annihilates the whole drug micromolecules, and then calculates the absolute free energy between the drug and the target, so as to achieve the purpose of greatly improving the application range in the free energy perturbation re-drug design.
Drawings
FIG. 1 is a schematic flow diagram of the process;
FIG. 2 is a diagram of the thermodynamic cycle involved in calculating absolute binding free energy;
FIG. 3 is a schematic diagram of a sampled converted probability distribution P (Δ U) of Δ U;
FIG. 4 is a graph showing the results of a series of Gaussian function fits to the probability distribution P (Δ U) obtained from the sampling;
FIG. 5 shows exp (- β Δ U) P during energy calculation0(Delta U) and P0(Δ U) a graph showing the relationship;
FIG. 6 is a schematic diagram of the difference in method accuracy before and after fitting a Gaussian function to a sampling probability distribution;
fig. 7 is a schematic diagram of the accuracy comparison between the absolute bound free energy calculation method enhanced by the gaussian function and the relative bound free energy calculation method of schrodinger formula in the present invention.
Detailed Description
The drawings are for illustrative purposes only and are not to be construed as limiting the patent;
for the purpose of better illustrating the embodiments, certain features of the drawings may be omitted, enlarged or reduced, and do not represent the size of an actual product;
it will be understood by those skilled in the art that certain well-known structures in the drawings and descriptions thereof may be omitted.
The technical solution of the present invention is further described below with reference to the accompanying drawings and examples.
Example 1
As shown in fig. 1, an absolute free energy perturbation method for predicting drug-target binding strength comprises the following steps:
s1: constructing a drug small molecular ligand Lig system and a drug-target Rec-Lig system, and preparing for intermediate perturbation state parameters;
s2: performing molecular dynamics simulation on parameters of each state of the Lig system and the Rec-Lig system by using a molecular dynamics simulation program to obtain a molecular dynamics track of each state of the Lig system and the Rec-Lig system;
s3: carrying out single-point energy calculation on the molecular dynamics track of each state of the obtained Lig system and Rec-Lig system through the force field files of adjacent states to obtain the potential energy difference delta U between the adjacent states;
s4: carrying out probability distribution statistics on potential energy difference between each two adjacent states, and fitting by utilizing linear combination of a plurality of Gaussian functions to obtain smooth probability distribution;
s5: and regenerating the potential energy difference of each state according to the fitted probability distribution, and calculating to obtain the binding free energy of the drug and the target.
Example 2
Based on example 1, further, the basic principle of free energy calculation using free energy perturbation is shown in fig. 2, wherein Rec stands for Receptor, i.e. drug-target Receptor; lig stands for Ligand, a pharmaceutical small molecule Ligand; dummy stands for Dummy, i.e. imaginary atoms, all of which do not have any interaction with the surrounding environment. The binding free energy delta A of the Lig system and the Rec-Lig system is obtained by calculationbindingRespectively obtaining the thermodynamic cycles as shown in FIG. 2Andafter the size of (2), the binding free energy Delta A can be obtainedbinding。
Further, for example, 16 small drug molecules against CDK2 drug target have the following specific steps:
s1: constructing a drug small molecular ligand Lig system and a drug-target Rec-Lig system, and preparing for intermediate perturbation state parameters, specifically:
s11: preparing a structure of a receptor-ligand complex, and butting 16 drug small-molecule ligands into a binding pocket of a macromolecular CDK2 receptor by using a molecular butting method;
s12: calculating ESP charges of the small molecular ligand at HF/6-31G × level by Gaussian software, fitting RESP charges by using Antechamber, and assigning GAFF force field parameters;
s13: addition of receptor-assigned AMBER03 force field in receptor-ligand complexes to the entire complexAdding counter ions into the water tank of the chamfered octahedron to keep the whole system electrically neutral, and storing the coordinates and initial force field parameters of the finally generated system as a force field file as initial state parameters of free energy perturbation of the Rec-Lig system;
s14: adding separate small molecule ligandsThe water tank of the chamfer octahedron saves the coordinate of the system finally generated and the initial force field parameter as a force field file, and the force field file is used as the parameter of the initial state of free energy perturbation of the Lig system;
s15: annihilate the charges of small molecules in the Rec-Lig system and the Lig system. For both systems, the charge of all atoms in the small molecule is reduced to 0 by 5 steps, respectively. Namely, the charge of each atom in the micromolecule is modified to be 0.8, 0.6, 0.4, 0.2 and 0 in each step, and the parameter file after the charge is modified in each step is respectively saved as an independent force field file to be used as the parameter of the intermediate perturbation state;
s16: after annihilation of the small molecule charge, the small molecule also retains van der Waals interactions with the outside world. It is therefore necessary to annihilate the Rec-Lig system and the van der waals interactions of small molecules in Lig systems. The intermolecular van der Waals interactions are calculated as the Lennard-Jones (LJ) potential energy by the formula:
wherein ,εijIs original toLJ potential energy parameter between child i and atom j, rijIs the distance between atom i and atom j; in this embodiment, by decreasing εijTo attenuate the van der waals interactions between atoms; since sampling of these non-true states becomes more difficult with annihilation of van der waals interactions, as van der waals interactions decrease to a smaller scale, closer proximity between the different states is required to improve the convergence of the sampling. Thus epsilon for each atom is in the range of (1-0.2i)6Is reduced. After 5 steps of modification of van der Waals interaction are carried out in this way, all the effects of the small molecules on the outside are completely eliminated, and the effect of annihilation of the small molecules is achieved. Storing all parameter files with modified epsilon as independent force field files as parameters of the intermediate perturbation state;
s2: and (4) performing molecular dynamics simulation on the Lig system and the Rec-Lig system by using an original force field parameter file and a force field file in an intermediate perturbation state of the past life of the step S1 respectively by using a molecular dynamics simulation program AMBER, a Gromacs program or a Namd program. The system cannot be limited by using the shift function in the simulation, because the intermediate state is unstable due to the characteristics of the shift function, and therefore the step size of the simulation is selected to be 1 fs. In addition, the periodic boundary condition of the system needs to be calculated by adopting a PME method for the long-range electrostatic interaction, and a cutoff value of 8A is selected for the short-range van der Waals interaction; molecular dynamics simulations were performed for 4ns at 300K for each state of each system. In the simulation process of the last 2ns, extracting the conformation of one frame every 100fs, wherein the track of each state comprises 20000 frame conformations;
s3: and (3) carrying out single-point energy calculation on the molecular dynamics track of each state of the obtained Lig system and Rec-Lig system through the force field files of adjacent states to obtain the potential energy difference delta U between the adjacent states, which specifically comprises the following steps:
solving the system average of the potential energy difference between the state i +1 and the state i on the basis of the molecular dynamics track of the ith state, thereby solving forward perturbation potential energy difference; solving the system average of the potential energy difference between the state i-1 and the state i, wherein the solved potential energy difference is a backward perturbation potential energy difference; the results of the forward perturbation potential energy difference and the backward perturbation potential energy difference are combined, a BAR method is used for calculation, and the accuracy of final energy calculation is effectively improved.
S4: carrying out probability distribution statistics on potential energy difference between each adjacent state, fitting by utilizing linear combination of a plurality of Gaussian functions to obtain smooth probability distribution, specifically:
after step S3, the trajectory of each state contains 20000 constellations, where you sum the probability distributions of all potential energy differences Δ U. Therefore, 20000 data points are obtained after the track of each state is calculated; as shown in fig. 3, specifically:
s41: firstly, the potential energy difference delta U is converted into corresponding probability distribution P (delta U)data；
S42: using a plurality of linear combinations of Gaussian functions to match probability distribution P (Δ U)dataFitting is carried out, and a specific calculation formula is as follows:
wherein, P (Delta U)fitIs probability distribution after Gaussian function fitting, ciIs a parameter of the ith Gaussian function, muiIs the mean value, σ, of the ith Gaussian functioniIs the standard deviation of the ith Gaussian function; contains 15 parameters c in totali，μi，σiWhere i is 1, 2, 3, 4, 5, thereby defining a cost function of the form:
wherein, P (Delta U)dataProbability distribution for kinetic sampling, P (Δ U)fitIs the probability distribution after fitting by a gaussian function,ciis a parameter of the ith Gaussian function, muiIs the mean value, σ, of the ith Gaussian functioniIs the standard deviation of the ith Gaussian function; optimizing the cost function by using a random gradient descent method, specifically optimizing the 15 parameters c by using an SGD methodi，μi，σiWherein i is 1, 2, 3, 4, 5, specifically:
wherein ,(ci，μi，σi) t is a parameter value after t iterations; obtained by continuously losing the band until the t iteration (c)i，μi，σi)tI.e. the most significant parameters of the fitted function. After fitting, the probability distribution of the results is shown in fig. 4.
S5: according to the fitted probability distribution, the potential energy difference of each state is regenerated, and the binding free energy of the drug and the target is obtained through calculation, specifically:
s51: regenerating 500000 to 1000000 new potential energy differences delta U' by using probability distribution after Gaussian function fitting;
s52: solving an equation set based on the newly generated potential energy difference delta U' to obtain the difference delta A of the combined free energy between the adjacent statesi，i+1The specific equation is as follows:
wherein ,ΔAi，i+1Is the combined free energy difference, delta U, of two adjacent states i and i +1i，i+1For the newly generated potential energy difference based on the fitted probability distribution, β ═1/kBT，kBIs the Boltzmann constant, niAnd ni+1Respectively sampling sample quantities contained in two adjacent states i and i +1, wherein C is an implicit variable between equation sets;
s53: adding potential energy differences of all states of the Rec-Lig system to obtain energyAnd adding energy differences of all states of the Lig system to obtain energyCalculating to obtain the binding free energy delta A of the drug and the targetbindingThe specific calculation formula is as follows:
wherein ,represents the energy required to perturb all small drug molecules in solution to an annihilation state;represents the energy that perturbs all the small drug molecules in the drug-target complex to annihilate the bricks; the annihilation state is that the charge and van der waals force L-J potential energy of the drug small molecule are 0, and the surrounding yards cannot sense the existence of the drug small molecule.
In the specific implementation process, according to the steps S1-S5, 16 small molecules of CDK2 system are calculated, the structures of the 16 small molecules are known, and IC50The formula for the experimental value of the binding free energy, which can be determined experimentally, is:
ΔG＝-RTlnIC50；
the obtained result is shown in fig. 4, and by contrast, the absolute binding freedom calculated by the method of the invention has very high correlation with the experimental value, so that the method has high accuracy of the obtained result.
In the specific implementation process, the free energy perturbation method for predicting the enhancement of the Gaussian function of the absolute binding free energy between ligand receptors has higher accuracy, the existing free energy perturbation method can only be used for calculating the relative binding free energy, and the method can calculate the absolute binding free energy, so that the method can be more widely applied to the work of rational drug design.
In a specific implementation, a sample probability distribution P (Δ U) is fitteddataThe method for improving the accuracy and the convergence of the free energy perturbation comprises the following steps:
the basic formula perturbed by free energy:
ΔA＝-βln∫exp(-βΔU)P0(ΔU)dΔU；
β is a constant, delta A represents the free energy difference of two states, sampled P (delta U) needs to be multiplied by exp (- β delta U), the result is shown in FIG. 5, the product is integrated, the area of the obtained curve is proportional to the free energy difference delta A, when P (delta U) is multiplied by exp (- β delta U), the part of P (delta U) in smaller delta U is amplified by the product of exp (- β delta U), the sampling accuracy of the part has great influence on the calculation of energy calculation, the part is often insufficient due to the small probability of data point occurrence, the convergence and the accuracy are insufficient when large perturbation is faced, when the sampled P (delta U) is simulated, the property of Gaussian distribution like P (delta U) is considered, the fitting is carried out by the addition and the addition of a series of Gaussian functions, the overall property of P (delta U) is fully considered, the property of P (delta U) in smaller delta U is also considered, the calculation accuracy of the energy calculation is greatly improved by the extrapolation method.
In the specific implementation process, the accuracy of absolute free energy calculation by using free energy perturbation is greatly improved. Taking 16 small drug molecules aiming at the CDK2 drug target as an example, fitting calculation is performed on the probability distribution P (Δ U) by using the existing free energy perturbation method, and the obtained linear regression result of the experimental value and the calculated value is shown in fig. 6, and the regression coefficient between the experimental value and the calculated value is 0.40; after the method is used for fitting and calculating probability distribution, the regression coefficient between the obtained experiment and the calculated value is 0.88, and the comparison shows that the accuracy of the free energy perturbation method in the invention is greatly improved when the absolute combination free energy of the medicine and the target is calculated compared with the existing free energy perturbation method, and the accuracy of the existing method is not enough to effectively guide medicine design.
In the specific implementation process, the absolute free energy calculation is performed by using a gaussian function enhanced free energy perturbation method, the result is compared with the existing bound free energy calculation method, the comparison result is shown in fig. 7 for several drug targets, namely CDK2, TYK2, JNK1 and thrombobin, and corresponding 63 drug small molecules, and it can be known that the accuracy of the result of the absolute free energy calculation performed by using the gaussian function enhanced free energy perturbation method is equivalent to that of the relative bound free energy calculation method, but the absolute free energy calculation does not need known reference drug molecules, and the application range is far larger than that of the existing relative bound free energy calculation method.
It should be understood that the above-described embodiments of the present invention are merely examples for clearly illustrating the present invention, and are not intended to limit the embodiments of the present invention. Other variations and modifications will be apparent to persons skilled in the art in light of the above description. And are neither required nor exhaustive of all embodiments. Any modification, equivalent replacement, and improvement made within the spirit and principle of the present invention should be included in the protection scope of the claims of the present invention.
Claims (6)
1. A method of predicting absolute free energy perturbation of drug-target binding strength comprising the steps of:
s1: constructing a drug small molecular ligand Lig system and a drug-target Rec-Lig system, and preparing for intermediate perturbation state parameters;
s2: performing molecular dynamics simulation on parameters of each state of the Lig system and the Rec-Lig system by using a molecular dynamics simulation program to obtain a molecular dynamics track of each state of the Lig system and the Rec-Lig system;
s3: carrying out single-point energy calculation on the molecular dynamics track of each state of the obtained Lig system and Rec-Lig system through the force field files of adjacent states to obtain the potential energy difference delta U between the adjacent states;
s4: carrying out probability distribution statistics on potential energy difference between each two adjacent states, and fitting by utilizing linear combination of a plurality of Gaussian functions to obtain smooth probability distribution;
s5: and regenerating the potential energy difference of each state according to the fitted probability distribution, and calculating to obtain the binding free energy of the drug and the target.
2. The method according to claim 1, wherein the step S1 specifically comprises the following steps:
s11: preparing a structure of a receptor-ligand complex;
s12: calculating the charge of the small molecule ligand and assigning a force field to the ligand;
s13: assigning a force field to a receptor in the receptor-ligand complex, soaking the complex in a water tank, preparing an initial system for molecular simulation of the complex, recording the initial system as a Rec-Lig system, and exporting a force field file as an initial state parameter of the Rec-Lig system;
s14: soaking an individual small molecule ligand in a water tank, preparing an individual small molecule molecular simulation initial system, recording the system as a Lig system, and exporting a force field file as an initial state parameter of the Lig system;
s15: modifying a force field corresponding to the Rec-Lig system, reducing the charges of all atoms in a micromolecular ligand in a force field file to 0 through multiple modifications, and storing a new force field file in each step as a middle perturbation state parameter of the Rec-Lig system; the same modification is carried out on the force field file of the Lig system, the charges of all atoms of the small molecular ligand are reduced to 0, and the parameters of the middle perturbation state of the Lig system are stored;
s16: after the small molecule charge is 0, reducing Van der Waals parameters of all atoms in the small molecule ligand in a Rec-Lig system force field file to 0, and storing a new force field file as an intermediate perturbation state parameter of the Rec-Lig system in each step; and (3) carrying out the same treatment on the force field file of the Lig system, reducing the Van der Waals parameters of all atoms in the micromolecular ligand in the force field file to 0, and saving a new force field file as an intermediate perturbation state parameter of the Lig system in each step.
3. The absolute free energy perturbation method for predicting drug-target binding strength of claim 2 wherein in step S2, the molecular dynamics simulation program is one of AMBER program, gromas program, and Namd program.
4. The method according to claim 3, wherein the step S3 is specifically performed by the following steps:
s31: respectively calculating the single point energy of each frame by using the force fields of the (i + 1) th state and the (i-1) th state based on the molecular dynamics track of the (i) th state;
s32: and calculating potential energy difference delta U of the i-th to i + 1-th states and the i-th to i-1-th states based on the track of each frame according to the single point energy of each frame.
5. The method according to claim 4, wherein the step S4 is specifically as follows:
s41: carrying out probability distribution statistics on potential energy difference between each two adjacent states to obtain probability distribution P (delta U)data；
S42: will probability distribution P (Δ U)dataFitting by utilizing the linear combination of a plurality of Gaussian functions to obtain smooth probability distribution, wherein the specific calculation formula is as follows:
wherein, P (Delta U)fitIs the probability distribution after fitting by a gaussian function,ciis a parameter of the ith Gaussian function, muiIs the mean value, σ, of the ith Gaussian functioniIs the standard deviation of the ith gaussian function.
6. The method according to claim 5, wherein the step S5 specifically comprises the following steps:
s51: regenerating potential energy difference delta U' of each state according to the probability distribution obtained by fitting;
s52: based on the new potential energy difference delta U', solving an equation set to obtain the difference delta A of the combined free energy between the adjacent statesi,i+1The specific equation is as follows:
wherein ,ΔAi,i+1Is the combined free energy difference, delta U, of two adjacent states i and i +1i,i+1For the newly generated potential energy difference based on the fitted probability distribution, β ═ 1/kBT，kBIs the Boltzmann constant, niAnd ni+1Respectively sampling sample quantities contained in two adjacent states i and i +1, wherein C is an implicit variable between equation sets;
s53: adding potential energy differences of all states of the Rec-Lig system to obtain energyAnd adding energy differences of all states of the Lig system to obtain energyCalculating to obtain the binding free energy delta A of the drug and the targetbindingThe specific calculation formula is as follows:
wherein ,represents the energy required to perturb all small drug molecules in solution to an annihilation state;represents the energy that perturbs all the small drug molecules in the drug-target complex to annihilate the bricks; the annihilation state is that the charge and van der waals force L-J potential energy of the drug small molecule are 0, and the surrounding yards cannot sense the existence of the drug small molecule.
Priority Applications (2)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201910045003.XA CN109859806B (en) | 2019-01-17 | 2019-01-17 | Absolute free energy perturbation method for predicting drug-target binding strength |
| PCT/CN2020/071663 WO2020147668A1 (en) | 2019-01-17 | 2020-01-13 | Absolute binding free energy perturbation method for predicting drug target binding strength |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201910045003.XA CN109859806B (en) | 2019-01-17 | 2019-01-17 | Absolute free energy perturbation method for predicting drug-target binding strength |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN109859806A true CN109859806A (en) | 2019-06-07 |
| CN109859806B CN109859806B (en) | 2023-09-22 |
Family
ID=66894914
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201910045003.XA Active CN109859806B (en) | 2019-01-17 | 2019-01-17 | Absolute free energy perturbation method for predicting drug-target binding strength |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN109859806B (en) |
Cited By (9)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN110767262A (en) * | 2019-09-18 | 2020-02-07 | 复旦大学 | Structure-based optimized design method for aptamer |
| CN111161810A (en) * | 2019-12-31 | 2020-05-15 | 中山大学 | Free energy perturbation method based on constraint probability distribution function optimization |
| CN111341391A (en) * | 2020-02-25 | 2020-06-26 | 深圳晶泰科技有限公司 | Free energy perturbation computing and scheduling method used in heterogeneous cluster environment |
| WO2020147668A1 (en) * | 2019-01-17 | 2020-07-23 | 中山大学 | Absolute binding free energy perturbation method for predicting drug target binding strength |
| CN112199909A (en) * | 2020-10-22 | 2021-01-08 | 深圳晶泰科技有限公司 | Method for accurately calculating absolute free energy of gas molecules |
| WO2021031545A1 (en) * | 2020-02-25 | 2021-02-25 | 深圳晶泰科技有限公司 | Free energy perturbation calculation scheduling method used in heterogeneous cluster environment |
| CN114187971A (en) * | 2021-12-10 | 2022-03-15 | 上海智药科技有限公司 | Molecular free energy calculation and stability analysis method, device, equipment and storage medium |
| CN114550844A (en) * | 2022-01-28 | 2022-05-27 | 厦门大学 | Method for accelerating acidity constant and oxidation-reduction potential based on machine learning potential energy |
| CN115116537A (en) * | 2022-08-29 | 2022-09-27 | 香港中文大学（深圳） | Method and system for calculating multiple transformation paths of biomolecule functional dynamics |
Citations (5)
| 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 |
| WO2005029351A1 (en) * | 2003-09-22 | 2005-03-31 | Algodign, Llc | Method of modeling and predicting binding of ligand molecules to target molecules by methods of quantum mechanics taking effect of solvent into account |
| CN102222177A (en) * | 2011-07-08 | 2011-10-19 | 上海生物信息技术研究中心 | Auxiliary prediction method for molecular modification of antibody protein |
| US20140249787A1 (en) * | 2011-10-06 | 2014-09-04 | Universitat De Barcelona | Method of exploring the flexibility of macromolecular targets and its use in rational drug design |
| US20150095007A1 (en) * | 2012-03-20 | 2015-04-02 | University Of Maryland, Baltimore | Site-specific fragment identification guided by single-step free energy perturbation calculations |
-
2019
- 2019-01-17 CN CN201910045003.XA patent/CN109859806B/en active Active
Patent Citations (5)
| 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 |
| WO2005029351A1 (en) * | 2003-09-22 | 2005-03-31 | Algodign, Llc | Method of modeling and predicting binding of ligand molecules to target molecules by methods of quantum mechanics taking effect of solvent into account |
| CN102222177A (en) * | 2011-07-08 | 2011-10-19 | 上海生物信息技术研究中心 | Auxiliary prediction method for molecular modification of antibody protein |
| US20140249787A1 (en) * | 2011-10-06 | 2014-09-04 | Universitat De Barcelona | Method of exploring the flexibility of macromolecular targets and its use in rational drug design |
| US20150095007A1 (en) * | 2012-03-20 | 2015-04-02 | University Of Maryland, Baltimore | Site-specific fragment identification guided by single-step free energy perturbation calculations |
Non-Patent Citations (1)
| Title |
|---|
| 曹冉等: "计算化学方法在基于受体结构的药物分子设计中的基础理论及应用", 《药学学报》 * |
Cited By (13)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| WO2020147668A1 (en) * | 2019-01-17 | 2020-07-23 | 中山大学 | Absolute binding free energy perturbation method for predicting drug target binding strength |
| CN110767262A (en) * | 2019-09-18 | 2020-02-07 | 复旦大学 | Structure-based optimized design method for aptamer |
| CN111161810A (en) * | 2019-12-31 | 2020-05-15 | 中山大学 | Free energy perturbation method based on constraint probability distribution function optimization |
| CN111161810B (en) * | 2019-12-31 | 2022-03-22 | 中山大学 | Free energy perturbation method based on constraint probability distribution function optimization |
| CN111341391B (en) * | 2020-02-25 | 2023-12-01 | 深圳晶泰科技有限公司 | Free energy perturbation calculation scheduling method for heterogeneous cluster environment |
| CN111341391A (en) * | 2020-02-25 | 2020-06-26 | 深圳晶泰科技有限公司 | Free energy perturbation computing and scheduling method used in heterogeneous cluster environment |
| WO2021031545A1 (en) * | 2020-02-25 | 2021-02-25 | 深圳晶泰科技有限公司 | Free energy perturbation calculation scheduling method used in heterogeneous cluster environment |
| CN112199909A (en) * | 2020-10-22 | 2021-01-08 | 深圳晶泰科技有限公司 | Method for accurately calculating absolute free energy of gas molecules |
| CN114187971A (en) * | 2021-12-10 | 2022-03-15 | 上海智药科技有限公司 | Molecular free energy calculation and stability analysis method, device, equipment and storage medium |
| CN114550844A (en) * | 2022-01-28 | 2022-05-27 | 厦门大学 | Method for accelerating acidity constant and oxidation-reduction potential based on machine learning potential energy |
| CN114550844B (en) * | 2022-01-28 | 2024-04-30 | 厦门大学 | Method for accelerating acidity constant and oxidation-reduction potential based on machine learning potential energy |
| CN115116537A (en) * | 2022-08-29 | 2022-09-27 | 香港中文大学（深圳） | Method and system for calculating multiple transformation paths of biomolecule functional dynamics |
| CN115116537B (en) * | 2022-08-29 | 2022-12-06 | 香港中文大学（深圳） | Method and system for calculating multiple transformation paths of biomolecule functional dynamics |
Also Published As
| Publication number | Publication date |
|---|---|
| CN109859806B (en) | 2023-09-22 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN109859806A (en) | A kind of prediction drug-target bond strength absolute freedom energy perturbation method | |
| Homeyer et al. | Binding free energy calculations for lead optimization: assessment of their accuracy in an industrial drug design context | |
| Duarte et al. | Resolving apparent conflicts between theoretical and experimental models of phosphate monoester hydrolysis | |
| KR102071491B1 (en) | Breast cancer prognosis prediction method and system based on machine learning using next generation sequencing | |
| Merz Jr | Using quantum mechanical approaches to study biological systems | |
| Zheng et al. | Adaptive quantum mechanics/molecular mechanics methods | |
| Nam et al. | Specific reaction parametrization of the AM1/d Hamiltonian for phosphoryl transfer reactions: H, O, and P atoms | |
| Scheff et al. | Assessment of pharmacologic area under the curve when baselines are variable | |
| Florián et al. | Calculations of hydration entropies of hydrophobic, polar, and ionic solutes in the framework of the Langevin dipoles solvation model | |
| Tornøe et al. | Stochastic differential equations in NONMEM®: implementation, application, and comparison with ordinary differential equations | |
| Jena et al. | Mechanisms of formation of 8-oxoguanine due to reactions of one and two OH• radicals and the H2O2 molecule with guanine: a quantum computational study | |
| CN111161810B (en) | Free energy perturbation method based on constraint probability distribution function optimization | |
| Qi et al. | Water 16-mers and hexamers: assessment of the three-body and electrostatically embedded many-body approximations of the correlation energy or the nonlocal energy as ways to include cooperative effects | |
| Prado-Prado et al. | Entropy multi-target QSAR model for prediction of antiviral drug complex networks | |
| Kim et al. | Inferring gene regulatory networks from temporal expression profiles under time-delay and noise | |
| Chiappino-Pepe et al. | Integration of metabolic, regulatory and signaling networks towards analysis of perturbation and dynamic responses | |
| Kim et al. | Population pharmacokinetics of unbound mycophenolic acid in pediatric and young adult patients undergoing allogeneic hematopoietic cell transplantation | |
| Ling et al. | Population pharmacokinetics of mycophenolic acid and its main glucuronide metabolite: a comparison between healthy Chinese and Caucasian subjects receiving mycophenolate mofetil | |
| WO2008116495A1 (en) | Method and apparatus for the design of chemical compounds with predetermined properties | |
| Green et al. | Evaluation of concomitant antiretrovirals and CYP2C9/CYP2C19 polymorphisms on the pharmacokinetics of etravirine | |
| Close et al. | Ionization energies of the nucleotides | |
| Nandy et al. | Reconstruction of the yeast protein-protein interaction network involved in nutrient sensing and global metabolic regulation | |
| Wang et al. | The QSAR study of flavonoid-metal complexes scavenging OH free radical | |
| Agnihotri et al. | Formation of 8-nitroguanine due to reaction between guanyl radical and nitrogen dioxide: catalytic role of hydration | |
| LaPointe et al. | A review of density functional theory quantum mechanics as applied to pharmaceutically relevant systems |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |