CN102930152B - A kind of ligand molecular of simulating reacts with target receptor and calculates the method and system of thermodynamics and dynamics parameter predicting this reaction - Google Patents

Abstract

The present invention relates to a kind of react with target receptor based on computer simulation ligand molecular and calculate predict the thermodynamics (Thermodynamics) of this reaction and the method for kinetics (Kinetics) parameter and application system.Specifically include that 1) design the method calculating ligand molecular with target receptor Conjugated free energy, devise multiple target molecular docking method based on the method and develop its application system；2) use this system that simulation and sampling, the acquisition Conjugated free energy of ligand molecular with target receptor Interactions Mode are marked, build the method system of ligand molecular and the Conjugated free energy potential energy level of target receptor according to analog sampling result；3) method system of ligand molecular and target receptor identification response path is obtained according to potential energy level；4) determine that the avtive spot of reaction and transition state site are to calculate thermodynamics and dynamics parameter.The present invention provides important evidence and prediction and evaluation means for medicament research and development, is expected to improve drug development efficiency and significantly save experimental expenses, becomes effective aid of medicament research and development.

Description

A kind of ligand molecular of simulating reacts with target receptor and calculates the method for the thermodynamics and dynamics parameter predicting this reaction and be System

Technical field

The invention belongs to medicine, calculate chemistry, calculation biology and drug design field, relate to a kind of based on computer simulation Ligand molecular (such as, drug molecule) and the reaction of target receptor (such as disease targets such as albumen, nucleic acid), calculate to obtain and be somebody's turn to do The new method of the thermodynamics and dynamics parameter of reaction.Described method can accurate simulation ligand molecular and target receptor at molecule The model of action of level and reaction (in conjunction with/dissociate) mode, and calculate the thermodynamics obtaining this reaction quickly and accurately with dynamic Mechanics parameter, obtains combination or Dissociation path that ligand molecular reacts with target receptor, and then excavates target receptor-ligand molecular Pivotal role site (avtive spot and transition state site), it is adaptable to quantitatively, qualitative analysis ligand molecular and target receptor Between Interactions Mode and effect power.Foundation is provided for designing or transform higher active drug molecule, novel for screening Compound provides quantitative thermodynamics and dynamics to characterize, and improves the success rate finding reactive compound.

Background technology

Molecular recognition has vital effect in biosystem.Ligand molecular and the Interactions Mode meeting of target receptor Cause a series of physiology and pharmacological effect (Charles R.Cantor PRS (1980) Biophysical Chemistry:Part III:The Behavior of Biological Macromolecules(W.H.Freeman,New York),pp 849.Held M, NoéF(2012)Calculating kinetics and pathways of protein–ligand association.Eur J Cell Biol 91:357-364.), this identification process includes ligand molecular and the combination of target receptor or dissociation reaction.Particularly, due to many The drug effect of number medicines depend on part and corresponding target directly in conjunction with process, then drug design is then with part and target Mark receptor binding isotherm is as basic foundation.In consideration of it, drug design with the binding mechanism of receptor is then by understanding part in detail Vital (Taft CA, DA Silva VB, Da Silva CH (2008) Current topics in computer-aided drug design.J Pharm Sci 97:1089-1098.).Traditional drug design process is simply paid close attention to ligand molecular and target receptor The size of binding affinity, generally believe the Molecule Action that molecule that affinity is strong can be more weak than affinity higher (Zhang R, Monsma F(2010)Binding kinetics and mechanism of action:toward the discovery and development of better and best in class drugs.Expert Opin Drug Discov 5:1023-1029.).But this The drug design method success rate that kind is based purely on affinity is extremely low.Trace it to its cause, be owing to this method for designing have ignored instead Answer the kinetic property of ligand molecular and target receptor cohesive process.Substantially, drug effect is not only with thermodynamic parameter (such as: knot Close free energy (Δ G_binding), binding constant (K_A) and dissociation constant (K_D)) relevant, also with binding kinetics parameter (such as: combine speed Rate constant (k_on), kinetic constant (k_off), the holdup time (1/k of the combination of drug molecule and target receptor_off)) relevant (Kenakin T,Jenkinson S,Watson C(2006)Determining the potency and molecular mechanism of action of Insurmountable antagonists.J Pharmacol Exp Ther 319:710-723), the holdup time is i.e. regarded as medicine and divides The onset persistent period of son.A lot of examples prove that drug effect and drug molecule there is no direct phase with the binding affinity of target receptor Guan Xing, and present stronger dependency (Swinney DC (2004) with the holdup time of drug molecule Yu the combination of target receptor Biochemical mechanisms of drug action:what does it take for success？Nat Rev Drug Discov 3: 801-808.), therefore, drug design method based on affinity is not comprehensive, and focuses on the work of drug molecule and receptor simultaneously Then obtain paying attention to widely (Zhang R, Monsma F (2010) Binding by the thermodynamics and dynamics method of process kinetics and mechanism of action:toward the discovery and development of better and best in class drugs.Expert Opin Drug Discov 5:1023-1029.)。

Although the kinetic property that medicine (part) molecule is combined with target receptor is at the importance quilt of medicine (part) design field Disclose out, but in the last thirty years, the most do not develop can quickly, Accurate Prediction part and receptor binding kinetics The method of matter.In addition, although in the past research in, emphasize always affinity (i.e. macroscopic property) calculating prediction Method, existing thermodynamic calculation method there is also certain defect.Although have developed such as free energy perturbation (FEP) The method that this calculating accuracy is higher, but owing to its calculating time is oversize and make its application there is great limitation. On the contrary, the scoring functions used during molecular docking then has relatively quick advantage, but this method is due to greatly The means such as the employing approximation of amount, simplification so that accuracy in computation greatly reduces.Between both approaches MM-GBSA/PBSA has advantage (Wang, J.M. more accurately and faster；Hou,T.J.；Xu,X.J.(2006) Recent advances in free energy calculations with a combination of molecular mechanics and continuum models.Curr Comput-Aided Drug Des 2:287–306.).The method had obtained widely in the last few years Application, but this method is proved and often has system dependency, although can predict accurately for some system, but It is it cannot be guaranteed that all systems are predicted the most exactly.

In the last few years, researcher it is also proposed a kind of method based on molecular dynamics structure and combined potential energy level or Computational Thermodynamics With the scheme of kinetic property, the program first passes through employing molecular dynamics simulation sampling the most thousands of kinetics rails up to a hundred Mark, then by these tracks, obtains ligand receptor kinetic parameter according to corresponding theory computational methods.This method by In need to be by a large amount of sample track of dynamic method, it is therefore desirable to expend substantial amounts of calculating resource and time, prediction simultaneously obtains Being also inaccurate of thermodynamic parameter and kinetic parameter.As Buch et.al. carrys out structure by 495 kinetic locus of sampling Build zamidine and tryptic interaction potential energy face, thus predict macroscopic property more accurately according to potential energy level (with examination Test phase ratio 1kcal/mol), but calculated kinetic property then has at least an order of magnitude compared with test Deviation (Buch I, Giorgino T, De Fabritiis G (2011) Complete reconstruction of an enzyme-inhibitor binding process by molecular dynamics simulations.Pro Natl Acad Sci USA 108: 10184-10189.)。

Therefore, it is developed to predict thermodynamic parameter and the kinetics that ligand molecular target receptor reacts accurately and rapidly The method of parameter is still that problem demanding prompt solution.

Summary of the invention

Technical problem

The technical problem to be solved in the present invention is: simulation ligand molecular and the interaction mode of target receptor, obtains part and divides Son and target receptor (in conjunction with or dissociate) mode reacted, calculate thermodynamic parameter and the kinetic parameter of this reaction, determine The key binding sites of experience when ligand molecular reacts with target receptor, and then instruct design or the transformation of drug molecule, simultaneously Overcome simulation method in prior art time-consumingly long and inaccurate defect, improve drug development efficiency, reduce R&D costs.

Specifically, the purpose of the present invention is: the nature interacted under field conditions (factors) according to ligand molecular and target receptor Rule, simulation and the Interactions Mode that systematically sampling ligand molecular is possible with target receptor, quantify each binding mode Active force (i.e. Conjugated free energy), accurately and comprehensively builds the combination of ligand molecular and target receptor Interactions Mode certainly By energy potential energy level (Binding Free Energy Landscape, BFEL), and calculate pre-quickly and accurately according to potential energy level Survey ligand molecular to be combined with target receptor or thermodynamic parameter, kinetics in dissociation reaction process and this course of reaction is joined Number；Obtain and determine that ligand molecular combines or the key binding sites on path that dissociates on target receptor, for ligand molecular, especially It is that design or the transformation of drug molecule provides new approaches and effective foundation.

Technical scheme and beneficial effect

Therefore, the invention provides the model of action of accurate simulation and exhaustive ligand molecular and target receptor, explore ligand molecular The reaction of target receptor (in conjunction with and dissociation reaction), and calculate the new method of the thermodynamics and kinetics parameter obtaining this reaction And application system.For reaching this purpose, the present invention according between ligand molecular and target receptor Interactions Mode from So rule, it is proposed that the simulation of improvement and assess the molecular docking method of this Interactions Mode mode with The method of MM-GBSA/PBSA, and merge that both approaches carrys out system build complete ligand molecular and target receptor phase The BFEL of interaction, and then be combined or dissociation reaction path with target receptor using the BFEL obtained as prediction ligand molecular Basis, use the response path method of sampling obtain combine or dissociation reaction path, thus be calculated corresponding thermodynamics ginseng Number and kinetic parameter.

Specifically, method of the present invention is to use computer program to realize the scheme that external data processes, external data For have art-recognized meanings technical data (as in the present invention process part (medicine) molecule, and target receptor molecule), Processing procedure is according to the natural law of intermolecular interaction pattern, and described method is capable of simulating quickly and accurately and commenting Valency drug molecule and disease targets model of action, thus calculate thermodynamic parameter and the kinetics of intermolecular interaction pattern Parameter and then instruct SARS drug design or transformation and reach beneficial effect.Therefore, method of the present invention is employing technology Means, realize the technical scheme of technique effect.

Described method calculates prediction ligand molecular by following five key steps (as shown in Figure 1) and is combined with target receptor Or thermodynamic parameter when dissociating and kinetic parameter.Wherein, described thermodynamic parameter is: Conjugated free energy (Δ G_binding)、 Binding constant (K_A), dissociation constant (K_D)；Described kinetic parameter is: combine and dissociate free energy of activation Association rate constant and dissociation rate constant (k_on, k_off):

