CN102930152A - Method and system for simulating ligand molecule and target receptor reaction and calculating and forecasting thermodynamics and kinetics parameters of reaction - Google Patents

Abstract

The invention relates to a method and application system for simulating ligand molecule and target receptor reaction based on a computer and calculating and forecasting thermodynamics and kinetics parameters of the reaction. The method mainly comprises the following steps of: designing a method for calculating combined free energy of a ligand molecule and a target receptor, and designing a multi-target molecule docking method and developing an application system of the docking method based on the calculating method; simulating and sampling the interaction mode of the ligand molecule and the target receptor by the system, grading the free energy combination of the ligand molecule and the target receptor, and building a method and system for the potential energy surface of the combined free energy of the ligand molecule and the target receptor according to the simulating and sampling result; building a method and system for identifying a reaction path of the ligand molecule and the target receptor according to the potential energy surface; and determining a active site and a transition site of the reaction to calculate the thermodynamics and kinetics parameters. The method and the system provide important basis and a forecast evaluation method for drug development, and are expected to improve the drug development efficiency and save experimental expenses greatly, thus serving as an effective auxiliary means for the drug development.

Description

A kind of method and system of simulating ligand molecular and target receptor reaction and calculating the thermodynamics and dynamics parameter of this reaction of prediction

Technical field

The invention belongs to medicine, chemistry, calculation biology and drug design field, relate to a kind of based on the computer simulation ligand molecular (for example, drug molecule) with target receptor (such as albumen, the disease targets such as nucleic acid) reaction, the new method of the thermodynamics and dynamics parameter of this reaction is obtained in calculating.Described method can accurately simulate ligand molecular and target receptor the mode of action of molecular level and reaction (in conjunction with/dissociate) mode, and calculate quickly and accurately and obtain the thermodynamics and dynamics parameter of this reaction, obtain combination or the Dissociation path of the reaction of ligand molecular and target receptor, and then excavate the key effect site (avtive spot and transition state site) of target receptor-ligand molecular, be applicable to quantitatively, Interactions Mode and effect between qualitative analysis ligand molecular and target receptor be strong and weak.Provide foundation for designing or transform more highly active drug molecule, characterize for the screening new compound provides quantitative thermodynamics and dynamics, improve the success ratio of discovery reactive compound.

Background technology

Molecular recognition has vital effect in biosystem.The Interactions Mode of ligand molecular and target receptor can cause a series of physiology and pharmacological effect (CharlesR.Cantor PRS (1980) Biophysical Chemistry:PartIII: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 Biol91:357-364.), this identifying comprises combination or the dissociation reaction of ligand molecular and target receptor.Especially because the drug effect of most medicines mainly depends on the direct cohesive process of part and corresponding target, so drug design then with part and target receptor binding isotherm as basic foundation.Given this, the binding mechanism of understanding in detail part and acceptor is vital (Taft CA to drug design then, DA Silva VB, Da Silva CH (2008) Current topics in computer-aided drugdesign.J Pharm Sci 97:1089-1098.).Just pay close attention to the size of the binding affinity of ligand molecular and target receptor in traditional drug design process, generally believe that the strong molecule of affinity can stronger than the Molecule Action a little less than the affinity (Zhang R, Monsma F (2010) Binding kinetics and mechanism of action:toward the discovery anddevelopment of better and best in class drugs.Expert Opin Drug Discov 5:1023-1029.).But this simple extremely low based on the drug design method success ratio of affinity.Tracing it to its cause, is because this method for designing has been ignored the kinetic property of reaction ligand molecular and target receptor cohesive process.In fact, drug effect not only with thermodynamic parameter (as: in conjunction with free energy (Δ G _Binding), binding constant (K _A) and dissociation constant (K _D)) relevant, also with binding kinetics parameter (as: association rate constant (k _On), kinetic constant (k _Off), the hold-up time (1/koff) of drug molecule and the combination of target receptor) relevant (KenakinT, Jenkinson S, Watson C (2006) Determining the potency and molecular mechanism of action ofinsurmountable antagonists.J Pharmacol Exp Ther 319:710-723), the hold-up time can be thought the onset duration of drug molecule.The binding affinity of a lot of examples prove drug effects and drug molecule and target receptor there is no direct correlativity, and present stronger correlativity (Swinney DC (2004) Biochemical mechanisms of drug action:what does it take for success Nat Rev Drug Discov 3:801-808.) with hold-up time of the combination of drug molecule and target receptor, therefore, drug design method based on affinity is not comprehensive, the thermodynamics and dynamics method of paying attention to simultaneously the mechanism of drug molecule and acceptor has then obtained paying attention to widely (Zhang R, Monsma F (2010) Bindingkinetics and mechanism of action:toward the discovery and development of better and best inclass drugs.Expert Opin Drug Discov 5:1023-1029.).

Although the kinetic property that medicine (part) molecule is combined with target receptor is out revealed in the importance of medicine (part) design field, can be fast but also do not develop in the last thirty years,, the method for Accurate Prediction part and receptors bind kinetic property.In addition, although in the research in the past, emphasize the calculating Forecasting Methodology of affinity (being macroscopic property), also there is certain defective in existing thermodynamic calculation method always.Although developed such as the higher method of this accuracy of computation of free energy perturbation (FEP), because its computing time is oversize and so that there is great limitation in its application.And on the contrary, the scoring functions that adopts in the molecular docking process then has relatively fast advantage, but this method is owing to means such as a large amount of employings is approximate, simplification, so that accuracy in computation reduces greatly.MM-GBSA/PBSA between these two kinds of methods has more accurately and faster advantage (Wang, J.M.; Hou, T.J.; Xu, X.J. (2006) Recent advances in free energy calculations with a combination of molecular mechanics andcontinuummodels.Curr Comput-Aided Drug Des 2:287 – 306.).The method was widely used in the last few years, and this method is proved often has the system dependence, though can predict accurately for some system, can not guarantee all exactly predictions of all systems.

In the last few years, the researcher has also proposed a kind of method based on molecular dynamics and has made up in conjunction with potential energy surface or calculate the scheme of thermodynamic and kinetic properties, this scheme is at first by adopting up to a hundred even thousands of the kinetic locus of molecular dynamics simulation sampling, then by these tracks, obtain the ligand receptor kinetic parameter according to the corresponding theory computing method.Therefore this method need to expend a large amount of computational resources and time owing to need by a large amount of sample track of dynamic method, also being inaccurate of the thermodynamic parameter that simultaneously prediction obtains and kinetic parameter.Make up zamidine and tryptic interaction potential energy face such as Buch et.al. by 495 kinetic locus of sampling, thereby predicted more accurately macroscopic property (with test phase ratio 1kcal/mol) according to potential energy surface, but comparing with test, the kinetic property that calculates then has deviation (the Buch I that has an order of magnitude at least, Giorgino T, De Fabritiis G (2011) Complete reconstruction of an enzyme-inhibitorbinding process by molecular dynamics simulations.Pro Natl Acad Sci USA 108:10184-10189.).

Therefore, developing the thermodynamic parameter that can predict accurately and rapidly the reaction of ligand molecular-target receptor and the method for kinetic parameter still is problem demanding prompt solution.

Summary of the invention

Technical matters

The technical problem to be solved in the present invention is: the interaction mode of simulation ligand molecular and target receptor, obtain ligand molecular and target receptor (in conjunction with or dissociate) reaction mode, calculate thermodynamic parameter and the kinetic parameter of this reaction, the crucial binding site of experience when determining the reaction of ligand molecular and target receptor, and then instruct design or the transformation of drug molecule, overcome simultaneously simulation method length consuming time and inaccurate defective in the prior art, improve medicament research and development efficient, reduce R﹠D costs.

Particularly, purpose of the present invention is: according to ligand molecular and the target receptor interactional natural law under field conditions (factors), simulation and the possible Interactions Mode of ligand molecular and target receptor of systematically sampling, quantize the acting force (namely in conjunction with free energy) of each binding mode, accurately, make up all sidedly ligand molecular and target receptor Interactions Mode in conjunction with free energy potential energy surface (Binding Free Energy Landscape, BFEL), and according to potential energy surface come fast, calculate exactly that the prediction ligand molecular is combined with target receptor or dissociation reaction process and this course of reaction in thermodynamic parameter, kinetic parameter; Obtain on the target receptor determine ligand molecular in conjunction with or the path that dissociates on crucial binding site, for the design of ligand molecular, particularly drug molecule or transformation provide new approaches and effective foundation.

Technical scheme and beneficial effect

Therefore, the invention provides the mode of action of accurate simulation and exhaustive ligand molecular and target receptor, explore the reaction (combination and the reaction of dissociating) of ligand molecular target receptor, and calculate new method and the application system thereof of the thermodynamics and kinetics parameter of obtaining this reaction.For reaching this purpose, the present invention is according to the natural law of ligand molecular and target receptor interaction pattern, propose improved simulation and assessed the molecular docking method of this Interactions Mode mode and the method for MM-GBSA/PBSA, and merge complete ligand molecular and the interactional BFEL of target receptor of structure that these two kinds of methods are come system, and then be combined with target receptor as the prediction ligand molecular or the basis in dissociation reaction path with the BFEL that obtains, adopt the response path method of sampling to obtain combination or dissociation reaction path, thereby calculate corresponding thermodynamic parameter and kinetic parameter.

Particularly, method of the present invention is for adopting computer program to realize the scheme that external data is processed, external data is to possess skills the technical data of implication (such as part (medicine) molecule of processing among the present invention, and target receptor molecule), processing procedure is according to the natural law of intermolecular interaction pattern, described method can realize fast, simulate exactly and estimate drug molecule and the disease targets mode of action, thereby calculate the thermodynamic parameter of intermolecular interaction pattern and kinetic parameter and then instruct SARS drug design or transformation and reach beneficial effect.Therefore, method of the present invention is for adopting the technical scheme of technological means, realization technique effect.