Step S1: from online target receptor crystal structure database, such as RCSB (http://www.pdb.org/pdb/home/home.do) download by laboratory facilities (such as X-ray diffraction and nuclear-magnetism diffraction) The three-dimensional structure data of the target receptor studied obtained, measures ligand molecular and the space structure data of target receptor, base The simulated system of ligand molecular and target receptor built by space structure in ligand molecular and target receptor；

Step S2: build the BFEL of ligand molecular and target receptor；

Step S3: obtain ligand molecular according to BFEL, with target receptor, combination or the response path dissociated occur；

Step S4: calculate corresponding thermodynamics and kinetics parameter according to the response path obtained；

Step S5: obtain the mechanism of action of molecule reaction by response path, determine the pass bond of decision reaction kinetic characteristics Close site.

Wherein, described ligand molecular is Medicine small molecule or compound, and affiliated target receptor is disease-associated protein, this Bright embodiment uses the huperzine A little molecule of part as an example, and acetylcholinesterase is as target receptor.

Wherein, described step S1 is: build ligand molecular and target receptor simulated system, i.e. to intend research ligand molecular, The pre-treatment step of target receptor and definition simulated system spatial dimension size etc..Described step S1 comprise the following steps (as Shown in Fig. 2):

Step S1.1: definition research system bulk scope.Particularly as follows: the three of the target receptor obtained by assessment experiment Dimension structure delimit the bulk scope intending sampling ligand molecular and target receptor Interactions Mode, need here to point out Be, for this size range division with define, depend on the three dimensional structure size of target receptor, simultaneously in order to ensure Can obtain ligand molecular binding mode in the solution and active force, therefore this spatial dimension should be fully big, such that it is able to carry For enough solution space for ligand molecular in the case of not disturbed by the active force of target receptor with solution reaction.Separately Outward, if this bulk crosses conference results in the most meaningless calculating.In view of the foregoing, it is proposed that this space Size should be bigger than target crystal three dimensional structure 10～20 angstroms, and this numeral there is no harsh defining certainly, other size, such as 20～ ∞ should be protected, the biggest, is unalterable solution region such as more than the spaces of 20 angstroms, there is no calculating It is worth；

Step S1.2: the research space delimited is carried out subspace stress and strain model.Particularly as follows: to described in the definition of step S1.1 Research system space carries out subspace stress and strain model, and is marked every sub spaces, and the present invention is with l_iCarry out labelling, represent I-th subspace, is divided into several little spaces, referred to as subspace by the whole research space delimited.This division purpose has Two: 1) whole simulation can be calculated procedure decomposition, both put into ligand molecular in every sub spaces, this subrange of sampling Interior ligand molecular and the binding mode of target receptor, this mode is so that task is decomposed, it is simple to parallel computation, thus Increase substantially computational efficiency；2) avoid because optimization method convergence and cause missing important ligand molecular and target receptor Interactions Mode, as determined the Interactions Mode of ligand molecular and the target receptor of the transition state of kinetic property.Cause Optimization method for using in docking sampling process can make to sample the Interactions Mode obtained converge to global optimum (in conjunction with Free energy is minimum) or local optimum (Conjugated free energy is relatively low), if not carrying out local limit sampling, the then phase obtained Interaction pattern will necessarily converge at the avtive spot of target receptor that (at this, part with receptor Conjugated free energy theory is substantially The overall situation is minimum), and determine that the Interactions Mode at the transition state that the Conjugated free energy of kinetics character is higher can be wrong Lose, so that the reactive kinetics parameters obtaining ligand molecular with target receptor cannot be calculated；

Described step S1.2 have the beneficial effect that Subspace partition strategy makes calculating process can increase substantially with parallelization Time, also can be prevented effectively from simultaneously and cause missing important Interactions Mode, as determined kinetics because of optimization method convergence The transition state Interactions Mode of character.

Step S1.3: the crucial subspace of screening.It is specially the space to obtaining in step S1.2 to screen, by this space The subspace bigger with ligand molecular effect probability is picked out, and the subspace eliminating that probability is less；Described step S1.3 Have the beneficial effect that and avoid insignificant calculating, waste calculates resource.

Described step S1.3 comprises the following steps:

Step S1.3.1: from online target receptor crystal structure database, such as RCSB (http://www.pdb.org/pdb/home/home.do) download and obtained by laboratory facilities (such as X-ray diffraction and nuclear-magnetism diffraction) The three-dimensional structure data of the target receptor studied obtained, described crystal structure includes that apo (is i.e. only target receptor monomer knot Structure, coexists without little molecule) and holo structure (the target receptor composite structure i.e. coexisted with little molecule).Owing to target is subject to Body structure has behavioral characteristics, the most under difficult environmental conditions, or in the case of combining different little molecule, and its three-dimensional conformation Can there is certain difference, then the channel configurations reacted from ligand molecular accordingly is likely to different, therefore all by enrichment Known crystal structure, can avoid the skewed popularity in the result of calculation brought by the three dimensional structure of target receptor；

Step S1.3.2: use and combine for the little molecule of part on target receptor or the third party software of dissociation channel calculating sampling, The passage t that each target receptor crystal structure is potential is calculated such as MOLE or CAVER etc._i, and it is enriched with all passage T=∑ t_i, It is characterized in that, use third party software, passage t potential on each target receptor of sampling successively_i, each target the most at last The passage obtained on receptor is integrated, i.e. T=∑ t_i, (here, i represents the i-th target receptor three dimensional structure calculated, I depends on the required three dimensional structure number calculating target receptor of user, for be more than 1 positive integer；t_iFor being obtained by calculating Passage potential in the i-th structure obtained；MOLE or CAVER that the present invention enumerates can be wished with self-defined each calculating Hope the most possible number of active lanes obtained, be defaulted as 3.).

Step S1.3.3: cluster all critical passage T and obtain representative passage T_r, (concrete: to analyze in S1.3.2 and obtain All passage T, the passage of each orientation only retains one, and other filter for redundancy.)；

Step S1.3.4: select all coverings or neighbouring T_rSubspace l_iEmpty as joint mode sampling study portion key Between.Particularly as follows: according to these passages T_rS.1.1 with S.1.2, key is selected in all subspaces of obtaining from step The space that subspace is implemented as next step Interactions Mode sampling operation.It is of course possible to select all spaces, and select Method can also suitably adjust, as comprised other subspaces, potential binding site place outside the passage on target receptor, these All should be included among this patent scope；

Step S1.3.5: described crucial subspace is carried out BORDER PROCESSING, it is characterised in that interleaving of adjacent subspace Enter to cover the extra subspace l on border_j, carry out labelling with j, insert Σ l altogether_jThe extra subspace of height covers border；

Step S1.3.6: the final crucial subspace L=Σ l obtaining joint mode sampling_i∪Σl_j。

Described step S1.3.1 provides the benefit that: in view of the different three dimensional structures of target receptor, it is to avoid because of architectural difference band The calculating skewed popularity come, misses important results during i.e. avoiding channel sample in the next step.

Described S1.3.5 has the beneficial effect that the ligand molecular-target receptor effect preventing from missing boundary key that may be present Pattern.

Described step S2 includes step in detail below, as shown in Figure 4:

Step S2.1: in every sub spaces l_i∈ L puts into a ligand molecular, docks, every sub spaces l_iObtain one Individual target receptor-ligand molecular Interactions Mode set P_i=∑ p_k(k is the positive integer more than or equal to 0 here, p_kRepresentative obtains Kth target receptor-part the Interactions Mode obtained), particularly as follows: use molecular docking method analog sampling ligand molecular At all of Interactions Mode of this subspace Yu target receptor, for the interaction mould obtained accurately and multiformity is abundant Formula, the present invention designs and develops efficient, a novel molecular docking method, this docking calculation 1) right in sampling process The prediction of binding mode and evaluation calculate function (MM-GBSA/MM-PBSA) with the Conjugated free energy improved herein and measure Degree, each energy term in this function is recombinated according to feature and considers target receptor conformation change and the relative free energy change that causes Change；2) simultaneously according to these free energy computational methods, multiple target molecular docking method is developed in conjunction with Multipurpose Optimal Method；3) For taking into account accuracy and the computation rate of docking result, target receptor is carried out the division of rigid and flexible region and processes.Beneficial effect For: improve the accuracy of prediction ligand molecular and target receptor binding mode, improve and evaluate each Interactions Mode The accuracy of Conjugated free energy；

Specifically, 1) improvement to free energy function is characterised by: each energy term is reassembled as four or multiple according to feature Display items: such as Δ E_vdw,ΔE_es,ΔG_gb,ΔG_sa,；Describing the problem for apparent, the present invention illustrates as a example by above-mentioned four, Reorganization Energy quantifier also can be not only limited to above-mentioned four kinds, as long as split based on energy term, and the multiple target docking calculation designed Belong to scope

${\mathrm{\ΔG}}_{binding}^o={\mathrm{\ω}}_1{\mathrm{\ΔE}}_{vdw}+{\mathrm{\ω}}_2{\mathrm{\ΔE}}_{es}+{\mathrm{\ω}}_3{\mathrm{\ΔG}}_{gb}+{\mathrm{\ω}}_4{\mathrm{\ΔG}}_{sa}+...\left(1\right)$

${\mathrm{\ΔE}}_{vdw}={\mathrm{\ΔE}}_{vdw,RL}+{\mathrm{\ΔE}}_{vdw,L}^{conf}+{\mathrm{\ΔE}}_{vdw,R}^{conf}---\left(2\right)$

${\mathrm{\ΔE}}_{es}={\mathrm{\ΔE}}_{es,RL}+{\mathrm{\ΔE}}_{es,L}^{conf}+{\mathrm{\ΔE}}_{es,R}^{conf}---\left(3\right)$

${\mathrm{\ΔG}}_{gb}={\mathrm{\ΔG}}_{POL}+{\mathrm{\ΔG}}_{pol,R}^{conf}=(G_{pol,RL}-G_{pol,R}-G_{pol,L})+{\mathrm{\ΔG}}_{pol,R}^{conf}---\left(4\right)$

${\mathrm{\ΔG}}_{sa}=\[{\mathrm{\σ}}_1\mathrm{\Δ}\left({\mathrm{SA}}_{hp}\right)-{\mathrm{\σ}}_2\mathrm{\Δ}\left(SA\right)\]+\[{\mathrm{\σ}}_1\mathrm{\Δ}\left({\mathrm{SA}}_{hp}^{conf}\right)-{\mathrm{\σ}}_2\mathrm{\Δ}\left({\mathrm{SA}}^{conf}\right)\]---\left(5\right)$

ΔE_vdwFor van der Waals energy quantifier, Δ E_esFor electrostatic energy quantifier, Δ G_gbFor solvent polarization free energy, Δ G_saNon-polarized for solvent Free energy；ω₁-ω_nWeight for each energy term；RL is the interaction energy term of target receptor ligand molecular, and R is The energy term of target receptor self, L is the energy term of ligand molecular self；ΔE_vdw,RLFor target receptor and ligand molecular Conformation changes when identifying and causes the variable quantity of the interaction Van der Waals energy term of target receptor ligand molecular, ΔE_es,RLConformation for target receptor with ligand molecular changes when identifying and causes the phase of target receptor ligand molecular The variable quantity of interaction electrostatic term；It is to be sent out when identifying by the conformation of target receptor with ligand molecular respectively Changing and the van der Waals energy quantifier of target receptor that causes and the variable quantity of electrostatic term；It is by target respectively The van der Waals energy quantifier of the ligand molecular that the conformation of mark receptor and ligand molecular changes when identifying and causes and electrostatic term Variable quantity；With Δ (SA^conf) it is to be occurred when identifying by the conformation of target receptor with ligand molecular The solvent polarity changed and causeWith nonpolar free energy (With), σ₁=0.025, σ₂=0.015 is the conversion constant of area and energy unit kcal/mol.

2) the multiple target docking calculation designed is characterised by: based on Multipurpose Optimal Method (such as multi-objective genetic algorithm, Multi-objective Evolutionary Algorithm and multi-objective particle swarm algorithm etc.) with free energy computational methods and the flexible docking method that designs.With Tradition docking calculation is compared, and the having the beneficial effect that of this docking calculation can increase substantially the Interactions Mode that sampling obtains Multiformity, and improve the accuracy of sampling configuration and the evaluation accuracy of corresponding Conjugated free energy.Wherein, for object function Structure include following two scheme:

By Δ E_vdw,ΔE_es,ΔG_gb,ΔG_sa... as docking multiple targets or fitness function (each energy term definition see Claim 6, can also have other energy term to expand this multiple objective function system further simultaneously)；Or

By Δ E_intra,ΔE_inter,ΔG_gb,ΔG_sa... as multiple targets or the fitness function of docking；

Wherein, Δ E_vdwFor van der Waals energy quantifier, Δ E_esFor electrostatic energy quantifier, Δ G_gbFor solvent polarization free energy, Δ G_saFor The non-polarized free energy of solvent；ΔE_intraFor target receptor, the energy term of ligand molecular self, Δ E_interFor target receptor part Non-bonded interaction between molecule, and the target numbers employing elasticity setting system that optimization process is used, i.e. select between 1 n Select.

3) rigidity is used to be characterised by with flexible region division process target receptor during docking: to exist at simulation ligand molecular In described subspace and during all Interactions Modes of target receptor, target receptor is used rigidity and flexible region Division processes, and wherein, flexible region considers target receptor part-structure transmutability, as protein Key residues or domain (as Loop district) motion change, and use energy scoring tactics based on mesh point Yu atom pair process described flexible region and Rigid region.Flexible region is it is contemplated that protein Key residues or the motion change of ring region, and uses based on mesh point and atom To energy scoring tactics process both regions.Have the beneficial effect that and improve calculating speed；Simultaneously for avoiding because of rigidity The calculating error processed with flexible region division and use two kinds of scoring tactics to be brought, Interactions Mode is unified the most at last adopts Process with the scoring again of lattice point scoring tactics and solve this hidden danger.

Step S2.2: collect the receptor-ligand Interactions Mode in all subspaces and obtain P=∑ P_i=∑ ∑ p_k.Particularly as follows: The ligand molecular of all subspaces-target receptor Interactions Mode enrichment is obtained Interactions Mode set P=∑ P_i=∑ ∑ p_k, as the data element building BFEL；The system of having the beneficial effect that is comprehensively reflected ligand molecular and target receptor All possible model of action；

Step S2.3: the definition coordinate axes that reacts with target receptor of ligand molecular, mutual by each ligand molecular-target receptor Binding mode p_k(x, y) carries out projection to the reaction coordinate axle reacted at ligand molecular-target receptor, each interaction mould Formula p_k(x, y) combines potential energy level as target receptor-ligand molecular and builds element, and (reaction coordinate axle can be various dimensions, simplifies For the sake of, the present invention is as a example by two dimension, and two-dimentional reaction coordinate axle x, y as defined herein, such as p_k(x y) represents S2.2 step The Interactions Mode p obtaining ligand molecular-target receptor calculated_kIn x, the projection of y-axis).Have the beneficial effect that reflection is joined Body molecule is combined with target receptor or the instantaneous conformation change of the progress dissociated and part and target receptor, such that it is able to clearly Reflect intermolecular mechanism of action；

Step S2.4: by the three-dimensional of third party's three-dimensional drawing Software on Drawing ligand molecular Yu target receptor Interactions Mode BFEL, wherein, x-axis and y-axis are reaction coordinate axle, and z-axis is Conjugated free energy coordinate axes.

Described step S3 farther includes two kinds of strategies:

A () does not consider that the reaction backtracking phenomenon of ligand molecular and target receptor, path based on response path minimum energy principle are searched Rope strategy；

B () considers ligand molecular and the reaction backtracking phenomenon of target receptor and the Direction of Reaction random disturbance, based on response path energy Measure the route searching strategy of minimum principle.

Strategy (a) comprises the following steps, as shown in Figure 5:

Step S3.1.1: coarseness divides potential energy levelΩ represents whole BFEL,Represent whole respectively Individual potential energy level is according to the lattice point number of the upper division of reaction coordinate axle x, y defined in step S2.3, and n, m are respectively lattice point interval Size, determines reaction coordinate axle x, y upper lattice point number；Particularly as follows: the most more complicated owing to obtaining potential energy level, thus add Big calculating intractability, is therefore quick obtaining response path, first original potential energy level is carried out simplification process, the most right BFEL carries out coarseness divisionMake potential energy dimensionality reduction, thus reduce computation complexity, have the beneficial effect that and carry Computationally efficient；

Step S3.1.2: selected response path original positionx_iRepresent in reaction coordinate x I-th lattice point, y_jRepresent the jth lattice point in reaction coordinate y.Initialisation reactions path S:=S ∪ S₀；Concrete: choosing Determine ligand molecular and the dissociation reaction of target receptor or the original position of association reaction, if the original position of dissociation reaction can be gesture Energy face overall situation minimum point, association reaction original position can be part and the certain point in solution effects region；

Step S3.1.3:S_p=min{S_p(x_i-1,y_j), S_p(x_i-1,y_j+1), S_p(x_i,y_j+1),S_p(x_i+1,y_j+1),S_p(x_i+1,y_j)}；S:=S ∪ S_p+1, On the BFEL of coarseness, it is iterated search according to minimum energy principle and obtains the reaction road of ligand molecular and target receptor Footpath, wherein each present position of ligand molecular all faces 5 advanceable directions, next step moving direction of ligand molecular Being the minimum direction of Conjugated free energy in 5 directions, wherein S is for finally to obtain path, S_pRepresent ligand molecular to be subject to target The pth step of body node location in the potential energy level that coarseness divides, p is the positive integer more than or equal to 0；

Step S3.1.4: judge present position S_PWhether arrive user-defined reaction end (if as research ligand molecular with The association reaction of target receptor, then reaction end should be the catalysis on target receptor or the avtive spot of binding partner molecule；If Research ligand molecular and the dissociation reaction of target receptor, then reaction end should be the region away from target receptor, generally solvent Region.), it is then to enter step S3.1.5, otherwise proceeds to step S3.1.3；

Step S3.1.5: the ligand molecular that iteration is obtained and the coordinate information of target receptor response path S (the most corresponding x, y The position of axle) backtracking is to the BFEL not carrying out coarseness division, it is thus achieved that ligand molecular and target receptor response path and combination Functional relationship f (S) of free energy, according to response path coordinate, revert to totipotency face and (does not i.e. carry out coarseness division BFEL), thus Regeneration response path corresponding potential energy information: f (S)；

Strategy (b) comprises the following steps, as shown in Figure 6:

Step S3.2.1: BFEL is carried out coarseness divisionWith step S3.1.1；

Step S3.2.2: with step S3.1.2, selected response path original position,Just Beginningization response path S:=S ∪ S₀；

Step S3.2.3: judge present position S_PWhether arrive defined reaction end, be then to enter step S3.2.12, no Then proceed to step S3.2.4；