Described method is calculated thermodynamic parameter and the kinetic parameter of prediction ligand molecular when target receptor is combined or dissociate by following five key steps (as shown in Figure 1).Wherein, described thermodynamic parameter is: in conjunction with free energy (Δ G _Binding), binding constant (K _A), dissociation constant (K _D); Described kinetic parameter is: combination and the free energy of activation that dissociates

Association rate constant and the rate constant (k that dissociates _On, k _Off):

Step S1: from online target receptor crystal structure database, as RCSB ( Http:// www.pdb.org/pdb/home/home.do) three-dimensional structure data of the target receptor of studying that obtains of the means (such as X-ray diffraction and nuclear-magnetism diffraction) by experiment downloaded, measure the space structure data of ligand molecular and target receptor, build the simulated system of ligand molecular and target receptor based on the space structure of ligand molecular and target receptor;

Step S2: the BFEL that makes up ligand molecular and target receptor;

Step S3: obtain the response path that ligand molecular and target receptor combination occur or dissociate according to BFEL;

Step S4: calculate corresponding thermodynamics and kinetics parameter according to the response path that obtains;

Step S5: obtain the mechanism of action of molecular reaction, definite crucial binding site that determines the reaction kinetics feature by response path.

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

Wherein, described step S1 is: build ligand molecular and target receptor simulated system, the ligand molecular of namely plan being studied, the pre-treatment step of target receptor and definition simulated system spatial dimension size etc.Described step S1 may further comprise the steps (as shown in Figure 2):

Step S1.1: definition research system bulk scope.Be specially: the three-dimensional structure of testing the target receptor that obtains by assessment delimited the bulk scope of intending sampling ligand molecular and target receptor Interactions Mode, here be pointed out that, for the division of this range of size with define, mainly depend on the three-dimensional structure size of target receptor, simultaneously can obtain binding mode and the acting force of ligand molecular in solution in order to guarantee, therefore this spatial dimension should be fully large, thus can provide enough solution space for ligand molecular in the situation that the acting force that is not subject to target receptor disturbs with solution reaction.In addition, if crossing conference, this bulk causes producing a large amount of meaningless calculating.In view of the foregoing, the present invention advises that this bulk should be than large 10 ~ 20 dusts of target crystal three-dimensional structure, certainly this numeral there is no harsh defining, other size, all should be protected such as 20 ~ ∞, just too large, be unalterable solution zone such as the space that surpasses 20 dusts, there is no and calculate value;

Step S1.2: the subspace grid is carried out in the research space of delimiting divide.Be specially: the subspace grid is carried out in the described research system space of step S1.1 definition divide, and every sub spaces is carried out mark, the present invention is with l _iComing mark, represent the i sub spaces, is several little spaces with the whole research spatial division of delimiting, and is called the subspace.This division purpose has two: 1) can be with whole analog computation Process Decomposition, both put into ligand molecular in every sub spaces, the sample binding mode of ligand molecular and target receptor in this subrange, this mode can so that task decompose, be convenient to parallel computation, thereby increase substantially counting yield; 2) avoided because optimization method convergence causes missing important ligand molecular and the Interactions Mode of target receptor, such as the ligand molecular of the transition state of decision kinetic property and the Interactions Mode of target receptor.Because the optimization method meeting of adopting in the docking sampling process is so that the Interactions Mode that sampling obtains converges to global optimum's (minimum in conjunction with free energy) or local optimum (lower in conjunction with free energy), if do not carry out the partial restriction sampling, the Interactions Mode that then obtains will inevitably converge on the avtive spot place (this place's part and receptors bind free energy are theoretical for overall minimum) of target receptor substantially, and the Interactions Mode in conjunction with the higher transition state place of free energy that determines reaction kinetics character can be missed, thereby so that can't calculate the reaction power mathematic(al) parameter that obtains ligand molecular and target receptor;

The beneficial effect of described step S1.2 is: the subspace partition strategy is so that computation process can parallelization, increase substantially the time, also can effectively avoid simultaneously because the optimization method convergence causes missing important Interactions Mode, as determining the transition state Interactions Mode of kinetic property.

Step S1.3: screen crucial subspace.Be specially the space that obtains among the step S1.2 is screened, will pick out with the larger subspace of ligand molecular effect possibility in this space, and the less subspace of possibility is got rid of; The beneficial effect of described step S1.3 is: avoid insignificant calculating, the waste computational resource.

Described step S1.3 may further comprise the steps:

Step S1.3.1: from online target receptor crystal structure database, as RCSB ( Http:// www.pdb.org/pdb/home/home.do) three-dimensional structure data of the target receptor of studying that obtains of the means (such as X-ray diffraction and nuclear-magnetism diffraction) by experiment downloaded, described crystal structure comprises apo(namely only for the target receptor monomer structure, without little molecule coexistence) and holo structure (the target receptor composite structure that namely coexists with little molecule).Because the target receptor structure has behavioral characteristics, therefore under the varying environment condition, or in conjunction with in the different micromolecular situations, can there be certain difference in its three-dimensional conformation, then corresponding channel configurations from the ligand molecular reaction also may be different, therefore by all known crystal structures of enrichment, the skewed popularity on the result of calculation that can avoid being brought by the three-dimensional structure of target receptor;

Step S1.3.2: adopt the third party software that supplies the little molecule combination of part or dissociation channel calculating sampling on the target receptor, calculate the potential passage t of each target receptor crystal structure such as MOLE or CAVER etc. _i, and all passage T=∑ t of enrichment _i, it is characterized in that, adopt third party software, potential passage t on each target receptor of sampling successively _i, the passage that obtains on each target receptor is the most at last integrated T=∑ t _i, (here, i the target receptor three-dimensional structure that i representative is calculated, i depends on the three-dimensional structure number of the required calculating target receptor of user, for greater than 1 positive integer; t _iBe potential passage on i the structure that obtains by calculating; The MOLE that the present invention enumerates or CAVER can the desired most possible number of active lanes of obtaining of self-defined each calculating, are defaulted as 3.）。

Step S1.3.3: all critical passage T of cluster obtain representative passage T _r, (concrete: analyze and obtain all passage T among the S1.3.2, the passage of each orientation only keeps one, and other filter for redundant.）；

Step S1.3.4: select all coverings or contiguous T _rSubspace l _iAs the crucial subspace of joint mode sampling study portion.Be specially: according to these passages T _rCome from step S.1.1 with S.1.2 select the space that implement as next step Interactions Mode sampling operation crucial subspace in all subspaces of obtaining.Certainly can select institute to have living space, and system of selection also can suitably adjust, as comprise other the outer subspaces, potential binding site place of passage on the target receptor, these all should be included among this patent scope;

Step S1.3.5: boundary treatment is carried out in described crucial subspace, it is characterized in that, between adjacent subspace, insert the extra subspace l that covers the border _j, come mark with j, insert altogether Σ l _jThe extra subspace of height covers the border;

Step S1.3.6: the crucial subspace L=Σ l that finally obtains the joint mode sampling _i∪ Σ l _j

The beneficial effect of described step S1.3.1 is: consider the different three-dimensional structures of target receptor, avoid the calculating skewed popularity that brings because of textural difference, namely avoid missing important results in the channel sample process in next step.

The beneficial effect of described S1.3.5 is: prevent from missing the ligand molecular of the key that boundary may exist-target receptor binding mode.

Described step S2 comprises following concrete steps, as shown in Figure 4:

Step S2.1: at every sub spaces l _i∈ L puts into a ligand molecular, docks every sub spaces l _iObtain a 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 k target receptor-ligand interaction pattern), be specially: adopt molecular docking method analog sampling ligand molecular at all Interactions Modes of this subspace and target receptor, for obtaining accurately and the abundant Interactions Mode of diversity, the present invention's design is also developed one efficiently, novel molecular docking method, this docking calculation 1) in the sampling process to the prediction of binding mode with estimate and measure in conjunction with free energy computing function (MM-GBSA/MM-PBSA) so that this paper is improved, each energy term in this function is according to the feature restructuring and consider the target receptor conformation change and the relative Gibbs free that causes; 2) simultaneously according to these free energy computing method, develop multiple goal molecular docking method in conjunction with Multipurpose Optimal Method; 3) for taking into account docking result's accuracy and computation rate, target receptor is carried out the rigid and flexible zone divide processing.Beneficial effect is: improved the accuracy of prediction ligand molecular and target receptor binding mode, improved the accuracy in conjunction with free energy of estimating each Interactions Mode;

Particularly, 1) improvement of free energy function is characterised in that: each energy term is reassembled as four or a plurality of display items according to feature: such as Δ E _Vdw, Δ E _Es, Δ G _Gb, Δ G _Sa; Say something for more clear, the present invention illustrates that the Reorganization Energy quantifier also can not only be confined to above-mentioned four kinds as example take above-mentioned four, so long as split based on energy term, and the multiple goal docking calculation of design all belongs to protection domain of the present invention

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

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

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

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

Δ E _VdwBe the van der Waals energy quantifier, Δ E _EsBe the electrostatic energy quantifier, Δ G _GbBe solvent polarization free energy, Δ G _SaBe the non-polarized free energy of solvent; ω ₁-ω _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, and L is the energy term of ligand molecular self; Δ E _{Vdw, RL}For change when the identification interaction model that causes target receptor-ligand molecular of the conformation of target receptor and ligand molecular gets the variable quantity of Huaneng Group quantifier, Δ E _{Es, RL}Be change when the identification variable quantity of the interaction electrostatic term that causes target receptor-ligand molecular of the conformation of target receptor and ligand molecular; Respectively to be changed when the identification and the van der Waals energy quantifier of the target receptor that causes and the variable quantity of electrostatic term by the conformation of target receptor and ligand molecular; Respectively to be changed when the identification and the van der Waals energy quantifier of the ligand molecular that causes and the variable quantity of electrostatic term by the conformation of target receptor and ligand molecular;

With

To be changed when the identification and the solvent polarity that causes by the conformation of target receptor and ligand molecular

With nonpolar free energy

With

σ ₁=0.025 σ ₂=0.015 is the conversion constant of area and the kcal/mol of energy unit.

2) the multiple goal docking calculation of design is characterised in that: the flexible docking method that designs based on Multipurpose Optimal Method (such as multi-objective genetic algorithm, multi-objective Evolutionary Algorithm and multi-objective particle swarm algorithm etc.) and free energy computing method.Compare with traditional docking calculation, the beneficial effect of this docking calculation is: the diversity of the Interactions Mode that can increase substantially samples obtains, and accuracy and the corresponding evaluation accuracy in conjunction with free energy of raising sampling pattern.Wherein, the structure for objective function comprises following two schemes:

With Δ E _Vdw, Δ E _Es, Δ G _Gb, Δ G _Sa... a plurality of targets or fitness function (claim 6 is seen in each energy term definition, also can have other energy term further to expand this multiple objective function system simultaneously) as docking; Perhaps

With Δ E _Intra, Δ E _Inter, Δ G _Gb, Δ G _Sa... a plurality of targets or fitness function as docking;

Wherein, Δ E _VdwBe the van der Waals energy quantifier, Δ E _EsBe the electrostatic energy quantifier, Δ G _GbBe solvent polarization free energy, Δ G _SaBe the non-polarized free energy of solvent; Δ E _IntraBe the energy term of target receptor, ligand molecular self, Δ E _InterBe the non-bonded interaction between target receptor-ligand molecular, and the target numbers that adopts of optimizing process adopts elasticity to set system, namely between 1n, select.

3) when docking adopts the rigidity and tenderness zone to divide to process to target receptor to be characterised in that: the simulation ligand molecular in described subspace with the process of all Interactions Modes of target receptor in, adopting the rigidity and tenderness zone to divide to target receptor processes, wherein, flexible region is considered target receptor part-structure changeability, such as the motion change of protein Key residues or domain (such as the loop district), and adopt based on the energy scoring strategy of net point and atom pair and process described flexible region and rigid region.Flexible region can be considered the motion change in protein Key residues or ring district, and adopts based on the energy scoring strategy of net point and atom pair and process this two kinds of zones.Beneficial effect is: improved computing velocity; Process the errors of calculation that adopt two kinds of scoring strategies to bring for avoiding because of the division of rigidity and tenderness zone simultaneously, the unified employing of Interactions Mode lattice point scoring strategy is marked to process again and is solved this hidden danger the most at last.

Step S2.2: the receptor-ligand binding pattern of collecting in all subspaces obtains P=∑ P _i=∑ ∑ p _kBe specially: 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 that makes up BFEL; Beneficial effect is: system has reflected ligand molecular and all possible mode of action of target receptor all sidedly;

Step S2.3: the coordinate axis of definition ligand molecular and target receptor reaction, with the Interactions Mode p of each ligand molecular-target receptor _kReaction coordinate axle (x, y) in ligand molecular-target receptor reaction carries out projection, each Interactions Mode p _k(x, y) makes up element as target receptor-ligand molecular in conjunction with potential energy surface, and (the reaction coordinate axle can be various dimensions, and for simplicity, the present invention is take two dimension as example, and such as the two-dimentional reaction coordinate axle x that defines herein, y is such as p _k(x, y) represents the Interactions Mode p that obtains ligand molecular-target receptor that the S2.2 step is calculated _kAt x, the projection of y axle).Beneficial effect is: the reflection ligand molecular is combined with target receptor or the progress of dissociating and the instantaneous conformation change of part and target receptor, thereby can clearly reflect intermolecular mechanism of action;

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

Described step S3 further comprises two kinds of strategies:

(a) do not consider the reaction backtracking phenomenon of ligand molecular and target receptor, based on the route searching strategy of response path minimum energy principle;

(b) reaction backtracking phenomenon and the Direction of Reaction random perturbation of consideration ligand molecular and target receptor are based on the route searching strategy of response path minimum energy principle.

Strategy (a) may further comprise the steps, as shown in Figure 5:

Step S3.1.1: coarseness is divided potential energy surface

Ω represents whole BFEL,

Represent respectively whole potential energy surface according to the reaction coordinate axle x that defines among the step S2.3, the lattice point number of the upper division of y, n, m are respectively the lattice point size of space, determine reaction coordinate axle x, the upper lattice point number of y; Be specially: owing to obtain often more complicated of potential energy surface, thus strengthened the computing difficulty, be the quick obtaining response path therefore, at first original potential energy surface is simplified processing, namely BFEL is carried out coarseness and divide

So that the potential energy dimensionality reduction, thereby reducing computation complexity, beneficial effect is for improving counting yield;

Step S3.1.2: selected response path reference position S ₀(x _i, y _j),

x _iRepresent i lattice point on the reaction coordinate x, y _jRepresent j lattice point on the reaction coordinate y.Initialisation reactions path S:=S ∪ S ₀Concrete: selected ligand molecular and the dissociation reaction of target receptor or the reference position of association reaction, can be potential energy surface overall situation minimum point such as the reference position of dissociation reaction, the association reaction reference position can be in part and the solution effects zone certain a bit;

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 _i); S:=S ∪ S _P+1Obtain the response path of ligand molecular and target receptor in the enterprising row iteration search of the BFEL of coarseness according to the minimum energy principle, wherein each present position of ligand molecular all faces 5 advanceable directions, next step moving direction of ligand molecular is in 5 directions in conjunction with the minimum direction of free energy, wherein S is final way to acquire, S _pRepresent the p node location of step on the potential energy surface that coarseness is divided of ligand molecular and target receptor, p is the positive integer more than or equal to 0;

Step S3.1.4: judge present position S _p(if such as the association reaction of research ligand molecular and target receptor, then reaction end should be catalysis on the target receptor or the avtive spot of binding partner molecule whether to arrive user-defined reaction end; If the dissociation reaction of research ligand molecular and target receptor, then reaction end should be the zone away from target receptor, is generally the solvent zone.), be then to enter step S3.1.5, otherwise change step S3.1.3 over to;

Step S3.1.5: the ligand molecular that iteration is obtained and the coordinate information of target receptor response path S (are corresponding x, the position of y axle) backtracking is to the BFEL that does not carry out the coarseness division, obtain ligand molecular and target receptor response path and in conjunction with the funtcional relationship f (S) of free energy, according to the response path coordinate, revert to totipotency face (namely not carrying out the BFEL that coarseness is divided), thus the corresponding potential energy information in Regeneration response path: f (S);

Strategy (b) may further comprise the steps, as shown in Figure 6:

Step S3.2.1: BFEL is carried out coarseness divide

With step S3.1.1;

Step S3.2.2: with step S3.1.1, selected response path reference position, S ₀(x _i, y _j),

Initialisation reactions path S:=S ∪ S ₀

Step S3.2.3: judge present position S _pWhether arriving defined reaction end (with step S3.1.4), is then to enter step S3.2.12, otherwise changes step S3.2.4 over to;

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 is with step S2.3), return node ordered queue min_S _Vec, prepare for the subsequent reactions track search, enter step S3.2.5;

Step S3.2.5: each the present position S that checks the part target _p(x _i, y _i) site, border of potential energy surface whether, if not then change step S3.2.6 over to, if then change step S3.2.7 over to;

Step S3.2.6: the present position of judging part is which kind of the morphologic characteristics site in the potential energy surface, these morphologic characteristicss mainly are divided three classes: the potential well position, position, mountain peak and climbing position, if instant site is potential well position or position, mountain peak, the random number that we produce according to program decides next step moving direction of ligand molecular, random number adopts [0,1] random function produces, if value is in [0,0.9], be in the large probability situation, ligand molecular is according to the orderly sequence node min_S that returns among the step S3.2.4 _VecSelect, change step S3.2.11 over to; When if the random number value is in [0.9,1.0], namely in the small probability situation, the reaction working direction of ligand molecular namely changes step S3.2.12 over to for selecting at random;

Step S3.2.7: judge S _p(x _i, y _i) adjacent node in potential energy value minimum position whether be included among ligand molecular and the target receptor response path S, if this position then changes step S3.2.8 over to not in response path S; If being included in, this position explored in the response path that obtains, and this position is more than the appearance once in response path, then for avoiding the track search process to be absorbed in endless loop because pacing up and down in this position, the number of times that this position is occurred in response path limits.If some positions occur M time in response path, we just think that this position can cause endless loop, (M is the positive integer greater than 0, can look tolerance by the user defines, more close to 1, then redundant position is then fewer in the response path that obtains of whole exploration, reflects " making a circulation " feature that occurs in ligand molecular and the target receptor cohesive process greater than 1, is set as in embodiments of the present invention 3).If the above-mentioned definition that causes endless loop is not satisfied in this position, then change step S3.2.8 over to; Change step S3.2.9 over to if satisfy the endless loop definition;

Step S3.2.8: with present node S _pNeighbor node ordered queue min_S _VecIn the minimum node location min_S of potential energy _Vec[j] is as next step moving direction (S:=S ∪ min_S _Vec[j]), more new route present position (p=p+1) changes step S3.2.3 over to;

Step S3.2.9: facing in order the min_S of node queue at random _VecWorking direction (min_S of middle selection _Vec[j]), new route present position (S:=S ∪ min_S more _Vec[j], p=p+1) and change step S3.2.3 over to;