Step S3.2.4: according to neighbor node potential energy numerical value ascending sort S_pNeighbor node { S_p(x_i-1,y_j),S_p(x_i-1,y_j+1), S_p(x_i-1,y_j-1),S_p(x_i,y_j-1),S_p(x_i+1,y_j-1),S_p(x_i+1,y_j),S_p(x_i+1,y_j+1),S_p(x_i,y_j+1) (x, y are with step S2.3), return joint Point ordered queue min_S_vec, prepare for subsequent reactions track search, enter step S3.2.5；

Step S3.2.5: check each present position S of part target_p(x_i,y_i) it is whether the site, border of potential energy level, if not Then proceed to step S3.2.6, if then proceeding to step S3.2.7；

Step S3.2.6: the present position judging part is which kind of the geomorphic feature site in potential energy level, these geomorphic features master It is divided three classes: potential well position, position, mountain peak and climbing position, if instant site is potential well position or position, mountain peak, I The random number that produces according to program determine next step moving direction of ligand molecular, random number uses [0,1] random function Producing, if value is in [0,0.9], in the case of the biggest probability, ligand molecular is according to the ordered nodes returned in step S3.2.4 Sequence min_S_vecSelect, proceed to step S3.2.11；If random number value is in [0.9,1.0], i.e. small probability situation Under, the reaction direction of advance of ligand molecular, for randomly choosing, i.e. proceeds to step S3.2.12；

Step S3.2.7: judge S_p(x_i,y_i) adjacent node in potential energy value minimum position whether already contained in ligand molecular and target In mark receptor response path S, if this position is not in response path S, then proceed to step S3.2.8；If this position is It is included in and explores in the response path obtained, and this position occurs once, then for avoiding path in response path incessantly Heuristic process is absorbed in endless loop because hovering in this position, is limited by the number of times that this position occurs in response path.As Really some position occurs M time in response path, and we are considered as this position can cause endless loop, and (M is for being more than The positive integer of 0, can be defined depending on tolerance by user, be closer to 1, superfluous in the most whole response path exploring acquisition Remaining position is the fewest, reflects ligand molecular and " roundabout " feature of generation in target receptor cohesive process more than 1, It is set as 3 in embodiments of the present invention).If this position is also unsatisfactory for the above-mentioned definition causing endless loop, then proceed to step S3.2.8；If meeting endless loop definition to proceed to step S3.2.9；

Step S3.2.8: by present node S_pNeighbor node ordered queue min_S_vecIn the minimum node location of potential energy min_S_vec[j] is as next step moving direction (S:=S ∪ min_S_vec[j]), more new route present position (p=p+1), proceed to step Rapid S3.2.3；

Step S3.2.9: having neighbor node ordered queue min_S at random_vecOne direction of advance (min_S of middle selection_vec[j]), More new route present position (S:=S ∪ min_S_vec[j], p=p+1) and proceed to step S3.2.3；

Step S3.2.10: if S_p(x_i,y_j) be the climbing node in potential energy level, then with S_p(x_i,y_j) it is that boundary is by orderly for neighbor node team Row min_S_vecTwo points, free energy numerical value is more than S_p(x_i,y_j) node composition sequence be L_S_vec；Less than S_p(x_i,y_j) node Composition sequence is S_S_vec, we determine next step moving direction of ligand molecular according to the random number that program produces, at random Number uses [0,1] random function to produce, and is in all response paths that may be present as far as possible for ensureing the response path obtained Free energy is minimum, in the case of the present invention is set in greater probability (when the random number i.e. obtained is in the range of less one, this Invention is when being set as [0.3,1]) response path direction still according to minimum energy path direction, i.e. proceed to step S3.2.11, S_S_vecSequence selects；Otherwise, produce the effect of the random perturbation of response path, i.e. proceed to step S3.2.12, at L_S_vec In select；

Step S3.2.11: at S_pThe ordered queue min_S of neighbor node_vecOr S_S_vecThe middle position min_S selecting minimum_vec[j] Or S_S_vec[j] as next step the Direction of Reaction of ligand molecular, with step S3.2.7, for avoid track search process because of hover in This position and be absorbed in endless loop, the number of times that this position occurs in response path is limited.Can be come regarding tolerance by user Defining, test case of the present invention is set to 8；If this position is unsatisfactory for the above-mentioned definition causing endless loop, proceed to step S3.2.12, Otherwise using this position as next step reaction site, Regeneration response path S (S:=S ∪ min_S_vec[j] or S:=S ∪ S_S_vec[j], P=p+1)；

Step S3.2.12: at the neighbor node ordered queue min_S of current location_vecOr L_S_vecIn randomly choose a node min_S_vec[j] as next step the Direction of Reaction of path, more new route present position S (S:=S ∪ min_S_vec[j] or S:=S ∪ L_S_vec[j], p=p+1), and proceed to step S3.2.3；

Step S3.2.13: the ligand molecular that iteration is obtained and the coordinate information of the response path S of target receptor revert to into The BFEL that row coarseness divides, it is thus achieved that ligand molecular and target receptor response path and functional relationship f (S) of Conjugated free energy: That is, the response path S that will determine, according to its response path coordinate information, revert to the BFEL not having coarseness to process, Obtain accurate response path and obtain functional relationship: f (S) with Conjugated free energy, carry for calculating the thermodynamic power of reaction generation For foundation；

Step S4: calculate corresponding thermodynamics, kinetic parameter according to the response path obtained.In this step, first root According to the response path of the ligand molecular obtained with target receptor, determine the avtive spot (f (S) in course of reaction_{Minimum position}), transition State site (f (S)_{Maximum position}) and reaction terminating site (being defined by step S3.1.4).Conjugated free energy according to these three site Data △ G_bind,,△G_trans,△G_unbind, obtain Conjugated free energy in conjunction with transition state theory Association reaction activation energyWith dissociation reaction activation energy Thus obtain and obtain Conjugated free energy according to transition state theoryAssociation reaction activation energyWith dissociate Reaction activityAnd then according to Aileen's formula (such as formula 6-8) calculating acquisition thermodynamic parameter: binding constant K_DWith Kinetic parameter: association rate constant k_onWith rate constants k of dissociating_off；

$K_D=\frac{k_{off}}{k_{on}}---\left(6\right)$

${\mathrm{\ΔG}}_{binding}^o={\mathrm{RTlnK}}_D---\left(7\right)$

${\mathrm{\ΔG}}_{on}^{\≠}=-RT\ln \left(\frac{k_{on}h}{k_{B}T}\right),{\mathrm{\ΔG}}_{off}^{\≠}=-RT\ln \left(\frac{k_{off}h}{k_{B}T}\right)---\left(8\right)$

K_BFor Boltzmann constant, R is gas constant, and T is absolute temperature, and h is planck constant.

Having the beneficial effects that of described step S4: Accurate Prediction have rated the combination of drug molecule and target receptor or dissociates instead Answer the difficulty of mode, bond strength and association reaction and dissociation reaction, indirectly obtain drug molecule and target receptor Action time, i.e. (1/k_off)。

Step S5: obtain molecule response mechanism, verify the key binding sites determining reaction kinetic characteristics.Specifically: logical Cross the response path obtained to revert to free energy potential energy level and obtain each present position of ligand molecular on this path S_p(x_p,y_p) reaction coordinate (x_p,y_p), the part being projected in this position i.e. can be found in step S2.3 to divide by reaction coordinate Son and the binding pattern P of target receptor_k(x_i,y_i), x here_p=x_i,y_p=y_i, therefore can obtain ligand molecular active binding site, The binding pattern with target receptor in transition state site, learns the avtive spot on corresponding target receptor and transition state site (binding site on target receptor passage that i.e. ligand molecular is mentioned in step S1.3.4), and then can design and target Receptor has relatively low Conjugated free energy (i.e. at avtive spotLower), and have in transition state site with target receptor The compound of higher Conjugated free energy is (i.e.Higher), this newly-designed compound not only has higher with target receptor Affinity, have simultaneously and there is longer action time, thus there is higher compound activity.

Having the beneficial effects that of described step S5: break through in conventional medicament design and be concerned only with medicine (part) molecule and target Receptor must find or design the one-sidedness of medicine at avtive spot in conjunction with situation, it is provided that determines the transition state position of medicinal effectiveness Point, provides new key message and design considerations for drug design.Specifically, it is thus achieved that the avtive spot on target receptor And transition state site is for instructing medicine (part) design to play the effect of key；Show: can be by rough evaluation Method probes into other ligand moleculars with target receptor in the Conjugated free energy situation of avtive spot, can also evaluate part simultaneously Molecule and target receptor are in the Conjugated free energy situation of transition state, such that it is able to by calculating drug molecule in both sites Binding ability thus obtain roughly the kinetic property of ligand molecular, and then can reduce it with receptor in active sites by design The Conjugated free energy that point combines, improves its Conjugated free energy being combined with receptor in transition state site to improve the work of compound Property.Thus overcome and conventional medicament (part) design is confined to avtive spot Conjugated free energy and ignores its kinetic property Shortcoming.Obtain the mechanism of action of molecule reaction, verify the key binding sites determining reaction kinetic characteristics, (join as medicine Body) design key reference factor.

According to said method, the present invention is its application system based on computer development, and described system includes:

First module, it performs described step S1, i.e. builds ligand molecular and target receptor simulated system, according to by experiment The crystal three dimensional structure of the target receptor parsed is to define space structure and the size model of the possible effect of ligand molecular and its Enclose the research range being used as simulating ligand molecular with target receptor Interactions Mode, be used for preparing ligand molecular and target The method pretreatment of mark receptor and the ligand molecular and the target receptor that prepare pre-simulated；

Second module, it performs described step S2, builds the BFEL of ligand molecular and target receptor, simulate exhaustive part and divide Son and the Interactions Mode of target receptor also assess the corresponding Conjugated free energy of these binding modes, thus systematically build Ligand molecular and target receptor interaction BFEL；

Three module, it performs described step S3, obtains ligand molecular according to BFEL and occur to combine or dissociate with target receptor Response path (Pathway)；

4th module, it performs described step S4, calculates corresponding thermodynamics and kinetics parameter according to the response path obtained；

5th module, it performs described step S5, determines that ligand molecular is in whole response path by the response path obtained The ligand molecular configuration of each present position, revert to dock the Interactions Mode obtained, it is thus achieved that ligand molecular is in activity Binding site, the structural information in transition state site and the Constellation information of receptor when part is combined in corresponding site, determine target Avtive spot on mark receptor and transition state site, and then when designing drug molecule, can be designed to living with receptor Property site has of a relatively high Conjugated free energy when combining, and has relatively low Conjugated free energy when transition state site combines Molecule, the most not only improve the binding affinity of drug molecule and target receptor, and extend drug molecule and make with target receptor Holdup time, thus improve the biological activity of drug molecule designed.

Accompanying drawing explanation

In order to be illustrated more clearly that the embodiment of the present invention or technical scheme of the prior art, Figure of description will be carried out below Briefly describe.

Fig. 1 is the trunk flow chart of the method for the invention.

Fig. 2 is the flow chart of step S1 (building ligand molecular and target receptor simulated system).

Fig. 3 is the flow chart of step S1.3 (selecting crucial subspace).

Fig. 4 is that (simulation ligand molecular and the model of action of target receptor, build the knot of ligand molecular and target receptor to step S2 Close free energy potential energy level) flow chart.

The backtracking phenomenon existed when Fig. 5 is the reaction not considering ligand molecular and target receptor, based on response path minimum energy The response path searching method flow process (step S3, strategy (a)) of principle.

Backtracking phenomenon that may be present when Fig. 6 is the reaction considering ligand molecular and target receptor, and consider the Direction of Reaction simultaneously Random disturbance, based on the method for searching path flow process (step S3, strategy (b)) that steepest descent is theoretical.

Fig. 7 is that the space interacting the ligand molecular studied and target receptor carries out Subspace partition and possible part divides The combination of son or the schematic diagram of dissociation channel sampled result, say as a example by the case study on implementation (HupA-TcAChE) in the present invention Bright.

The combination or dissociation channel that Fig. 8 obtains according to sampling, selects to obtain the result schematic diagram of crucial subspace, with the present invention In case study on implementation as a example by explanation.

Fig. 9 is that the present invention plants calculated all ligand moleculars and the binding mode spatial distribution feelings of target receptor in embodiment Condition, is distributed as the sign schematic diagram represented with the configuration space of ligand molecular.

Figure 10 is to calculate the Conjugated free energy potential energy level three-dimensional diagram that the HupA-TcAChE obtained interacts.

Figure 11-A is the two dimension in the association reaction path by the HupA-TcAChE using response path search strategy (a) to obtain Schematic diagram, (owing to the present invention does not has artificial disturbance in whole simulation process to molecular action mode, thus it is believed that part Molecule is identical with dissociation reaction path with acceptor molecule association reaction path).

Figure 11-B is the Conjugated free energy change curve on the HupA-TcAChE association reaction path obtained in Figure 11-A.

Figure 12 is combined with acceptor molecule for using the ligand molecular that in route searching strategy (b), response path heuristic approach obtains/solves From response path schematic diagram.

Figure 13 is the path that the ligand molecular obtained by prediction is reacted with target receptor, it is thus achieved that target receptor on determine anti- Answer critical sites (avtive spot and transition state site) schematic diagram, wherein the Lycoperdon polymorphum Vitt cartoon of macroscopic property and kinetic property Figure is target proteins TcAChE, and Dark grey articulated model is the Key residues of molecular action little with part, black mallet shape mould Type is ligand molecular HupA.

The critical sites obtained according to Figure 14, and carry out drug design strategies schematic diagram by the two site, i.e. pass through Design so has the most highly active compound: its Conjugated free energy being combined at avtive spot with receptor is relatively low, with The Conjugated free energy that receptor combines in transition state site is higher, thus increases the affinity of drug molecule, increase drug molecule with The target receptor effect holdup time, improve drug effect.

It is embodied as and application mode

Below in conjunction with the embodiment of the present invention and accompanying drawing, the technical scheme in the embodiment of the present invention is clearly and completely retouched State, it is clear that described embodiment is only one embodiment of the invention rather than whole embodiments.Based on the present invention In embodiment, the every other enforcement that those of ordinary skill in the art are obtained under not making creative work premise Example, broadly falls into the scope of protection of the invention.

The invention discloses a kind of heating power being calculated by simulation ligand molecular and target receptor reactive mode and predicting this reaction Learning and the method for kinetic parameter, it is according to the natural law structure of Interactions Mode between simulation ligand molecular target receptor Building BFEL, search ligand molecular and the response path of target receptor, determine possible response path, thus obtain and join accordingly The thermodynamics and kinetics parameter of body molecule and target receptor effect, it is thus achieved that determine on receptor ligand molecular intensity in connection or The avtive spot of action time and transition state site, therefore obtained by design and so have the most highly active compound: Its Conjugated free energy being combined at avtive spot with receptor is relatively low, the Conjugated free energy being combined in transition state site with receptor Higher, thus increase the affinity of drug molecule, increase drug molecule and target receptor effect holdup time.The present embodiment with California electricity precious jade acetylcholinesterase (TcAchE) as target receptor, huperzine A ((-)-HuperizineA, HupA) As ligand molecular.It will be appreciated by those skilled in the art that method of the present invention can be equally used for quickly, Assess exactly or predict other any ligand moleculars and the thermodynamics of reactions of any target receptor and kinetic parameter.It is concrete Embodiment is as described below:

The exhaustive ligand molecular of a kind of simulation disclosed in the embodiment of the present invention and target receptor model of action, it was predicted that ligand molecular and target Mark receptor combination/dissociation reaction, and predict trunk flow process such as Fig. 1 institute of the thermodynamics and kinetics parametric technique calculating reaction Show, including:

Step S1: according to the crystal structure of the target receptor that experiment obtains, definition obtains ligand molecular and may make with target receptor Spatial dimension, and build simulation counting system according to the crystal structure of target receptor and the action space scope of definition.Tool The build process of body is as in figure 2 it is shown, include:

Step S1.1: the bulk scope of definition counting system, definition size should be noted: sufficiently large, not only can cover Receptor three dimensional structure, also provide for simultaneously enough solvent spaces for ligand molecular away from receptor protein not by the feelings of its function influence Under condition, freely-movable in the solution, the waste being simultaneously also unlikely to the biggest and meaningless calculates resource.To grind for this example The TcAChE system studied carefully, usesCube, as shown in Figure 7, owing to TcAChE system is less, Therefore, this size range is enough, can also be other sizes certainly for other systems, depends on research system and calculating needs Depending on Yaoing, also let loose in cube simultaneously, be equally other space forms such as spheroid.

Step S1.2: the calculating space of definition is carried out subspace stress and strain model, as it is shown in fig. 7, be parallelization joint mode Sampling process, avoids optimization method global convergence to miss the pivotal role pattern of intermediates simultaneously, the present invention propose by The space of the research system of definition carries out subspace division and processes, and original cube is divided into a series of small cubes, And each small cubes is a sub spaces, as the enforcement space of a docking sampling.Here, for concrete division Method, will in this exampleCube is divided into isopyknic1000 small cubes, when So concrete small cubes size being not limited to the size of this example, being also not required for small cubes size wants one simultaneously Cause.

Step S1.3: select crucial subspace.Due to above-mentioned steps definition research space in exist some ligand moleculars without Method combine scope (such as receptor residues (base) tightly packed region) or distance receptor solvent space too far, these regions for Acquisition ligand molecular is combined with acceptor molecule or the thermodynamic parameter of dissociation reaction there is no significant contribution, therefore with kinetic parameter Can ignore thus computing cost.Concrete selection scheme is following (such as Fig. 3):

Step S1.3.1: collect all of target receptor crystal structure, for same target receptor due to binding partner molecule Front and back, even combine different ligand moleculars and all may induce its conformation that different change occurs, and different receptor structures The passage for ligand molecular Interactions Mode or binding site as existing then are not quite similar, and therefore, the present invention proposes to lead to The method crossing the many receptor conformation passages of research is avoided missing critical passage or binding site.

Step S1.3.2: (present case obtains the crystal structure of 54 TcAChE altogether for the different target receptor structure that obtains For research, such as table 1.), use third party software, as MOLE (M,P,J,Otyepka M(2007) MOLE:A Voronoi diagram-based explorer of molecular channels,pores,and tunnels.Structure 15:1357-1363.), CAVER (Petrek, M., Otyepka, M., Banas, P., Kosinova, P., Koca, J., Damborsky,J.(2006).CAVER:a new tool to explore routes from protein clefts,pockets and Cavities.BMC Bioinformatics 7:316.) etc. search for important channel t present on each structure_i(represent i-th structure The passage obtained, number is indefinite), it is enriched with all critical passage T=∑ t_i；

Table 1

Table 1 obtains for downloading from albumin crystal data base RCSB when sampling ligand-receptor reaction channel in present example The code name of the crystal structure of all TcAChE, totally 54.

Step S1.3.3: cluster all passage T and obtain representative passage T_r.Owing to there is certain weight between the passage of acquisition Folded property, calculates for simplifying, and obtains passage in this step for these and clusters according to similarity, thus finally gives selection Passage T_r=8, such as accompanying drawing 7；

Step S1.3.4: according to the passage T obtained_rSelect crucial subspace.In the research range of definition, cover or face Nearly passage T_rSubspace l_iThe most optional do crucial subspace, as follow-up Interactions Mode sample space.

Step S1.3.5: process the border between the crucial subspace obtained.For avoiding missing the pass of adjacent subspace boundary Key Interactions Mode, the present invention proposes the extra method improved and process border, subspace, i.e. at the boundary of every pair of subspace Extra additional a series of subspace l_j(size can be different from the initial subspace set, for ease of operation, in present example The size subspaces such as employing) cover border, as extra sample space unit, as shown in Figure 8.

Step S1.3.6:L=Σ l_i∪Σl_j=72 are the subspace chosen as joint mode sampling.

For step S1.3, in order to build potential energy level more fully, can also ignore this step, whole subspaces are all as meter Calculate space.

Step S2: obtain its Interactions Mode according to the natural law that analog sampling ligand molecular and target receptor interact And the Conjugated free energy of Thermodynamic parameters pattern carries out calculating assessment, thus build ligand molecular and target receptor interaction mould The BFEL of formula；This step mainly includes flow process (such as Fig. 4) in detail below:

Step S2.1: at every sub spaces (that is, lattice) L_i∈ L puts into a ligand molecular and docks, and every sub spaces obtains Obtain a target receptor ligand molecular Interactions Mode set P_i=∑ p_k；

For build potential energy level, the multiformity of target receptor ligand molecular Interactions Mode, Interactions Mode accurate The evaluation accuracy of the Conjugated free energy of property and each Interactions Mode is most important.And common docking calculation is owing to beating The rough property of point function and evaluation add when combining energy and and the Interactions Mode multiformity obtained that causes sampling is poor, thus lead Cause above three aspect cannot ensure, thus cannot complete to build this operation of BFEL accurately.

The present invention is directed to these limitations, propose a series of corrective measure, be prevented effectively from these problems, be respectively as follows:

1) improvement to free energy evaluation methodology MM-GBSA or MM-PBSA, including for each energy term according to feature Restructuring and consider target receptor conformation change and the relative Gibbs free that causes, thus improve the most foreseeable accurately Property.In the present embodiment, being described as a example by Conjugated free energy is reassembled as four, energetic terms restructuring number can also lack In four or more than four, and unlimited number, carry out developer molecule docking calculation all using energy term as object function based on this Belong to protection scope of the present invention.