Step S3.2.10: if S _p(x _i, y _j) be the climbing node in the potential energy surface, then with S _p(x _i, y _j) be that the boundary is with neighbor node ordered queue min_S _VecTwo minutes, free energy numerical value was greater than S _p(x _i, y _j) the node composition sequence be L_S _VecLess than S _p(x _i, y _j) the node composition sequence be S_S _VecThe random number that we produce according to program decides next step moving direction of ligand molecular, random number adopts [0,1] random function produces, for the response path that guarantees to obtain as far as possible for free energy in the response path that might exist is minimum, the present invention is set in the greater probability situation that (random number that namely obtains is in a less scope time, the present invention is set as [0.3,1] time) the response path direction is still according to the minimum energy path direction, namely changes step S3.2.11 over to, at S_S _VecSelect in the sequence; Otherwise, produce the at random effect of perturbation of response path, namely change step S3.2.12 over to, 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 that selects minimum _Vec[j] or S_S _Vec[j] with step S3.2.7, is absorbed in endless loop because pacing up and down in this position for avoiding the track search process as next step the Direction of Reaction of ligand molecular, and the number of times that this position is occurred in response path limits.Can look tolerance by the user and define, test case of the present invention is made as 8; If this position is not satisfied and is above-mentionedly caused that the definition of endless loop then changes step S3.2.12 over to, otherwise with the reaction site of this position as next step, 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 select at random 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 changes step S3.2.3 over to;

Step S3.2.13: the ligand molecular that iteration is obtained revert to the BFEL that does not carry out the coarseness division with the coordinate information of the response path S of target receptor, obtain ligand molecular and target receptor response path and in conjunction with the funtcional relationship f (S) of free energy: namely, with the response path S that determines, according to its response path coordinate information, revert to the BFEL that does not have the coarseness processing, obtain accurate response path and get funtcional relationship in conjunction with free energy: f (S), for thermodynamic power of calculating the reaction generation provides foundation;

Step S4: calculate corresponding thermodynamics, kinetic parameter according to the response path that obtains.In this step, at first according to the ligand molecular that obtains and the response path of target receptor, determine the avtive spot (f (S) minimum position) in the course of reaction, transition state site (f (s) maximum position) and reaction terminating site (by step S3.1.4 definition).According to these three sites in conjunction with free energy data Δ G _Bind, Δ G _Trans, Δ G _Unbind,, obtain in conjunction with free energy in conjunction with transition state theory $({\mathrm{\ΔG}}_{\mathrm{binding}}^o=\mathrm{\Δ}$ $G_{\mathrm{bind}}-{\mathrm{\ΔG}}_{\mathrm{unbind}}),$ Association reaction energy of activation $({\mathrm{\ΔG}}_{\mathrm{on}}^{\≠}={\mathrm{\ΔG}}_{\mathrm{trans},}-{\mathrm{\ΔG}}_{\mathrm{unbind}})$ With dissociation reaction energy of activation $({\mathrm{\ΔG}}_{\mathrm{off}}^{\≠}=\mathrm{\Δ}$ $G_{\mathrm{trans},}-{\mathrm{\ΔG}}_{\mathrm{bind}}),$ Thereby obtain to obtain in conjunction with free energy according to transition state theory

Association reaction energy of activation With dissociation reaction energy of activation

And then calculate to obtain thermodynamic parameter according to Aileen's formula (such as formula 6-8): binding constant KD and kinetic parameter: association rate constant k _OnWith the rate constants k of dissociating _Off

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

ΔG _binding＝RTlnK _D (7)

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

K _BBe Boltzmann constant, R is gas law constant, and T is absolute temperature, and h is Planck constant.

The beneficial effect of described step S4 is: Accurate Prediction has been estimated the combination of drug molecule and target receptor or the difficulty of dissociation reaction mode, bond strength and association reaction and dissociation reaction, indirectly obtained the action time of drug molecule and target receptor, i.e. (1/k _Off).

Step S5: obtain molecular reaction mechanism, verify the crucial binding site that determines the reaction kinetics feature.Particularly: can revert to each present position S that the free energy potential energy surface is obtained ligand molecular on this path by the response path that obtains _p(x _p, y _p) reaction coordinate (x _p, y _p), namely can find the ligand molecular that is projected in this position among the step S2.3 and the binding pattern P of target receptor by reaction coordinate _k(x _i, y _i), x here _p=x _i, y _p=y _iTherefore can obtain ligand molecular at active binding site, the transition state site and binding pattern target receptor, learn avtive spot on the corresponding target receptor and transition state site (being the binding site on the target receptor passage in step S1.3.4, mentioned of ligand molecular), so can design with target receptor avtive spot have lower in conjunction with free energy (namely

Lower), and higher compound in conjunction with free energy is arranged (namely in the transition state site with target receptor

Higher), this newly-designed compound and target receptor not only have stronger affinity, have simultaneously to have long action time, thereby have stronger compound activity.

The beneficial effect of described step S5 is: only pay close attention to medicine (part) molecule and target receptor must find or design in conjunction with situation medicine at avtive spot one-sidedness in the design of breakthrough conventional medicament, provide the transition state site that determines medicinal effectiveness, for drug design provides new key message and design considerations.Particularly, the avtive spot on the acquisition target receptor and transition state site are for instructing medicine (part) design to play crucial effect; Show: can by rough Appraising Methods probe into other ligand moleculars and target receptor avtive spot in conjunction with the free energy situation, also can estimate simultaneously ligand molecular and target receptor transition state in conjunction with the free energy situation, thereby thereby can obtain roughly by calculating drug molecule the kinetic property of ligand molecular in the binding ability in these two kinds of sites, and then can by design can reduce its with acceptor avtive spot be combined in conjunction with free energy, improve it and acceptor is combined in the transition state site improves the activity of compound in conjunction with free energy.Thereby overcome in conventional medicament (part) design and be confined to avtive spot is ignored its kinetic property in conjunction with free energy shortcoming.Obtain molecular reaction mechanism of action, verify the crucial binding site that determines the reaction kinetics feature, as medicine (part) design key reference factor.

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

The first module, it carries out described step S1, namely build ligand molecular and target receptor simulated system, define ligand molecular and its space structure that may act on and range of size according to the crystal three-dimensional structure of the target receptor that parses by experiment and be used as simulating the research range of ligand molecular and target receptor Interactions Mode, the ligand molecular and the target receptor that adopt the method pre-service that is used for preparing ligand molecular and target receptor and prepare to simulate in advance;

The second module, it carries out described step S2, make up the BFEL of ligand molecular and target receptor, simulate the Interactions Mode of exhaustive ligand molecular and target receptor and assess these binding modes accordingly in conjunction with free energy, thereby systematically make up ligand molecular and target receptor interaction BFEL;

The 3rd module, it carries out described step S3, obtains the response path (Pathway) that ligand molecular and target receptor combination occur or dissociate according to BFEL;

Four module, it carries out described step S4, calculates corresponding thermodynamics and kinetics parameter according to the response path that obtains;

The 5th module, it carries out described step S5, determine the ligand molecular configuration of each present position of ligand molecular in whole response path by the response path that obtains, revert to the Interactions Mode that docking obtains, obtain ligand molecular at active binding site, the structural information in transition state site and in the conformation information of ligand binding acceptor during in corresponding site, determine avtive spot and transition state site on the target receptor, and then when the design drug molecule, can design and when avtive spot is combined, to have relatively high in conjunction with free energy with acceptor, and in the transition state site in conjunction with the time have a lower molecule in conjunction with free energy, namely not only improve the binding affinity of drug molecule and target receptor, and the hold-up time of prolong drug molecule and target receptor effect, thereby the biologically active of the drug molecule that raising is designed.

Description of drawings

In order to be illustrated more clearly in the embodiment of the invention or technical scheme of the prior art, the below will be briefly described Figure of description.

Fig. 1 is the trunk process flow diagram of the method for the invention.

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

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

Fig. 4 is step S2(simulation ligand molecular and the mode of action of target receptor, structure ligand molecular and target receptor in conjunction with the free energy potential energy surface) process flow diagram.

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

The backtracking phenomenon that Fig. 6 may exist when considering the reaction of ligand molecular and target receptor, and consider simultaneously the random perturbation of the Direction of Reaction, based on the method for searching path flow process (step S3, strategy (b)) of steepest descent theory.

Fig. 7 is the combination of ligand molecular and the interactional space of target receptor of research being carried out subspace division and possible ligand molecular or the synoptic diagram of dissociation channel sampled result, and the case study on implementation (HupA-TcAChE) in the present invention is as the example explanation.

Combination or dissociation channel that Fig. 8 obtains according to sampling select to obtain the result schematic diagram of crucial subspace, and the case study on implementation in the present invention is as the example explanation.

Fig. 9 is that the present invention plants all ligand moleculars of calculating among the embodiment and the binding mode space distribution situation of target receptor, is distributed as the sign synoptic diagram of representative with the configuration space of ligand molecular.

Figure 10 is interactional in conjunction with free energy potential energy surface three-dimensional representation for calculating the HupA-TcAChE that obtains.

Figure 11-A is the two-dimensional representation by the association reaction path of the HupA-TcAChE that adopts response path search strategy (a) acquisition, (because the present invention does not have artificial interference to the molecular action mode in whole simulation process, thereby can think that ligand molecular is identical with the dissociation reaction path with acceptor molecule association reaction path).

Figure 11-B be on the HupA-TcAChE association reaction path that obtains among Figure 11-A in conjunction with the Gibbs free curve.

Figure 12 is for adopting the ligand molecular that the response path heuristic approach obtains in the route searching strategy (b) to be combined with acceptor molecule/dissociation reaction path synoptic diagram.

The path that Figure 13 reacts for the ligand molecular that obtains by prediction and target receptor, determine critical sites (avtive spot and the transition state site) synoptic diagram of thermodynamics of reactions character and kinetic property on the target receptor that obtains, wherein grey cartoon figure is target proteins TcAChE, Dark grey club model is and the Key residues of the little molecular action of part that black mallet shape model is ligand molecular HupA.

Figure 14 is according to the critical sites that obtains, and carry out drug design strategy synoptic diagram by these two sites, namely have so relatively highly active compound by designing: its with acceptor avtive spot be combined relatively low in conjunction with free energy, with acceptor be combined in the transition state site higher in conjunction with free energy, thereby increase the affinity of drug molecule, increase drug molecule and target receptor effect hold-up time, improve drug effect.

Implementation and application mode