Each energy term is reassembled as four according to feature: Δ E_vdw,ΔE_es,ΔG_gb,ΔG_sa；

${\mathrm{\ΔG}}_{binding}^o={\mathrm{\ω}}_1{\mathrm{\ΔE}}_{vdw}+{\mathrm{\ω}}_2{\mathrm{\ΔE}}_{es}+{\mathrm{\ω}}_3{\mathrm{\ΔG}}_{gb}+{\mathrm{\ω}}_4{\mathrm{\ΔG}}_{sa}---\left(1\right)$

${\mathrm{\ΔE}}_{vdw}={\mathrm{\ΔE}}_{vdw,RL}+{\mathrm{\ΔE}}_{vdw,L}^{conf}+{\mathrm{\ΔE}}_{vdw,R}^{conf}---\left(2\right)$

${\mathrm{\ΔE}}_{es}={\mathrm{\ΔE}}_{es,RL}+{\mathrm{\ΔE}}_{es,L}^{conf}+{\mathrm{\ΔE}}_{es,R}^{conf}---\left(3\right)$

${\mathrm{\ΔG}}_{gb}={\mathrm{\ΔG}}_{POL}+{\mathrm{\ΔG}}_{pol,R}^{conf}=(G_{pol,RL}-G_{pol,R}-G_{pol,L})+{\mathrm{\ΔG}}_{pol,R}^{conf}---\left(4\right)$

${\mathrm{\ΔG}}_{sa}=\[{\mathrm{\σ}}_1\mathrm{\Δ}\left({\mathrm{SA}}_{hp}\right)-{\mathrm{\σ}}_2\mathrm{\Δ}\left(SA\right)\]+\[{\mathrm{\σ}}_1\mathrm{\Δ}\left({\mathrm{SA}}_{hp}^{conf}\right)-{\mathrm{\σ}}_2\mathrm{\Δ}\left({\mathrm{SA}}^{conf}\right)\]---\left(5\right)$

ΔE_vdwFor van der Waals energy quantifier, Δ E_esFor electrostatic energy quantifier, Δ G_gbFor solvent polarization free energy, Δ G_saNon-for solvent Polarization free energy；ω₁-ω₄For the weight of each energy term, obtain according to matching empirical value；RL, R, L are represented as respectively Interactions Mode between target receptor ligand molecular complex, R is the energy term of receptor self, and L is part self Energy term.WithIt is that the conformation by target receptor and ligand molecular is when identifying mutually Change, thus the Van der Waals caused changes with electrostatic energy quantifier；Equally,With Δ (SA^conf) be Changed when identifying mutually by the conformation of target receptor with ligand molecular, thus the solvent polarity caused and nonpolar freedom Energy；For weights omega₁-ω₄Acquisition, the present invention propose in view of it has been reported that thermodynamic data (K_A, K_D) and Crystal structural data (such as table 2) carrys out the Fitting Calculation result and obtains, respectively ω₁=0.2245, ω₂=-0.4818, ω₃= -0.6288 and ω₄=-0.6288, σ₁=0.025, σ₂=0.015；

Table 2

Table 2 is the experiment of different ligand molecular collected by matching MM-GBSA in present example with TcAChE Data, and by the scoring details of original MM-GBSA.PDB represents crystal structure database,^aRepresent from research The Medicine small molecule obtained in document and the dissociation constant of receptor,^bRepresent the combination arrived by obtaining dissociation constant conversion certainly By energy Δ G_binding-exp=RT ln K_D。

2) directly use the free energy computational methods improved to instruct docking sampling process as scoring functions, thus ensure to obtain The accuracy of Interactions Mode；

3) multiformity of horn of plenty joint mode, prevent tradition docking calculation in, due to combine energy add with unicity and wrong (shown in equation 1, the free energy that ligand molecular is combined with receptor protein is often by a series of energy for the multiformity of mistake joint mode Project superposition and obtain, therefore final energy is close, often but there are differences between each energetic terms, and this species diversity is past Toward being to cause due to the difference of Interactions Mode, traditional docking calculation is due to only using final combination energy as evaluating mark Standard, can be close thus tend to lose a lot of combination, and the Interactions Mode that structure there are differences), therefore, the present invention carries Go out each energy term respectively as object function, use the strategy of multiple-objection optimization docking to carry out joint mode sampling work, Thus avoid missing these Interactions Modes；