Below in conjunction with the embodiment of the invention and accompanying drawing, the technical scheme in the embodiment of the invention is clearly and completely described, obviously, described embodiment only is one embodiment of the invention, rather than whole embodiment.Based on the embodiment among the present invention, those of ordinary skills belong to the scope of protection of the invention not making the every other embodiment that obtains under the creative work prerequisite.

The invention discloses a kind of method of calculating the thermodynamics and kinetics parameter of this reaction of prediction by simulation ligand molecular and target receptor reactive mode, its natural law according to simulation ligand molecular-target receptor interaction pattern makes up BFEL, search for accordingly the response path of ligand molecular and target receptor, determine possible response path, thereby obtain the thermodynamics and kinetics parameter of ligand molecular and target receptor effect, obtain to determine on the acceptor ligand molecular and the avtive spot of its bond strength or action time and transition state site, therefore obtain to have so relatively highly active compound by design: its with acceptor avtive spot be combined relatively low in conjunction with free energy, with acceptor be combined in the transition state site higher in conjunction with free energy, thereby increase the affinity of drug molecule, increase drug molecule and target receptor effect hold-up time.As target receptor, huperzine ((-)-HuperizineA, HupA) is as ligand molecular with California electricity precious jade acetylcholinesterase (TcAchE) for present embodiment.But it will be appreciated by persons skilled in the art that method of the present invention can be used for assessing or predicting quickly and accurately thermodynamics of reactions and the kinetic parameter of other any ligand moleculars and any target receptor equally.Its embodiment is as described below:

The disclosed exhaustive ligand molecular of a kind of simulation of the embodiment of the invention and the target receptor mode of action, the prediction ligand molecular is combined with target receptor/dissociation reaction, and the trunk flow process of the thermodynamics and kinetics parametric technique of prediction and calculation reaction is as shown in Figure 1, comprising:

Step S1: according to the crystal structure of testing the target receptor that obtains, define the spatial dimension that acquisition ligand molecular and target receptor may act on, and build the analog computation system according to the crystal structure of target receptor and the action space scope of definition.Concrete build process comprises as shown in Figure 2:

Step S1.1: the bulk scope of definition counting system, the definition size should be noted: enough large, not only can cover the acceptor three-dimensional structure, also provide simultaneously enough solvents space not being subjected to away from receptor protein in the situation of its function influence for ligand molecular, free movement in solution also is unlikely to too large and meaningless waste computational resource simultaneously.TcAChE system for this example will be studied adopts Cube, as shown in Figure 7, because the TcAChE system is less, therefore, this range of size is enough, also can be other sizes for other systems certainly, depending on research system and need to and decide with calculating, also let loose in cube simultaneously, can be other space forms such as spheroid equally.

Step S1.2: the computer memory to definition carries out subspace grid division, as shown in Figure 7, be parallelization joint mode sampling process, avoid simultaneously the optimization method global convergence and miss the key effect pattern of intermediates, the present invention proposes that the subspace division is carried out in the space of the research system of definition and processes, 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.At this, for concrete division methods, will in this example

Cube is divided into isopyknic 1000 small cubes, certainly be not limited to the size of this example for concrete small cubes size, simultaneously also and do not require that the small cubes size wants consistent.

Step S1.3: select crucial subspace.Since exist in the research space of above-mentioned steps definition some ligand moleculars can't in conjunction with scope (such as the tightly packed zone of acceptor residue (base)) or the solvent space excessively far away apart from acceptor, ligand molecular is combined with acceptor molecule or thermodynamic parameter and the kinetic parameter of dissociation reaction there is no significant contribution for obtaining in these zones, thereby therefore can ignore computing cost.Concrete selection scheme following (such as Fig. 3):

Step S1.3.1: collect all target receptor crystal structures, for same target receptor since the binding partner molecule before and after, or even all may induce its conformation that different variations occurs in conjunction with different ligand moleculars, passage or binding site for the ligand molecular Interactions Mode that different receptor conformations exists then are not quite similar, therefore, the present invention proposes to avoid missing critical passage or binding site by the method for studying many receptor conformation passages.

Step S1.3.2: (present case obtains the crystal structure of 54 TcAChE altogether for research, such as table 1 for the different target receptor structure of obtaining.), adopt third party software, such as MOLE

Otyepka M (2007) MOLE:A Voronoi diagram-based explorer of molecular channels, pores, and tunnels.Structure15: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 andcavities.BMC Bioinformatics 7:316.) etc. search for the important channel ti (represent the passage that i structure obtains, number is indefinite) that exists on each structure, all critical passage T=∑ ti of enrichment;

Table 1

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

Step S1.3.3: all passage T of cluster obtain representative passage T _rOwing to have certain plyability between the passage that obtains, be to simplify to calculate, obtain passages in this step for these and carry out cluster according to similarity, thereby finally obtain the passage T that selects _r=8, such as accompanying drawing 7;

Step S1.3.4: according to the passage T that obtains _rSelect crucial subspace.In the research range of definition, cover or adjacent passages T _rSubspace l _iAll optionally do crucial subspace, as follow-up Interactions Mode sample space.

Step S1.3.5: process the border between the crucial subspace that obtains.For avoiding missing the key interactions pattern of adjacent subspace boundary, the present invention proposes the extra method of processing the border, subspace of improving, and namely the boundary every pair of subspace additionally appends a series of subspace l _i(size can be different from the subspace of initial setting, for ease of operation, the big or small subspace such as adopt in the example of the present invention) cover the border, as extra unit, sample space, as shown in Figure 8.

Step S1.3.6:L=Σ l _i∪ Σ l _i=72 are the subspace of choosing as the joint mode sampling.

For step S1.3, in order to make up more fully potential energy surface, can also ignore this step, all the subspace is all as computer memory.

Step S2: according to analog sampling ligand molecular and the interactional natural law of target receptor obtain its Interactions Mode and to Interactions Mode calculate assessment in conjunction with free energy, thereby make up the BFEL of ligand molecular and target receptor Interactions Mode; This step mainly comprises following idiographic flow (such as Fig. 4):

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 a target receptor-ligand molecular Interactions Mode set P _i=∑ p _k

For making up potential energy surface, the evaluation accuracy in conjunction with free energy of the target receptor-diversity of ligand molecular Interactions Mode, the accuracy of Interactions Mode and each Interactions Mode is most important.And common docking calculation adds owing to the rough property of scoring functions and when estimating in conjunction with energy and the Interactions Mode diversity that obtains of causing sampling is poor, thereby causes above-mentioned three aspects to guarantee, makes up accurately this operation of BFEL thereby can't finish.

The present invention is directed to these limitations, propose a series of innovative approach, effectively avoid these problems, be respectively:

1) to the improvement of free energy evaluation method MM-GBSA or MM-PBSA, comprises for each energy term according to feature restructuring and consider the target receptor conformation change and the relative Gibbs free that causes, thereby improved free foreseeable accuracy.In the present embodiment; be described as example to be reassembled as four in conjunction with free energy; energy project restructuring number also can be less than four or more than four, and number do not limit, and comes developer molecule docking calculation all to belong to protection scope of the present invention with energy term as objective function based on this.

Each energy term is reassembled as four according to feature: Δ E _Vdw, Δ _Ee _s, Δ G _Gb, Δ G _Sa;

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

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

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

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

Δ E _VdwBe the van der Waals energy quantifier, Δ E _EsBe the electrostatic energy quantifier, Δ G _GbBe solvent polarization free energy, Δ G _SaBe the non-polarized free energy of solvent; ω ₁-ω ₄Be the weight of each energy term, obtain according to the match empirical value; RL, R, L are represented as respectively the Interactions Mode between target receptor-ligand molecular compound, and R is the energy term of acceptor self, and L is the energy term of part self. With

To be changed when the mutual identification by the conformation of target receptor and ligand molecular, thus the Van der Waals that causes and the variation of electrostatic energy quantifier; Equally,

With

To be changed when the mutual identification by the conformation of target receptor and ligand molecular, thus the solvent polarity that causes and nonpolar free energy; For weights omega ₁-ω ₄Acquisition, the present invention proposes in view of the thermodynamic data of having reported

K _A, K _D) and crystal structure data (such as table 2) come the Fitting Calculation result and obtain, be respectively ω ₁=0.2245, ω ₂=-0.4818, ω ₃=-0.6288 and ω ₄=-0.6288, σ ₁=0.025, σ ₂=0.015;

Table 2

Table 2 is for being the collected different ligand molecular of match MM-GBSA and the experimental data of TcAChE in the example of the present invention, and the scoring details by original MM-GBSA.PDB represents crystal structure database, ^aThe Medicine small molecule that representative obtains from Research Literature and the dissociation constant of acceptor, ^bRepresentative by arriving of obtaining that dissociation constant transforms in conjunction with free energy Δ G _Binding-exp=RTlnK _D

2) directly adopt improved free energy computing method to instruct the docking sampling process as scoring functions, thus the accuracy of the Interactions Mode that assurance is obtained;

3) diversity of horn of plenty joint mode, prevent in traditional docking calculation, owing to missing the diversity of joint mode (shown in equation 1 in conjunction with adding with unicity of energy, the free energy that ligand molecular is combined with receptor protein is often obtained by a series of energy project stacks, therefore final energy is close, often but there are differences between each energy project, and this species diversity causes owing to the difference of Interactions Mode often, traditional docking calculation since only with final in conjunction with energy as evaluation criterion, thereby tend to lose a lot close in conjunction with energy, and the Interactions Mode that structure there are differences), therefore, the present invention's proposition respectively as objective function, is adopted each energy term the strategy of multiple-objection optimization docking carry out the joint mode sampling work, thereby is avoided missing these Interactions Modes;

4) in addition, losing crucial local superiority for the convergence of further avoiding optimization method separates, the present invention is in traditional Multipurpose Optimal Method, integrated elite's solution archives preservation method, namely in optimizing process, set up archives " container ", be used for preserving all local superiorities' solutions that whole optimizing process obtains, thereby further avoid missing the advantage solution that causes because of convergence.