4) it addition, lose crucial local superiority in order to avoid the convergence of optimization method further and solve, the present invention is traditional In Multipurpose Optimal Method, incorporate elite solution archives preservation method, i.e. during optimizing, set up archives " container ", be used for protecting The all local superiorities obtained during depositing whole optimization solve, thus avoid further missing the advantage solution caused because of convergence.

Except solving the problems referred to above, the present invention, for improving docking speed, merges based on lattice point (Zou XQ, Sun in docking operation YX,Kuntz ID(1999)Inclusion of Solvation in Ligand Binding Free Energy Calculations Using the Generalized-Born Model.J Am Chem Soc 121:8033-8043.Meng EC,Shoichet BK,Kuntz ID(1992)Automated docking with grid-based energy evaluation.J Comput Chem 13:505-524.) With based on atom pair (Hawkins GD, Cramer CJ, Truhlar DG (1995) Pairwise solute descreening of solute charges from a dielectric medium.Chem Phys Lett 246:122-129.Wang JM,Cieplak P, Kollman PA(2000)How well does a restrained electrostatic potential(RESP)model perform in calculating conformational energies of organic and biological molecules？J Comput Chem 21: 1049-1074.) the method that both evaluates ligand molecular target receptor interaction force, for target receptor conformation mutability Part is defined as flexibility, and other parts are defined as rigidity, for flexible portion in docking operation with based on atom pair evaluation Method calculate the interaction force between itself and ligand molecular and the energy variation self caused because of conformation change；And it is firm Property part calculates it with ligand molecular with receptor flexible portion mutual with evaluation methodology based on lattice point in docking operation Active force, is defined as flexible residue the residue on prediction reaction channel, determines with number of residues and residue in this example Justice is set based on experience value by user.

In Interactions Mode sampling process, putting into a ligand molecular in each lattice, ligand molecular is at each lattice In be free to carry out change and the translation rotation of conformation, the calculating between subspace with subspace do not intersect with impact, Therefore this sampling process can be task number equal numbers of with subspace according to the number assignment of subspace, carries out parallel Calculate, the significantly Reduction Computation time.To in the docking sampling process of every sub spaces, enter initially with flexible docking method The Interactions Mode sampling that row is initial, for improving the multiformity of Interactions Mode further, simultaneously in order to avoid above-mentioned two Kind of (i.e. based on atom pair with based on lattice point) method merges and makes to calculate entirety and be combined energy generation error, then uses receptor complete The method being entirely considered as rigidity carries out the docking sampling of multiple target rigidity, carries out for institute's selectively lattice in each receptor conformation Sampling operation, this process uses scoring method based on lattice point, the Interactions Mode obtained similarly for flexible docking completely Also this lattice point marking form is used again to mark.So, the Conjugated free energy of all joint modes just unified by Method evaluation based on lattice point obtains, and while being enriched Interactions Mode, turn avoid existence between distinct methods Error and the inaccuracy that may cause.

Step S2.2: the Interactions Mode of the ligand molecular-target receptor of all subspaces is enriched with and obtains Interactions Mode Set P=∑ P_i=∑ ∑ p_k, as the present embodiment obtains all data such as Fig. 9, P=127,371.

Step S2.3: definition reaction coordinate axle, by each target receptor ligand molecular Interactions Mode according to reaction coordinate Axle projection, each Interactions Mode p_k(x y) then builds unit as target receptor ligand molecular Interactions Mode potential energy level Element, in the present embodiment, definition reaction coordinate is collectively formed by two variablees, i.e. the conformational difference of ligand molecular and part The instant configuration of molecule with it at the centroidal distance of the configuration of avtive spot.Certainly for building multidimensional potential energy level, reaction coordinate Number visually require and increase.

Step S2.4: by third party mapping software (such as matlab, R etc.), draws ligand molecular and target receptor phase interaction With the BFEL of pattern, such as Figure 10.Wherein x-axis and y-axis are reaction coordinate axle, and z-axis is Conjugated free energy coordinate axes.

After obtaining BFEL, then can calculate, according to it, the thermodynamics and dynamics ginseng that ligand molecular reacts with target receptor Number.

Step S3: obtain ligand molecular according to BFEL and occur to combine or dissociation reaction path with target receptor；Part is divided Son and the combination of target receptor or the exploration of dissociation reaction, two kinds of strategies of this research design, strategy (a): part is anti-with receptor Answer process not consider the ligand molecular backtracking effect in course of reaction, and follow strictly the principle of response path minimum energy.As Shown in Fig. 5, concrete step is as follows:

Step S3.1.1: calculate for convenience, first carries out coarseness division to BFEL,At this example Middle n=m=100, x=[0.0,0.04], y=[0.0,35.0], this division yardstick can be taken the circumstances into consideration according to the complexity of potential energy level etc. certainly Set.

Step S3.1.2: selected response path original position,Response path is risen Beginning position selected, in view of the potential energy level obtained in the present invention is simulation ligand molecular and naturally the interacting of target receptor Journey, there is no the intervention of any anthropic factor, and therefore the Dissociation path for obtaining can also be regarded as simultaneously and combines path.By In ligand molecular and target receptor identify the indefinability of initial position in reaction, the design explores dissociating of ligand molecular Path is main research approach, and for dissociation reaction, then reaction initial point should be ligand molecular with target receptor at avtive spot Bonding state, therefore for there being the system of crystal complex, then using the position of ligand molecular in crystal complex as initially Position, if not having crystal complex information, then regarding overall situation potential energy in whole potential energy level can initiate as dissociation reaction in minimum site Position.After having set original position, then initialisation reactions path S:=S ∪ S₀。

Step S3.1.3: ligand molecular is when reacting with target receptor, and ligand molecular will not occur excessive within the continuous print time Conformation or the change of configuration, simultaneously do not consider the backtracking effect in ligand molecular course of reaction, ligand molecular was dissociating Journey always faces 5 directions, then a bit that during its direction of advance is these five directions, potential energy is minimum, i.e. S_p=min{S_p(x_i-1,y_j), S_p(x_i-1,y_j+1), S_p(x_i,y_j+1),S_p(x_i+1,y_j+1),S_p(x_i+1,y_j)}；S:=S ∪ S_p+1,.By this side Method iteration is found until arriving reaction end.

Step S3.1.4: reaction end, is the iteration termination locations of route searching.In this example, with ligand molecular HupA Present position distance avtive spotPlace is as the terminal of dissociation reaction, and this numeral can be according to target receptor size, part Bulks of molecule etc. determine, the numerical value not used with this example, as standard, as long as ligand molecular dissociates is the most fully Can, it is judged that standard is away from target receptor, acts on completely with solution.

Step 3.15: iteration is obtained the coordinate information backtracking of response path S to the BFEL not having coarseness, it is thus achieved that phase With change functional relationship f (S) of the Conjugated free energy under coordinate Yu path, the result obtained such as this example, such as Figure 11 B.

The combination of ligand molecular and target receptor or the strategy (b) of dissociation reaction:

Complement one another with strategy (a), this programme not only consider the backtracking effect that ligand-receptor combines in dissociation reaction, On the basis of based on response path energy variation minimum theoretical, add the simulation operations of the random perturbation of small probability simultaneously, i.e. join Body molecule is in course of reaction, if if the different directions faced in the next step of reaction is likely to be of the feelings that energy is more or less the same Condition, then may continue reaction to randomly choose the form in direction.Idiographic flow such as accompanying drawing 5.The most important step is as follows:

Step S3.2.1: BFEL is carried out coarseness divisionN=m=100, with step S3.1.1；

Step S3.2.2: selected response path original position, i.e. potential energy level overall situation extreme lower position,Initialisation reactions path S:=S ∪ S₀, with step S3.1.2；

Step S3.2.3: judge whether to arrive defined reaction end, is then to enter step S3.2.12, otherwise proceeds to step S3.2.4

Step S3.2.4: according to neighbor node potential energy numerical value ascending sort S_pNeighbor node { S_p(x_i-1,y_j),S_p(x_i-1,y_j+1), S_p(x_i-1,y_j-1),S_p(x_i,y_j-1),S_p(x_i+1,y_j-1),S_p(x_i+1,y_j),S_p(x_i+1,y_j+1),S_p(x_i,y_j+1), return node ordered queue min_S_vec, Explore for subsequent path and prepare, enter step S3.2.5；

Step S3.2.5: check each present position S of part target_p(x_i,y_j) it is whether the site, border of potential energy level, if not Then proceed to step S3.2.6, if then proceeding to step S3.2.7；

Step S3.2.6: the present position judging part is which kind of the geomorphic feature site in potential energy level, these geomorphic features master It is divided three classes: potential well position, position, mountain peak and climbing position.If instant site is potential well position or position, mountain peak, Under greater probability (being set as 90% as in present example, can be any bigger numerical), ligand molecular all should be according to S3.2.4 The sequence node of middle return selects, and specifically chosen process proceeds to step S3.2.11.And remain under the probability of 10%, part The reaction direction of advance of molecule, for randomly choosing, i.e. proceeds to step S3.2.12.

Step S3.2.7: judge S_p(x_i,y_j) face potential energy value minimum position in node and the most lived through in track search process, Or had been subjected to but approach number of times is not more than 3 times that (this numeral can be regarded oneself for the tolerance of path multiplicity by user Depending on), if above-mentioned both of these case then proceeds to step S3.2.8, otherwise proceed to step S3.2.9.

Step S3.2.8: by present node S_pNeighbor node ordered queue min_S_vecIn the minimum node location of potential energy min_S_vec[j] is as next step moving direction (S:=S ∪ min_S_vec[j]), more new route present position (p=p+1), proceed to step Rapid S3.2.3；

Step S3.2.9: facing node queue min_S in order at random_vecOne direction of advance (min_S of middle selection_vec[j]), update Path present position (S:=S ∪ min_S_vec[j], p=p+1) and proceed to step S3.2.3；

Step S3.2.10: if S_p(x_i,y_j) be the climbing node in potential energy level, then with S_p(x_i,y_j) it is that boundary is by orderly for neighbor node team Row min_S_vecTwo points, free energy numerical value is more than S_p(x_i,y_j) node composition sequence be L_S_vec；Less than S_p(x_i,y_j) node group One-tenth sequence is S_S_vec, we determine next step moving direction of ligand molecular, random number according to the random number that program produces Use [0,1] random function to produce, for ensure the response path obtained be as far as possible in all response paths that may be present certainly Minimum by energy, in the case of the present invention is set in greater probability (when the random number i.e. obtained is in the range of less one, this Bright when being set as [0.3,1]) response path direction still according to minimum energy path direction, i.e. proceed to step S3.2.11, at S_S_vec Sequence selects；Otherwise, produce the effect of the random perturbation of response path, i.e. proceed to step S3.2.12, at L_S_vecIn carry out Select；

Step S3.2.11: at S_pThe ordered queue min_S of neighbor node_vecOr S_S_vecThe middle position min_S selecting minimum_vec[j] Or S_S_vec[j] as next step the Direction of Reaction of ligand molecular, with step S3.2.7, for avoid track search process because of hover in This position and be absorbed in endless loop, the number of times that this position occurs in response path is limited.Can be come regarding tolerance by user Defining, test case of the present invention is set to 8；If this position is unsatisfactory for the above-mentioned definition causing endless loop, proceed to step S3.2.12, Otherwise using this position as next step reaction site, Regeneration response path S (S:=S ∪ min_S_vec[j] or S:=S ∪ S_S_vec[j], P=p+1)；

Step S3.2.12: face node ordered queue min_S in current location_vecOr L_S_vecIn randomly choose a node min_S_vec[j] as next step the Direction of Reaction of path, more new route present position S (S:=S ∪ min_S_vec[j] or S:=S ∪ L_S_vec[j], p=p+1), and proceed to step S3.2.3；

Step S3.2.13: the ligand molecular that iteration is obtained and the coordinate information of the response path S of target receptor revert to into The BFEL that row coarseness divides, it is thus achieved that ligand molecular and target receptor response path and functional relationship f (S) of Conjugated free energy: That is, the response path S that will determine, according to its response path coordinate information, revert to the BFEL not having coarseness to process, Obtain accurate response path and obtain functional relationship: f (S), such as accompanying drawing 12 with Conjugated free energy.

Step S4: calculate corresponding thermodynamics, kinetic parameter according to the response path obtained；In this step, first root According to acquisition ligand molecular and the response path of target receptor, determine the binding site in course of reaction, transition state site and reaction The Conjugated free energy numerical value of termination site and these three critical sites, thus obtain and obtain combination freely according to transition state theory EnergyAssociation reaction activation energyWith dissociation reaction activation energyThus according to Aileen's formula (such as formula Acquisition thermodynamic parameter i.e. binding constant K 6-8) can be calculated_DWith kinetic parameter, i.e. association rate constant k_on, speed of dissociating Rate constant k_off。

$K_D=\frac{k_{off}}{k_{on}}---\left(6\right)$

ΔG_binding=RT ln K_D (7)

${\mathrm{\ΔG}}_{on}^{\≠}=-RT\ln \left(\frac{k_{on}h}{k_{B}T}\right),{\mathrm{\ΔG}}_{off}^{\≠}=-RT\ln \left(\frac{k_{off}h}{k_{B}T}\right)---\left(8\right)$

K_BFor Boltzmann constant, R is gas constant, and T is absolute temperature (this example uses 298.15K), and h is Pu Lang Gram constant.

Step S5: obtain the mechanism of action of molecule reaction, verify the key binding sites determining reaction kinetic characteristics, as Medicine (part) design key reference factor.By obtaining response path, it may be determined that ligand molecular is in whole response path The configuration of each present position, revert to dock the Interactions Mode obtained and i.e. can obtain the structural information of course of reaction Details.Thus, it is possible to acquisition ligand molecular is in active binding site, transition state site structure information, learn on target receptor Avtive spot and transition state site, such as accompanying drawing 13 and Figure 14.

The acquisition of kinetic parameter is of great significance, because can obtain part (medicine) point by dissociation rate constant The holdup time of son, i.e. lasting medicine effect, the pharmacological effect for research part (medicine) molecule serves the most crucial Meaning.

In the present embodiment, with TcAChE (target receptor)-HupA (ligand molecular) as example, use proposed by the invention Example system is simulated calculating by method, it is thus achieved that HupA with TcAChE thermodynamics and dynamics parameter when reacting, Simulation calculates the numerical value obtained and is sufficiently close to the numerical value that experiment test obtains.Wherein, Conjugated free energyOnly with experiment Value differs about 1kJ/mol, and kinetic parameter, i.e.WithCompare with experiment value, the poorest 1.32kJ/mol and 3.61kJ/mol.This shows that method of the present invention can calculate the knot of prediction ligand molecular and target receptor quickly and accurately Close affinity and kinetic parameter, and then be expected to be generalized to medicine (part) design field and instruct excavation or design improvement leading Compound.

Claims (16)

1. simulate the side that ligand molecular reacts with target receptor and calculates the thermodynamics and dynamics parameter predicting this reaction for one kind Method, said method comprising the steps of:

Step S1: download the described target receptor obtained by laboratory facilities from online target receptor crystal structure database With the three-dimensional structure data of ligand molecular, three dimensional structure based on ligand molecular and target receptor is built ligand molecular and is subject to target The simulated system of body；

Step S2: build the BFEL of ligand molecular and target receptor；

Step S3: obtain ligand molecular according to BFEL, with target receptor, combination or the response path dissociated occur；

Step S4: calculate corresponding thermodynamics and kinetics parameter according to the response path obtained；

Step S5: obtain the mechanism of action of molecule reaction by response path, determine the pass bond of decision reaction kinetic characteristics Close site,

Wherein, described BFEL refers to that Conjugated free energy potential energy level, described step S2 comprise the following steps:

Step S2.1: in every sub spaces l_i∈ L puts into a ligand molecular, carries out docking of ligand molecular and target receptor Calculate, every sub spaces l_iObtain a target receptor-ligand molecular Interactions Mode set P_i=∑ p_k, wherein, k is big In the positive integer equal to 0, p_kRepresent the kth target receptor-part Interactions Mode obtained；

Step S2.2: collect the interaction of calculated target receptor-ligand molecular in all subspaces in step S2.1 Pattern obtains P=∑ P_i, it is characterised in that the ligand molecular of all subspaces-target receptor Interactions Mode enrichment is obtained To Interactions Mode set P=∑ P_i=∑ ∑ p_k, as the data element building BFEL；

Step S2.3: the definition coordinate axes that reacts with target receptor of ligand molecular, mutual by each ligand molecular-target receptor Binding mode p_k(x, y) carries out projection to the reaction coordinate axle reacted at ligand molecular-target receptor, each interaction mould Formula p_k(x y) combines potential energy level as target receptor-ligand molecular and builds element；

Step S2.4: by the three-dimensional of third party's three-dimensional drawing Software on Drawing ligand molecular Yu target receptor Interactions Mode BFEL, wherein, x-axis and y-axis are reaction coordinate axle, and z-axis is Conjugated free energy coordinate axes,

Wherein, described step S3 farther include strategy (a) or strategy (b):

Strategy (a) comprises the following steps:

Step S3.1.1: BFEL is carried out coarseness divisionΩ represents whole BFEL,Respectively Representing the lattice point number of the upper division of reaction coordinate axle x, y of whole potential energy level, n, m are respectively the lattice point size of space, determine anti- Answer coordinate axes x, y upper lattice point number；

Step S3.1.2: selected response path original positionx_iRepresent in reaction coordinate x I-th lattice point, y_jRepresenting the jth lattice point in reaction coordinate y, initialisation reactions path is S:=S ∪ S₀；

Step S3.1.3:S_p=min{S_p(x_i-1,y_j), S_p(x_i-1,y_j+1), S_p(x_i,y_j+1),S_p(x_i+1,y_j+1),S_p(x_i+1,y_j)}；S:=S ∪ S_p+1, On the BFEL of coarseness, it is iterated search according to minimum energy principle and obtains the reaction road of ligand molecular and target receptor Footpath, wherein each present position of ligand molecular all faces 5 advanceable directions, next step moving direction of ligand molecular Being the minimum direction of Conjugated free energy in 5 directions, wherein S is for finally to obtain path, S_pRepresent ligand molecular to be subject to target The pth step of body node location in the potential energy level that coarseness divides, p is the positive integer more than or equal to 0；

Step S3.1.4: judge described S_PWhether arrive user-defined reaction end, be then to enter step S3.1.5, otherwise Proceeding to step S3.1.3, wherein, when studying the association reaction of ligand molecular and target receptor, reaction end is target receptor On catalysis or the avtive spot of binding partner molecule；When studying the dissociation reaction of ligand molecular and target receptor, reaction is eventually Point is the region away from target receptor, and the described region away from target receptor is solvent area；

Step S3.1.5: ligand molecular iteration obtained and the coordinate information backtracking of target receptor response path S are to not carrying out The BFEL that coarseness divides, it is thus achieved that ligand molecular and target receptor response path and functional relationship f (S) of Conjugated free energy, root According to response path coordinate, it revert to totipotency face, thus Regeneration response path corresponding potential energy information: f (s), wherein, described S Coordinate information be corresponding x, the position of y-axis, described totipotency face is the BFEL not carrying out coarseness division；

Strategy (b) comprises the following steps:

Step S3.2.1: BFEL is carried out coarseness divisionWith step S3.1.1；

Step S3.2.2: with step S3.1.2, selected response path original position,Just Beginningization response path S:=S ∪ S₀；

Step S3.2.3: judge present position S_PWhether arrive defined reaction end, be then to enter step S3.2.12, no Then proceed to step S3.2.4；

Step S3.2.4: according to neighbor node potential energy numerical value ascending sort S_pNeighbor node { S_p(x_i-1,y_j),S_p(x_i-1,y_j+1), S_p(x_i-1,y_j-1),S_p(x_i,y_j-1),S_p(x_i+1,y_j-1),S_p(x_i+1,y_j),S_p(x_i+1,y_j+1),S_p(x_i,y_j+1), return node ordered queue min_S_vec, prepare for subsequent reactions track search, enter step S3.2.5；

Step S3.2.5: check each present position S of part target_p(x_i,y_i) it is whether the site, border of potential energy level, if not Then proceed to step S3.2.6, if then proceeding to step S3.2.7；

Step S3.2.6: the present position judging part is which kind of the geomorphic feature site in potential energy level, when instant site is gesture When trap position or position, mountain peak, the random number produced according to program determines next step moving direction of ligand molecular, random number [0,1] random function is used to produce, if in the case of value is in the big probability of [0,0.9], ligand molecular is in step S3.2.4 Ordered nodes sequence min_S returned_vecSelect successively, proceed to following steps S3.2.11；When random number value be in [0.9, 1.0], in the case of small probability, the reaction direction of advance of ligand molecular, for randomly choosing, i.e. proceeds to following steps S3.2.12；Its In, described geomorphic feature site includes potential well position, position, mountain peak and climbing position；

Step S3.2.7: judge S_p(x_i,y_i) adjacent node in potential energy value minimum position whether already contained in ligand molecular and target In mark receptor response path S, if this position is not in response path S, then proceed to step S3.2.8；If this position is It is included in and explores in the response path obtained, and this position occurs once, then for avoiding path in response path incessantly Heuristic process is absorbed in endless loop because hovering in this position, is limited by the number of times that this position occurs in response path；Its In, when some position occurs M time in response path, then it is assumed that this position can cause endless loop, M is more than 0 Positive integer, M is defined depending on tolerance by user, is closer to 1, the most whole explore obtain response path in redundancy function Put the fewest, reflect ligand molecular and " roundabout " feature of generation in target receptor cohesive process more than 1；If should Position occurrence number in response path less than M, then proceeds to following steps S3.2.8；If this position occurs in response path Number of times is more than or equal to M, then proceed to following steps S3.2.9；

Step S3.2.8: by present node S_pNeighbor node ordered queue min_S_vecIn the minimum node location of potential energy min_S_vec[j] is as next step moving direction S:=S ∪ min_S_vec[j], more new route present position p=p+1, proceeds to step S3.2.3；

Step S3.2.9: random at neighbor node ordered queue min_S_vecOne direction of advance min_S of middle selection_vec[j], more New route present position S:=S ∪ min_S_vec[j], p=p+1 also proceeds to step S3.2.3；

Step S3.2.10: if S_p(x_i,y_j) be the climbing node in potential energy level, then with S_p(x_i,y_j) it is that boundary is by orderly for neighbor node team Row min_S_vecTwo points, free energy numerical value is more than S_p(x_i,y_j) node composition sequence be L_S_vec；Less than S_p(x_i,y_j) node group One-tenth sequence is S_S_vec, wherein, the random number produced according to program determines next step moving direction of ligand molecular, at random Number use [0,1] random functions produce, be set in the random number obtained in the range of [0.3,1] time, response path direction is still pressed According to minimum energy path direction, i.e. proceed to step S3.2.11, at S_S_vecSequence selects；Otherwise, produce response path with The effect of machine perturbation, i.e. proceeds to step S3.2.12, at L_S_vecIn select；

Step S3.2.11: at S_pThe ordered queue min_S of neighbor node_vecOr S_S_vecThe middle position min_S selecting minimum_vec[j] Or S_S_vec[j] as next step the Direction of Reaction of ligand molecular, with step S3.2.7；Wherein, for avoiding track search process In, ligand molecular is absorbed in endless loop because hovering in this position, this position occurs according to user's tolerance in response path Number of times limit；If this position is unsatisfactory for the above-mentioned definition causing endless loop, proceed to step S3.2.12, otherwise by this position Put the reaction site as next step, more new route present position S_p, wherein S:=S ∪ min_S_vec[j] or S:=S ∪ S_S_vec[j], P=p+1；

Step S3.2.12: at the neighbor node ordered queue min_S of current location_vecOr L_S_vecIn randomly choose a node L_S_vec[j] is as next step the Direction of Reaction of path, more new route present position S_p, wherein S:=S ∪ min_S_vec[j] or S:=S ∪L_S_vec[j], p=p+1, and proceed to step S3.2.3；

Step S3.2.13: the ligand molecular that iteration is obtained and the coordinate information of the response path S of target receptor revert to into The BFEL that row coarseness divides, it is thus achieved that ligand molecular and target receptor response path and functional relationship f (S) of Conjugated free energy.

Method the most according to claim 1, wherein,

In step S3.2.7, M is 3.

Method the most according to claim 1, wherein,

In step S3.2.11, the number of times of described appearance is set to 8.

Method the most according to claim 1, wherein, described online target receptor crystal structure database is RCSB, Described laboratory facilities are X-ray diffraction and nuclear-magnetism diffraction, and described ligand molecular is Medicine small molecule or compound, described target Receptor is disease-associated protein or nucleic acid.

Method the most according to claim 1, wherein, described online target receptor crystal structure database is RCSB, Described laboratory facilities are X-ray diffraction and nuclear-magnetism diffraction, and described ligand molecular is huperzine A, and described target receptor is disease Related protein or nucleic acid.

Method the most according to claim 1, wherein, described online target receptor crystal structure database is RCSB, Described laboratory facilities are X-ray diffraction and nuclear-magnetism diffraction, and described ligand molecular is Medicine small molecule or compound, described target Receptor is acetylcholinesterase.

Method the most according to claim 1, wherein, described online target receptor crystal structure database is RCSB, Described laboratory facilities are X-ray diffraction and nuclear-magnetism diffraction, and described ligand molecular is huperzine A, and described target receptor is acetyl Acetylcholine esterase.

Method the most according to claim 1, wherein, described step S1 comprises the following steps:

Step S1.1: the size range in definition research system space, it is characterised in that come by the three dimensional structure of target receptor Delimit and intend the bulk scope that sampling ligand molecular interacts with target receptor；

The described bulk scope three-dimensional more than described target receptor intending sampling ligand molecular and target receptor interaction 10～20 angstroms of structure；

Step S1.2: the described research system space of step S1.1 definition is carried out subspace stress and strain model, and empty to every height Between be marked, with l_iCarry out labelling, represent i-th subspace；

Step S1.3: the crucial subspace of screening joint mode sampling；

Wherein, described step S1.3 comprises the following steps:

Step S1.3.1: download the described target obtained by laboratory facilities from online target receptor crystal structure database The crystal structural data of receptor, described online target receptor crystal structure database is RCSB, and described laboratory facilities are X-ray diffraction and nuclear-magnetism diffraction, the apo crystal knot that described crystal structure only target receptor monomer structure coexists without little molecule Structure and the holo structure of target receptor complex coexisted with little molecule；

Step S1.3.2: use combine for the little molecule of part on target receptor or dissociation channel calculate sampling MOLE or CAVER software, calculates the passage t that each target receptor crystal structure is potential_i, and it is enriched with all passage T=∑ t_i, its feature It is, uses described software, passage t potential on each target receptor of sampling successively_i, obtain on each target receptor the most at last The passage obtained is integrated, i.e. T=∑ t_i, wherein, i represents the i-th target calculated, and i depends on the required meter of user Calculate the three dimensional structure number of target receptor, for the positive integer more than 1；t_iFor potential in the i-th structure that obtained by calculating Passage；

Step S1.3.3: the described passage T clustering all keys obtains representative passage T_r, analyze S1.3.2 obtains and own Passage T, the passage of each orientation only retains one, and other filter for redundancy；

Step S1.3.4: select all coverings or neighbouring T_rSubspace l_iAs the crucial subspace of joint mode sampling, it is special Levy and be, according to described passage T_rFrom step S.1.2 all subspaces of obtaining filter out crucial subspace Σ l_i；

Step S1.3.5: described crucial subspace is carried out BORDER PROCESSING, it is characterised in that interleaving of adjacent subspace Enter to cover the extra subspace l on border_j, carry out labelling with j, insert Σ l altogether_jIndividual extra subspace covers border；

Step S1.3.6: the final crucial subspace L=Σ l obtaining joint mode sampling_i∪Σl_j。

Method the most according to claim 8, wherein,

In step S1.3.2, use the most possible of the desired acquisition of the self-defined each calculating of MOLE or CAVER software Number of active lanes, i is defaulted as 3.

Method the most according to claim 1, wherein,

Described step S2.1 is characterised by, simulation ligand molecular in described subspace with all phase interactions of target receptor With using the molecular docking method developed based on multi-objective optimization algorithm, described multiple target docking calculation during pattern it is Flexible docking method, in described flexible docking method, includes following two scheme for the structure of object function:

By Δ E_vdw,ΔE_es,ΔG_gb,ΔG_sa... as multiple targets or the fitness function of docking；Or

By Δ E_intra,ΔE_inter,ΔG_gb,ΔG_sa... as multiple targets or the fitness function of docking,

Wherein, Δ E_vdwFor van der Waals energy quantifier, Δ E_esFor electrostatic energy quantifier, Δ G_gbFor solvent polarization free energy, Δ G_saFor The non-polarized free energy of solvent；ΔE_intraFor target receptor, the energy term sum of ligand molecular self, including: ΔE_interFor the interaction force sum between target receptor ligand molecular, including: Δ E_vdw,RLWith ΔE_es,RL, Δ G_gb, Δ G_sa,

Wherein,It is to be changed when identifying by the conformation of target receptor with ligand molecular and cause respectively The van der Waals energy quantifier of target receptor and the variable quantity of electrostatic term；It is by target receptor and part respectively The van der Waals energy quantifier of the ligand molecular that the conformation of molecule changes when identifying and causes and the variable quantity of electrostatic term； ΔE_vdw,RLConformation for target receptor with ligand molecular changes when identifying and causes the phase of target receptor ligand molecular The variable quantity of interaction Van der Waals energy term, Δ E_es,RLConformation for target receptor and ligand molecular change when identifying and Cause the variable quantity of the interaction electrostatic term of target receptor ligand molecular；

Target numbers employed in described multi-objective optimization algorithm is the elastic system that sets, and i.e. selects between 1-n, and n is user The object function number of definition, n is the positive integer more than 1.

11. methods according to claim 1, wherein,

Described step S2.1 is characterised by, simulation ligand molecular in described subspace with all phase interactions of target receptor During pattern, using rigidity to process with flexible region division target receptor, wherein, flexible region considers that target is subject to Body portion structural variability, and use energy scoring tactics based on mesh point Yu atom pair to process described flexible region and just Property region.

12. methods according to claim 11, wherein,

Described target receptor part-structure transmutability is the motion change of protein Key residues or domain.

13. methods according to claim 12, wherein,

Described domain is ring region.

14. methods according to claim 1, wherein,

Described step S4 is characterised by, according to the response path of the ligand molecular obtained with target receptor, determines and reacted Avtive spot f (S) in journey_{Minimum position}, transition state site f (S)_{Maximum position}With reaction terminating site, wherein,

According to f (S)_{Minimum position}, transition state site and the Conjugated free energy data △ G in reaction terminating site_bind,△G_trans, △ G_unbind, in conjunction with transition state theory acquisition Conjugated free energy:Association reaction activation energy:With dissociation reaction activation energy:And then according to Aileen's formula, Equation below (6)-(8) calculate and obtain binding constant K_D, association rate constant k_onWith rate constants k of dissociating_off；Obtain further Obtain the action time of ligand molecular and target receptor: 1/k_off；

$K_D=\frac{k_{off}}{k_{on}}---\left(6\right)$

${\mathrm{\ΔG}}_{binding}^o=RT\ln K_D---\left(7\right)$

${\mathrm{\ΔG}}_{on}^{\≠}=-RTln\left(\frac{k_{on}h}{k_{B}T}\right),{\mathrm{\ΔG}}_{off}^{\≠}=-RTln\left(\frac{k_{off}h}{k_{B}T}\right)---\left(8\right)$

Wherein, K_BFor Boltzmann constant, R is gas constant, and T is absolute temperature, and h is planck constant.

15. methods according to claim 1, wherein,

Described step S5 is characterised by, revert to free energy potential energy level by the response path obtained and obtains and join on this path Each present position S of body molecule_p(x_p,y_p) reaction coordinate (x_p,y_p), found by reaction coordinate and be projected in this position Ligand molecular and the binding pattern P of target receptor_i(x_i,y_i), at this x_p=x_i, y_p=y_i, it is derived from ligand molecular and ties in activity Close site, the binding pattern with target receptor in transition state site, learn the avtive spot on corresponding target receptor and mistake Cross state site.

16. 1 kinds for simulating the thermodynamics and dynamics parameter that ligand molecular reacts with target receptor and calculates this reaction of prediction Computer based application system, including:

First module, it performs step S1 as claimed in claim 1；

Second module, it performs step S2 as claimed in claim 1；

Three module, it performs step S3 as claimed in claim 1；

4th module, it performs step S4 as claimed in claim 1；

5th module, it performs step S5 as claimed in claim 1.