Except addressing the above problem, the present invention is for improving docking speed, in docking operation, merge based on lattice point (Zou XQ, SunYX, Kuntz ID (1999) Inclusion of Solvation in Ligand Binding Free Energy Calculations Usingthe Generalized-Born Model.J Am Chem Soc 121:8033-8043.Meng EC, Shoichet BK, KuntzID (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 ofsolute 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 incalculating conformational energies of organic and biological molecules J Comput Chem 21:1049-1074.) these two kinds methods of estimating ligand molecular target receptor interaction force, part for target receptor conformation mutability is defined as flexibility, and other parts are defined as rigidity, use the energy variation of calculating the interaction force between itself and the ligand molecular and self causing because of conformation change based on the atom pair Appraising Methods in docking operation for flexible portion; And rigid element is used the interaction force that calculates itself and ligand molecular and acceptor flexible portion based on the evaluation method of lattice point in docking operation, in this example flexible residue is defined as at the residue of predicting on the reaction channel, the definition with residue number and residue rule of thumb is worth setting by the user.

In the Interactions Mode sampling process, in each lattice, put into a ligand molecular, ligand molecular is free to carry out variation and the translation rotation of conformation in each lattice, calculating between subspace and the subspace does not intersect and affects, therefore this sampling process can be the task number that equates with the subspace number according to the number assignment of subspace, carry out parallel computation, significantly reduce computing time.In the docking sampling process to every sub spaces, at first adopt the flexible docking method to carry out initial Interactions Mode sampling, for further improving the diversity of Interactions Mode, merge so that calculate whole in conjunction with the power generation error for fear of above-mentioned two kinds (namely based on atom pair and based on lattice point) methods simultaneously, adopt again the method that acceptor is considered as rigidity fully to carry out multiple goal rigidity docking sampling, carry out sampling operation for the selectable lattice of institute in each receptor conformation, this process adopts the scoring method based on lattice point fully, and the Interactions Mode that obtains for flexible docking equally also adopts this lattice point marking form again to mark.So, being obtained by the method evaluation based on lattice point in conjunction with free energy is just unified of all joint modes, in enrichment in the Interactions Mode, avoided again the error that exists between the distinct methods and the inaccuracy that may cause.

Step S2.2: the Interactions Mode enrichment of the ligand molecular-target receptor of all subspaces is obtained Interactions Mode set P=∑ P _i=∑ ∑ p _k, obtain all data such as Fig. 9, P=127,371 such as present embodiment.

Step S2.3: definition reaction coordinate axle, with each target receptor-ligand molecular Interactions Mode according to reaction coordinate axle projection, each Interactions Mode p _k(x, y) then make up element as target receptor-ligand molecular Interactions Mode potential energy surface, in the present embodiment, the definition reaction coordinate is made of jointly two variablees, i.e. the instant configuration of the conformational difference of ligand molecular and ligand molecular and its centroidal distance at the configuration of avtive spot.Certainly for making up the multidimensional potential energy surface, the visual requirement of the number of reaction coordinate and increasing.

Step S2.4: by third party mapping software (such as matlab, R etc.), draw the BFEL of ligand molecular and target receptor Interactions Mode, such as Figure 10.Wherein x axle and y axle are the reaction coordinate axle, and the z axle is in conjunction with the free energy coordinate axis.

Obtained after the BFEL, then can calculate according to it the thermodynamics and dynamics parameter of ligand molecular and target receptor reaction.

Step S3: obtain ligand molecular and target receptor generation combination or dissociation reaction path according to BFEL; For the combination of ligand molecular and target receptor or the exploration of dissociation reaction, this research and design two kinds of strategies, strategy (a): part and receptor response process are not considered ligand molecular in the backtracking effect of course of reaction, and strictly follow the principle of response path minimum energy.As shown in Figure 5, concrete step is as follows:

Step S3.1.1: calculate for convenient, at first BFEL is carried out coarseness and divide,

N=m=100 in this example, x=[0.0,0.04], y=[0.0,35.0], certainly should divide yardstick and can take the circumstances into consideration to set according to the complexity of potential energy surface etc.

Step S3.1.2: selected response path reference position, S ₀(x _i, y _j),

Selected for the reference position of response path, in view of the natural interaction process of the potential energy surface that obtains among the present invention for simulation ligand molecular and target receptor, there is no the intervention of any human factor, therefore also can be regarded as in conjunction with the path simultaneously for the Dissociation path that obtains.Because the indefinability to initial position in the recognition reaction of ligand molecular and target receptor, the Dissociation path that the design explores ligand molecular is main research approach, for dissociation reaction, then react initial point and should be ligand molecular and target receptor at the bonding state of avtive spot, therefore for the system that crystal complex is arranged, then with the position of ligand molecular in the crystal complex as initial position, if there is not crystal complex information, look then overall potential energy can minimum site be the dissociation reaction reference position in the whole potential energy surface.After having set reference position, initialisation reactions path S:=S ∪ S then ₀

Step S3.1.3: ligand molecular is when reacting with target receptor, excessive conformation or the variation of configuration can not occur in ligand molecular within the continuous time, do not consider simultaneously the backtracking effect in the ligand molecular course of reaction, ligand molecular always faces 5 directions in dissociation process, then its working direction be in these five directions potential energy minimum a bit, 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.Iteration is sought until arrive reaction end by this method.

Step S3.1.4: reaction end is the iteration termination locations of route searching.In this example, with ligand molecular HupA present position apart from avtive spot

The place is as the terminal point of dissociation reaction, and this numeral can be according to the target receptor size, and the size of ligand molecular etc. are determined, the numerical value that does not adopt take this example is as standard, enough fully get final product so long as ligand molecular dissociates, criterion is away from target receptor, acts on fully with solution.

Step 3.15: the coordinate information backtracking that iteration is obtained response path S obtains the variation funtcional relationship f (S) in conjunction with free energy and path under the same coordinate to the BFEL that does not have coarseness, such as the result of this example acquisition, such as Figure 11 B.

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

Complement one another with strategy (a), considered not only in this programme that ligand-receptor is in conjunction with the backtracking effect in the dissociation reaction, added simultaneously the at random simulation operations of perturbation of small probability on the basis based on response path energy variation minimum theoretical, be that ligand molecular is in course of reaction, if if the different directions that faces in the next step of reaction may have the situation that energy is more or less the same, then may continue with the form of at random choice direction reaction.Idiographic flow such as accompanying drawing 5.Wherein important step is as follows:

Step S3.2.1: BFEL is carried out coarseness divide

N=m=100 is with step S3.1.1;

Step S3.2.2: selected response path reference position, i.e. potential energy surface overall situation extreme lower position, S ₀(x _i, y _j),

Initialisation reactions path S:=S ∪ S ₀, with step S3.1.2;

Step S3.2.3: judging whether to arrive defined reaction end, is then to enter step S3.2.12, otherwise changes step S3.2.4 over to

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 the subsequent path exploration, enter step S3.2.5;

Step S3.2.5: each the present position S that checks the part target _p(x _i, y _j) site, border of potential energy surface whether, if not then change step S3.2.6 over to, if then change step S3.2.7 over to;

Step S3.2.6: the present position of judging part is which kind of the morphologic characteristics site in the potential energy surface, and these morphologic characteristicss mainly are divided three classes: potential well position, position, mountain peak and climbing position.If instant site is potential well position or position, mountain peak, (be set as 90% in such as example of the present invention in greater probability, can be any bigger numerical) under, ligand molecular all should be selected according to the sequence node that returns among the S3.2.4, and concrete selection course changes step S3.2.11 over to.And remain under 10% the probability, the reaction working direction of ligand molecular is to select at random, namely changes step S3.2.12 over to.

Step S3.2.7: judge S _p(x _i, y _j) face whether potential energy value minimum position is lived through in the track search process in the node, or lived through but the approach number of times is not more than (this numeral can be decided for the tolerance of path multiplicity on own by the user) 3 times, if above-mentioned both of these case then changes step S3.2.8 over to, otherwise change step S3.2.9 over to.

Step S3.2.8: with present node S _pNeighbor node ordered queue min_S _VecIn the minimum node location min_S of potential energy _Vec[j] is as next step moving direction (S:=S ∪ min_S _Vec[j]), more new route present position (p=p+1) changes step S3.2.3 over to;

Step S3.2.9: facing in order the min_S of node queue at random _VecWorking direction (min_S of middle selection _Vec[j]), new route present position (S:=S ∪ min_S more _Vec[j], p=p+1) and change step S3.2.3 over to;

Step S3.2.10: if S _p(x _i, y _j) be the climbing node in the potential energy surface, then with S _p(x _i, y _j) be that the boundary is with neighbor node ordered queue min_S _VecTwo minutes, free energy numerical value was greater than S _p(x _i, y _j) the node composition sequence be L_S _VecLess than S _p(x _i, y _j) the node composition sequence be S_S _VecThe random number that we produce according to program decides next step moving direction of ligand molecular, random number adopts [0,1] random function produces, for the response path that guarantees to obtain as far as possible for free energy in the response path that might exist is minimum, the present invention is set in the greater probability situation that (random number that namely obtains is in a less scope time, the present invention is set as [0.3,1] time) the response path direction is still according to the minimum energy path direction, namely changes step S3.2.11 over to, at S_S _VecSelect in the sequence; Otherwise, produce the at random effect of perturbation of response path, namely change step S3.2.12 over to, 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 that selects minimum _Vec[j] or S_S _Vec[j] with step S3.2.7, is absorbed in endless loop because pacing up and down in this position for avoiding the track search process as next step the Direction of Reaction of ligand molecular, and the number of times that this position is occurred in response path limits.Can look tolerance by the user and define, test case of the present invention is made as 8; If this position is not satisfied and is above-mentionedly caused that the definition of endless loop then changes step S3.2.12 over to, otherwise with the reaction site of this position as next step, 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 select at random 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 changes step S3.2.3 over to;

Step S3.2.13: the ligand molecular that iteration is obtained revert to the BFEL that does not carry out the coarseness division with the coordinate information of the response path S of target receptor, obtain ligand molecular and target receptor response path and in conjunction with the funtcional relationship f (S) of free energy: namely, with the response path S that determines, according to its response path coordinate information, revert to the BFEL that does not have the coarseness processing, obtain accurate response path and get funtcional relationship in conjunction with free energy: f (S), such as accompanying drawing 12.

Step S4: calculate corresponding thermodynamics, kinetic parameter according to the response path that obtains; In this step, at first according to the response path that obtains ligand molecular and target receptor, determine the binding site in the course of reaction, transition state site and reaction terminating site and these three critical sites in conjunction with free energy numerical value, thereby obtain to obtain in conjunction with free energy according to transition state theory Association reaction energy of activation

With dissociation reaction energy of activation

Thereby can calculate the acquisition thermodynamic parameter according to Aileen's formula (such as formula 6-8) is binding constant K _DWith kinetic parameter, i.e. association rate constant k _On, dissociation rate constant k _Off

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

ΔG _binding＝RTlnK _D (7)

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

K _BBe Boltzmann constant, R is gas law constant, and T is absolute temperature (this example adopts 298.15K), and h is Planck constant.

Step S5: obtain molecular reaction mechanism of action, verify the crucial binding site that determines the reaction kinetics feature, as medicine (part) design key reference factor.By obtaining response path, can determine the configuration of each present position of ligand molecular in whole response path, revert to and dock the structural information details that the Interactions Mode that obtains namely can obtain course of reaction.Thus, ligand molecular be can obtain in active binding site, transition state site structural information, avtive spot and transition state site on the target receptor learnt, such as accompanying drawing 13 and Figure 14.

Obtaining of kinetic parameter is of great significance, because can obtain the hold-up time of part (medicine) molecule by dissociation rate constant, be the lasting medicine effect, for the pharmacological effect of studying part (medicine) molecule has played very crucial meaning.

In the present embodiment, take the TcAChE(target receptor)-HupA (ligand molecular) is example, adopt method proposed by the invention that the example system is carried out analog computation, thermodynamics and dynamics parameter when the HupA that obtains and TcAChE reaction, the numerical value that the numerical value that analog computation obtains extremely obtains near experiment test.Wherein, in conjunction with free energy Only differ about 1kJ/mol with experiment value, and kinetic parameter, namely

With

Compare with experiment value, also only poor 1.32kJ/mol and 3.61kJ/mol.This shows that method of the present invention can calculate binding affinity and the kinetic parameter of prediction ligand molecular and target receptor quickly and accurately, and then is expected to be generalized to medicine (part) design field and instructs excavation or design improvement lead compound.

Claims (12)

1. simulate the method that ligand molecular and target receptor reaction and calculating are predicted the thermodynamics and dynamics parameter of this reaction for one kind, said method comprising the steps of:

Step S1: download the described target receptor of by experiment means acquisition and the three-dimensional structure data of ligand molecular from online target receptor crystal structure database, build the simulated system of ligand molecular and target receptor based on the three-dimensional structure of ligand molecular and target receptor;

Step S2: the BFEL that makes up ligand molecular and target receptor;

Step S3: obtain the response path that ligand molecular and target receptor combination occur or dissociate according to BFEL;

Step S4: calculate corresponding thermodynamics and kinetics parameter according to the response path that obtains;

Step S5: obtain the mechanism of action of molecular reaction, definite crucial binding site that determines the reaction kinetics feature by response path.

2. method according to claim 1, wherein, described online target receptor crystal structure database is RCSB( Http:// www.pdb.org/pdb/home/home.do), described laboratory facilities are X-ray diffraction and nuclear-magnetism diffraction, described ligand molecular is Medicine small molecule or compound, and preferred huperzine, described target receptor is disease-associated protein or nucleic acid, preferred acetylcholinesterase.

3. method according to claim 1, wherein, described step S1 may further comprise the steps:

Step S1.1: the range of size in definition research system space, it is characterized in that, delimit sampling ligand molecular and the interactional bulk scope of target receptor intended by the three-dimensional structure of target receptor;

The interactional bulk scope of described sampling ligand molecular and target receptor is greater than 10 ~ 20 dusts of the three-dimensional structure of described target receptor;

Step S1.2: the subspace grid is carried out in the described research system space of step S1.1 definition divide, and every sub spaces is carried out mark, with l _iCome mark, represent the i sub spaces;

Step S1.3: the crucial subspace of screening joint mode sampling, it is characterized in that, the subspace that obtains among the step S1.2 is screened, select the subspace larger with the possibility of described ligand molecular effect, get rid of the subspace less with the possibility of described ligand molecular effect.

4. method according to claim 3, wherein, described step S1.3 may further comprise the steps:

Step S1.3.1: download the by experiment crystal structure data of the described target receptor of means acquisition from online target receptor crystal structure database, described online target receptor crystal structure database is RCSB( Http:// www.pdb.org/pdb/home/home.do), described laboratory facilities are X-ray diffraction and nuclear-magnetism diffraction, described crystal structure only is the target receptor monomer structure reaches the target receptor compound that coexists with little molecule without the apo crystal structure of little molecule coexistence holo structure;

Step S1.3.2: adopt the third party software that supplies the little molecule combination of part or dissociation channel calculating sampling on the target receptor, preferred MOLE or CAVER software calculate the potential passage t of each target receptor crystal structure _i, and all passage T=∑ t of enrichment _i, it is characterized in that, adopt described third party software, potential passage t on each target receptor of sampling successively _i, the passage that obtains on each target receptor is the most at last integrated T=∑ t _i,, wherein, the three-dimensional structure of i the target receptor that i representative is calculated, i depends on the three-dimensional structure number of the required calculating target receptor of user, for greater than 1 positive integer; t _iBe potential passage on i the structure that obtains by calculating; Preferably, adopt MOLE or the desired most possible number of active lanes of obtaining of the self-defined each calculating of CAVER software, i is defaulted as 3;

Step S1.3.3: all critical passage T of cluster obtain representative passage T _r, to analyze and obtain all passage T among the S1.3.2, the passage of each orientation only keeps one, and other filter for redundant;

Step S1.3.4: select all coverings or contiguous T _rSubspace l _iAs the joint mode crucial subspace of sampling, it is characterized in that, according to described passage T _rFrom all subspaces that step obtains S.1.2, filter out crucial subspace Σ l _i

Step S1.3.5: boundary treatment is carried out in described crucial subspace, it is characterized in that, between adjacent subspace, insert the extra subspace l that covers the border _j, come mark with j, insert altogether Σ l _jIndividual extra subspace covers the border;

Step S1.3.6: the crucial subspace L=Σ l that finally obtains the joint mode sampling _i∪ Σ l _j

5. method according to claim 1, wherein, described step S2 may further comprise the steps:

Step S2.1: at every sub spaces l _i∈ L puts into a ligand molecular, carries out the calculating of docking of ligand molecular and target receptor, every sub spaces l _iObtain a target receptor-ligand molecular Interactions Mode set P _i=∑ p _k, wherein, k is the positive integer more than or equal to 0, p _kRepresentative obtains k target receptor-ligand interaction pattern;

Step S2.2: the Interactions Mode of collecting the target receptor-ligand molecular that calculates in all subspaces among the step S2.1 obtains P=∑ P _i, it is characterized in that, 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 that makes up BFEL;

Step S2.3: the coordinate axis of definition ligand molecular and target receptor reaction, with the Interactions Mode p of each ligand molecular-target receptor _kReaction coordinate axle (x, y) in ligand molecular-target receptor reaction carries out projection, each Interactions Mode p _k(x, y) makes up element as target receptor-ligand molecular in conjunction with potential energy surface;

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

6. method according to claim 5, wherein,

Described step S2.1 is characterised in that, the simulation ligand molecular in described subspace with the process of all Interactions Modes of target receptor in, each energy term of free energy function is recombinated according to feature:

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

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

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

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

Δ E _VdwBe the van der Waals energy quantifier, Δ E _EsBe the electrostatic energy quantifier, Δ G _GbBe solvent polarization free energy, Δ G _SaBe the non-polarized free energy of solvent; ω ₁-ω _nBe the weight of each energy term, wherein, n is the positive integer greater than 1, preferred 4; RL is the interaction energy term of target receptor ligand molecular, and R is the energy term of target receptor self, and L is the energy term of ligand molecular self; Δ E _{Vdw, RL}For change when the identification interaction model that causes target receptor-ligand molecular of the conformation of target receptor and ligand molecular gets the variable quantity of Huaneng Group quantifier, Δ E _{Es, RL}Be change when the identification variable quantity of the interaction electrostatic term that causes the target receptor ligand molecular of the conformation of target receptor and ligand molecular;

Respectively to be changed when the identification and the van der Waals energy quantifier of the target receptor that causes and the variable quantity of electrostatic term by the conformation of target receptor and ligand molecular;

Respectively to be changed when the identification and the van der Waals energy quantifier of the ligand molecular that causes and the variable quantity of electrostatic term by the conformation of target receptor and ligand molecular;

With

To be changed when the identification and the solvent polarity that causes by the conformation of target receptor and ligand molecular

With nonpolar free energy

With

σ ₁=0.025 σ ₂=0.015 is the conversion constant of area and the kcal/mol of energy unit.

7. method according to claim 5, wherein,

Described step S2.1 is characterised in that, in described subspace with in the process of all Interactions Modes of target receptor, adopt the molecular docking method of developing based on multi-objective optimization algorithm at the simulation ligand molecular, described multiple goal docking calculation is the flexible docking method, in described flexible docking method, comprise following two schemes for the structure of objective function:

With Δ E _Vdw, Δ E _Es, Δ G _Gb, Δ G _Sa... a plurality of targets or fitness function as docking; Perhaps

With Δ E _Intra, Δ E _Inter, Δ G _Gb, Δ G _Sa... as a plurality of targets or the fitness function of docking,

Wherein, Δ E _VdwBe the van der Waals energy quantifier, Δ E _EsBe the electrostatic energy quantifier, Δ G _GbBe solvent polarization free energy, Δ G _SaBe the non-polarized free energy of solvent; Δ E _IntraEnergy term sum for target receptor, ligand molecular self comprises:

Δ E _InterInteraction force sum between target receptor-ligand molecular comprises: Δ E _{Vdw, RL}With Δ E _{Es, RL}, Δ G _Gb, Δ G _Sa,

The target numbers that adopts in the described objective optimization algorithm is that elasticity is set system, namely selects between 1-n, and n is user-defined objective function number, and n is the positive integer greater than 1.

8. method according to claim 5, wherein,

Described step S2.1 is characterised in that, the simulation ligand molecular in described subspace with the process of all Interactions Modes of target receptor in, adopting the rigidity and tenderness zone to divide to target receptor processes, wherein, flexible region is considered target receptor part-structure changeability, and adopts based on the energy scoring strategy of net point and atom pair and process described flexible region and rigid region; Preferably described target receptor part-structure changeability is the motion change of protein Key residues or domain; More preferably, described domain is the loop district.

9. method according to claim 5, wherein, described step S3 further comprises strategy (a) or strategy (b):

Strategy (a) may further comprise the steps:

Step S3.1.1: BFEL is carried out coarseness divide

Ω represents whole BFEL,

Represent respectively the reaction coordinate axle x of whole potential energy surface, the lattice point number of the upper division of y, n, m are respectively the lattice point size of space, determine reaction coordinate axle x, the upper lattice point number of y;

Step S3.1.2: selected response path reference position S ₀(x _i, y _j),

Xi represents i lattice point on the reaction coordinate x, y _jRepresent j lattice point on the reaction coordinate y, the 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+1Obtain the response path of ligand molecular and target receptor in the enterprising row iteration search of the BFEL of coarseness according to the minimum energy principle, wherein each present position of ligand molecular all faces 5 advanceable directions, next step moving direction of ligand molecular is in 5 directions in conjunction with the minimum direction of free energy, wherein S is final way to acquire, S _pRepresent the p node location of step on the potential energy surface that coarseness is divided of ligand molecular and target receptor, p is the positive integer more than or equal to 0;

Step S3.1.4: judge present position S _pWhether arriving user-defined reaction end, is then to enter step S3.1.5, otherwise changes step S3.1.3 over to, and wherein, when the association reaction of research ligand molecular and target receptor, reaction end is catalysis on the target receptor or the avtive spot of binding partner molecule; When studying the dissociation reaction of ligand molecular and target receptor, reaction end is the zone away from target receptor, and described zone away from target receptor is the solvent zone;

Step S3.1.5: the ligand molecular that iteration is obtained and the coordinate information backtracking of target receptor response path S are to the BFEL that does not carry out the coarseness division, obtain ligand molecular and target receptor response path and in conjunction with the funtcional relationship f (S) of free energy, according to the response path coordinate, revert to the totipotency face, thereby the corresponding potential energy information in Regeneration response path: f (S), wherein, the coordinate information of described S is corresponding x, the position of y axle, the BFEL that described totipotency face is divided for not carrying out coarseness;

Strategy (b) may further comprise the steps:

Step S3.2.1: BFEL is carried out coarseness divide With step S3.1.1;

Step S3.2.2: with step S3.1.1, selected response path reference position, S ₀(x _i, y _j), Initialisation reactions path S:=S ∪ S ₀

Step S3.2.3: judge present position S _pWhether arriving defined reaction end, is then to enter step S3.2.12, otherwise changes step S3.2.4 over to, with step S3.1.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 the subsequent reactions track search, enter step S3.2.5;

Step S3.2.5: each the present position S that checks the part target _p(xi, yi) be the site, border of potential energy surface whether, if not then change step S3.2.6 over to, if then change step S3.2.7 over to;

Step S3.2.6: the present position of judging part is which kind of the morphologic characteristics site in the potential energy surface, when instant site is potential well position or position, mountain peak, the random number that produces according to program decides next step moving direction of ligand molecular, random number adopts [0,1] random function produces, if value is in the large probability situation of [0,0.9], the orderly sequence node min_S that ligand molecular returns in step S3.2.4 _VecSelect successively, change following steps S3.2.11 over to; Be in when the random number value in the small probability situation of [0.9,1.0], the reaction working direction of ligand molecular namely changes following steps S3.2.12 over to for selecting at random; Wherein, described morphologic characteristics site comprises 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 be included among ligand molecular and the target receptor response path S, if this position then changes step S3.2.8 over to not in response path S; If being included in, this position explored in the response path that obtains, and this position is more than the appearance once in response path, then for avoiding the track search process to be absorbed in endless loop because pacing up and down in this position, the number of times that this position is occurred in response path limits; Wherein, when some positions occur M time in response path, think that then this position can cause endless loop, M is the positive integer greater than 0, and M looks tolerance by the user and defines, more close to 1, then redundant position is then fewer in the response path of whole exploration acquisition, reflect " making a circulation " feature that occurs in ligand molecular and the target receptor cohesive process greater than 1, preferably, M is 3; If this position occurrence number in response path then changes following steps S3.2.8 over to less than M; If this position occurrence number in response path then changes following steps S3.2.9 over to more than or equal to M;

Step S3.2.8: with present node S _pNeighbor node ordered queue min_S _VecIn the minimum node location min_S of potential energy _Vec[j] is as next step moving direction S:=S ∪ min_S _Vec[j], more new route present position p=p+1 changes step S3.2.3 over to;

Step S3.2.9: facing in order the min_S of node queue at random _VecWorking direction min_S of middle selection _Vec[j], more new route present position S:=S ∪ min_S _Vec[j], p=p+1 also changes step S3.2.3 over to;

Step S3.2.10: if S _p(x _i, y _i) be the climbing node in the potential energy surface, then with S _p(x _i, y _i) be that the boundary is with neighbor node ordered queue min_S _VecTwo minutes, free energy numerical value was greater than S _p(x _i, y _i) the node composition sequence be L_S _VecLess than S _p(x _i, y _i) the node composition sequence be S_S _VecWherein, the random number that produces according to program decides next step moving direction of ligand molecular, random number adopts [0,1] random function to produce, preferably, be set in the random number that obtains among a small circle the time, the response path direction namely changes step S3.2.11 over to, at S_S still according to the minimum energy path direction _VecSelect in the sequence; Otherwise, produce the at random effect of perturbation of response path, namely change step S3.2.12 over to, at L_S _VecIn select; More preferably, described be set as among a small circle [0.3,1];

Step S3.2.11: at S _pThe ordered queue min_S of neighbor node _VecOr S_S _VecThe middle position min_S that selects minimum _Vec[j] or S_S _Vec[j] is as next step the Direction of Reaction of ligand molecular, with step S3.2.7; Preferably, be absorbed in endless loop because pacing up and down in this position for avoiding the track search process, according to user's tolerance the number of times that this position occurs in response path limited, more preferably, the number of times of described appearance is made as 8; If this position is not satisfied and is above-mentionedly caused that the definition of endless loop then changes step S3.2.12 over to, otherwise with the reaction site of this position as next step, 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 select at random 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 changes step S3.2.3 over to;

Step S3.2.13: the ligand molecular that iteration is obtained revert to the coordinate information of the response path S of target receptor and does not carry out the BFEL that coarseness is divided, and obtains ligand molecular and target receptor response path and in conjunction with the funtcional relationship f (S) of free energy.

10. method according to claim 1, wherein,

Described step S4 is characterised in that, according to the ligand molecular that obtains and the response path of target receptor, determines avtive spot f (S) minimum position in the course of reaction, transition state site f (S) maximum position and reaction terminating site, wherein,

According to f (S) minimum position, transition state site and reaction terminating site in conjunction with free energy data Δ G _Bind,, Δ G _Trans, △ G _Unbind, obtain in conjunction with free energy in conjunction with transition state theory:

Association reaction energy of activation: ${\mathrm{\ΔG}}_{\mathrm{on}}^{\≠}={\mathrm{\ΔG}}_{\mathrm{trans},}-{\mathrm{\ΔG}}_{\mathrm{unbind}}$ With dissociation reaction energy of activation: ${\mathrm{\ΔG}}_{\mathrm{off}}^{\≠}={\mathrm{\ΔG}}_{\mathrm{trans},}-{\mathrm{\ΔG}}_{\mathrm{bind}},$ And then according to Aileen's formula, following formula (6)-(8) are calculated and are obtained binding constant K _D, association rate constant k _OnWith the rate constants k of dissociating _OffFurther obtain the action time of ligand molecular and target receptor: 1/k _Off

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

ΔG _binding＝RTlnK _D (7)

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

Wherein, K _BBe Boltzmann constant, R is gas law constant, and T is absolute temperature, and h is Planck constant.

11. method according to claim 1, wherein,

Described step S5 is characterised in that, revert to each present position S that the free energy potential energy surface is obtained ligand molecular on this path by the response path that obtains _p(x _p, y _p) reaction coordinate (x _p, y _p), find the ligand molecular that is projected in this position and the binding pattern P of target receptor by reaction coordinate _i(x _i, y _i), at this x _p=x _i, y _p=y _i, obtain thus ligand molecular at active binding site, the transition state site and binding pattern target receptor, learn avtive spot and transition state site on the corresponding target receptor.

12. a computer based application system, it comprises:

The first module, it carries out step S1 as claimed in claim 1;

The second module, it carries out step S2 as claimed in claim 1;

The 3rd module, it carries out step S3 as claimed in claim 1;

Four module, it carries out step S4 as claimed in claim 1;

The 5th module, it carries out step S5 as claimed in claim 1.

Images