JP2621842B2 - Method for searching for stable structure of biopolymer-ligand molecule - Google Patents

Description

Description: TECHNICAL FIELD The present invention relates to a method for searching for a stable complex structure of a biopolymer-ligand molecule, which can be used for structural design of medicines, agricultural chemicals, and other physiologically active compounds.

BACKGROUND ART In order for drugs and biologically active substances to exhibit physiological activity, target biopolymers (pharmacological receptors involved in signal transduction between cells, as well as enzymes, cytokine proteins and other proteins, and It is necessary to bind strongly to complexes containing these as main components and nucleic acids). Since 1960, the three-dimensional structures of many biopolymers have been revealed at the atomic level by X-ray crystallography, and the results have been published in protein data banks.

These low molecular weight compound molecules (drug molecules, enzyme substrates / inhibitors / coenzyme molecules, etc.) that bind to biopolymers are called ligand molecules. The ligand molecule also includes a drug molecule candidate, for example, a drug molecule to be synthesized.

In order for such a stable bond between the biopolymer and the ligand molecule, the shape of the groove or dent (called a ligand binding site) that binds the ligand in the three-dimensional structure of the biopolymer and the shape of the surface of the ligand molecule are key holes. It is necessary that they are complementary like a key, and that an interaction that brings about a special affinity between both molecules works. It is known from the crystal analysis of a complex of a biopolymer and a ligand molecule that hydrogen bonding, electrostatic interaction, hydrophobic interaction, and the like are particularly important among these intermolecular interactions. Whether any given compound molecule can form a stable complex with the target biopolymer, and if so, what kind of binding mode (which functional group of the ligand binding site of the biopolymer corresponds to which functional group of the ligand molecule) It is extremely important in drug design and structure-activity relationship to know what kind of complex it is and how it interacts with the group). Exploring a stable binding mode means that, in the case of a ligand molecule having a degree of freedom in conformational change, it simultaneously forms a molecular conformation (active conformation) when bound to a target biopolymer, This is also important for drug design.

When taking the most stable structure in terms of free energy as a complex, the structure of the ligand molecule is determined by the crystal structure of the ligand molecule alone, the structure in solution, and the energy re-stability (or minimum) obtained by various energy calculations. It is known that the structure is often different.

It is impossible to determine the binding mode of all interesting ligand molecules to a biopolymer by an experimental method such as X-ray crystallography. The first reason is that it requires a great deal of effort and time to analyze the binding mode of each ligand molecule to a biopolymer by an experimental method, but the ligand molecule is an enzyme substrate. This is because, in such a case, the stable complex structure may not be detected in the experiment due to the progress of the oxygen reaction, and it may be difficult to obtain the compound as the ligand molecule even if it actually exists. Furthermore, it is sometimes required to predict the binding mode of the imaginary ligand molecule compound to the biopolymer and the stability of the complex, and to judge whether or not the synthesis is actually worthwhile. Such a prediction is extremely useful especially when a new drug or active substance is devised by molecular design.

Research and work to simulate a stable binding mode of a ligand molecule having an arbitrary structure to a biopolymer having a known three-dimensional structure is called a docking study. Conventionally, docking studies have been performed using molecular models.However, there are problems such as the time and labor required to assemble the molecular models, low model accuracy and low reproducibility of results, and lack of quantitativeness. Was. As a method for solving such a problem, a simulation method using a computer and computer graphics has recently been generally used.

In a simulation method using a computer and computer graphics, a rough initial complex structure is set interactively on the computer graphics screen by visual judgment, and then refined and quantified computationally. This is most commonly done.
However, in this method, it is easy for an operator to prejudice the binding mode and the molecular conformation, it is extremely difficult to reach a correct solution from among the enormous possibilities, and it takes time and effort. In addition, the results differ depending on the operator, and the results lack reliability and reproducibility.

In order to solve such problems, a few studies have been made on a method that does not require the operator's subjectivity to intervene. Kuntz and colleagues describe the approximate representation of a ligand molecule as a few spheres, and instead of changing the degree of freedom of the molecular conformation, change the relative relationship of the spheres to fit the ligand molecule binding site of the receptor. (Docking Flexible Ligands to Macromolecular R
eceptors by Molecular Shape, RLDes Julais, RPShe
ridan, J. Scott Dixon, IDKuntz, and R. Venkataraghav
an: J. Med Chem. (1986) 29 , 2149-2153). However, since this method cannot deal with specific intermolecular interactions including hydrogen engagement, it is not effective in reliably determining the binding mode or conformation that forms a stable complex.

In addition, F. Jiang et al. Devised an automatic docking method that qualitatively considers intermolecular interactions such as hydrogen bonding and electrostatic interaction in addition to complementation of the shape noted by Kuntz et al. (F. Ji
ang, S. Kim, J. Mol. Biol. 219 , 79 (1991)). However, this method does not take into account the degree of freedom of conformation of the ligand molecule, and is insufficient in quantitative evaluation of intermolecular interaction.

Other methods have been studied, but none of them consider the specific interaction between the molecules or the degree of freedom of the conformation of the ligand molecule. However, it has a drawback that the docking state can be searched only by the structure, etc., and is not an effective method. In general, the binding mode of the ligand molecule to the biopolymer and the active conformation of the ligand molecule are completely coupled, and it is meaningless without considering both. The binding mode and active conformation can each have a vast number of possibilities, and without considering all of the combinations, it is not possible to cover all possible complex structures, but all of the binding modes and active conformations A method for searching for a complex structure in consideration of the combination of the two has not yet been established.

Accordingly, an object of the present invention is to provide a method for searching for a stable complex structure of a biopolymer-ligand molecule, which solves the above problems.

Another object of the present invention is to provide a method capable of simultaneously searching for a binding mode of a ligand molecule to a biopolymer and an active conformation of the ligand molecule.

DISCLOSURE OF THE INVENTION As a result of intensive efforts, the present inventors have developed a binding mode of a stable complex in consideration of hydrogen bonding, electrostatic interaction, and Van der Waals force as interactions between a biopolymer and a ligand molecule. By searching for the active conformation of the ligand molecule at the same time as the search, a method for automatically docking any ligand molecule to the ligand binding region of the biopolymer was developed, and the above-mentioned problem was successfully solved. That is, the present invention provides: (1) The association between a dummy atom set at a position of a hetero atom that can be a hydrogen bond partner of a hydrogen bond functional group in a biopolymer and a hydrogen bond hetero atom in a ligand molecule. A first step of covering the hydrogen bonding mode between the biopolymer and the ligand molecule by covering in combination; (2) comparing the distance between the dummy atoms and the distance between the hydrogen-bonding heteroatoms; By doing so, biopolymers-
A second step of simultaneously estimating the hydrogen bonding mode between the ligand molecules and the conformation of the hydrogen bonding portion of the ligand molecule; and (3) the hydrogen bonding mode and the conformation obtained in the second step are used in the ligand molecule. A third step of obtaining a complex structure of a biopolymer-ligand molecule by replacing the coordinates of all atoms of the ligand molecule with a coordinate system of a biopolymer based on the correspondence between the hydrogen-bonding heteroatom and the dummy atom. The present invention provides a method for searching for the structure of a stable biopolymer-ligand molecule complex containing the same.

Further, the present invention provides: (1) a method comprising: A first step of covering the hydrogen bonding mode between the biopolymer and the ligand molecule by covering the association in combination; (2) the distance between the dummy atoms and the distance between the hydrogen bonding heteroatoms; By comparing
A second step of simultaneously estimating the hydrogen bonding mode between the ligand molecules and the conformation of the partial structure of the ligand molecule; (3) the part of the ligand molecule based on the hydrogen bonding mode and the conformation obtained in the second step A third step of preserving a combination of a dummy atom and a hydrogen-bonding heteroatom that gives a hydrogen bonding mode impossible in the structure, and a hydrogen-bonding heteroatom that cannot form a hydrogen bond with the dummy atom; (4) third step Except for the combination containing the hydrogen-bonding heteroatom conserved in the step and the combination of the dummy atom and the hydrogen-bonding heteroatom conserved in the third step, the combination of the dummy atom and the hydrogen-bonding heteroatom in the ligand molecule A fourth step of covering the hydrogen bonding mode between the biopolymer and the ligand molecule by covering the association in combination; (5) the distance between the dummy atoms and the hydrogen bonding heteroatom; By comparing the distance, biopolymers -
A fifth step of simultaneously estimating the hydrogen bonding mode between the ligand molecules and the conformation of the hydrogen bonding portion of the ligand molecule; and (6) for each hydrogen bonding mode and conformation obtained in the fifth step, A sixth step of obtaining a complex structure of a biopolymer-ligand molecule by replacing the coordinates of all atoms of the ligand molecule with a biopolymer coordinate system based on the correspondence between the hydrogen-bonding hetero atom and the dummy atom. It is intended to provide a method for searching for a structure of a stable complex of a biopolymer-ligand molecule. In this way,
The search for a stable structure of a biopolymer-ligand molecule can be speeded up, and the structure of a stable complex can be reduced in a short time even when the biopolymer and / or the ligand molecule have a complicated structure. It becomes possible to search.

Furthermore, the present invention relates to (1) associating a dummy atom set at a position of a heteroatom which can be a hydrogen bond partner of a hydrogen bondable functional group of a biopolymer with a hydrogen bondable hetero atom in a ligand molecule. A first step of covering the hydrogen bonding mode between the biopolymer and the ligand molecule by covering in combination, (2) comparing the distance between the dummy atoms and the distance between the hydrogen bonding heteroatoms A second step of simultaneously estimating the hydrogen bonding mode between the biopolymer and the ligand molecule and the conformation of the hydrogen bonding portion of the ligand molecule, (3) maintaining the hydrogen bonding mode obtained in the second step. While
Optimize the conformation of the ligand molecule so that the position of the hydrogen-bonding heteroatom in the dummy atom and the ligand molecule coincide,
Next, a third step of removing the conformation of the ligand molecule having a high energy in the molecule, (4) for each conformation not removed in the third step, the difference between the hydrogen-bonding hetero atom and the dummy atom in the ligand molecule A fourth step of obtaining the complex structure of the biopolymer-ligand molecule by replacing the coordinates of all atoms of the ligand molecule with the coordinate system of the biopolymer based on the relative relationship; (5) the complex obtained in the fourth step From the body structure, the complex structure in which the intramolecular energy of the hydrogen bonding portion in the ligand molecule and the intermolecular interaction energy of the biopolymer-hydrogen bonding portion in the ligand molecule are high is removed, and then the remaining complex is removed. A fifth step of optimizing the structure of the body structure, (6) for each complex structure obtained in the fifth step, a new complex structure is generated by generating a conformation of a non-hydrogen bonding portion of the ligand molecule. 6th to get And (7) removing, from the complex structure obtained in the sixth step, a complex structure having a high intramolecular energy of the entire ligand molecule and a high intermolecular interaction energy of the biopolymer-ligand molecule, and then removing the remaining complex structure. It is intended to provide a method for searching for the structure of a biopolymer-ligand molecule stable complex including a seventh step of optimizing the structure of the complex structure. By this method, an appropriate complex structure can be selected even if the number of conformations to be generated is reduced, and a stable complex with high accuracy can be obtained.

In the present invention, a biopolymer is a concept that includes not only a polymer found in a living body but also a molecule simulating a polymer found in a living body.

The hydrogen bonding functional group is a concept including a functional group and an atom that are considered to be involved in hydrogen bonding.

The term "hydrogen-bonding heteroatom" refers to a heteroatom constituting a hydrogen-bonding functional group present in a ligand molecule.

The term “hydrogen-bonding moiety” refers to a structural part of a ligand molecule that includes a hydrogen-bonding heteroatom associated with a dummy atom, and the non-hydrogen-bonding moiety refers to a structural part other than the hydrogen-bonding moiety. Shall be referred to.

BRIEF DESCRIPTION OF THE DRAWINGS FIG. 1 schematically shows the main steps in the method of the invention.

FIG. 2 shows a flow chart of the method of the present invention. In the figure, S indicates a step.

FIG. 3 shows the molecular structures of methotrexate (MTX) and dihydrofolate (DHF). In the figure, the circled atoms are hydrogen-bonding heteroatoms, and a to d indicate bonds that are systematically rotated.

FIG. 4 shows the structure of the most stable model of the dihydrofolate reductase-methotrexate complex obtained by the method of the present invention. In the figure, the thick line shows the molecular structure of methotrexate, and the thin line shows a part of the molecular structure of dihydrofolate reductase.

FIG. 5 shows the conformation (a) of methotrexate in the crystal structure used alone as the input structure, and the methotrexate in the most stable model of dihydrofolate reductase-methotrexate complex obtained by the method of the present invention. Conformation (b) and conformation of methotrexate in the crystal structure of dihydrofolate reductase-methotrexate complex (c)
Is shown.

FIG. 6 shows the stereospecificity of the dihydrofolate reductase reaction.

FIG. 7 shows the structure of the most stable model of the dihydrofolate reductase-dihydrofolate complex obtained by the method of the present invention. In the figure, the bold line shows the molecular structure of dihydrofolate, and the dotted line shows a part of the molecular structure of dihydrofolate reductase.

FIG. 8 is a schematic diagram (a) of the binding mode of methotrexate to dihydrofolate reductase in the crystal structure of the dihydrofolate reductase-methotrexate complex, and is estimated from the stereospecificity of the dihydrofolate reductase reaction product. (B) shows a schematic diagram of the mode of binding of dihydrofolate to dihydrofolate reductase.

BEST MODE FOR CARRYING OUT THE INVENTION A preferred embodiment of the present invention will be described with reference to the flowchart of FIG. In FIG. 2,
S indicates each step.

First, the number and atomic coordinates of the atoms constituting the biopolymer (excluding the atomic coordinates of hydrogen atoms) are input (S
1).

The atomic coordinates of a biopolymer can be calculated from the three-dimensional structure obtained by X-ray crystallographic analysis or NMR analysis, obtained using a protein crystal database, etc. The atomic coordinates of the polymer model can be used. This atomic coordinate is preferably represented by a three-dimensional coordinate system. In addition, a cofactor that plays an important structural and functional role by binding to a biopolymer is regarded as a part of the biopolymer, and the atoms of the biopolymer in a state where the cofactor is bound are considered. The subsequent steps may be performed by inputting the coordinates. Examples of the cofactor include a coenzyme, a water molecule, and an iron ion.

Next, the atomic coordinates of the hydrogen atoms constituting the amino acid residues present in the biopolymer are calculated (S2).

In general, the positions of hydrogen atoms in biopolymers cannot be determined by experimental methods such as X-ray crystallography and NMR analysis, and information on hydrogen atoms cannot be obtained from protein crystal databases or the like. Based on the structure of the amino acid residue existing in the above, the atomic coordinates of the hydrogen atom constituting the amino acid residue are calculated. A hydrogen atom whose atomic coordinates cannot be unambiguously determined because it is bonded to a rotatable amino acid residue is preferably calculated on the assumption that it is present in the trans position.

Next, a charge is added to the atoms constituting the amino acid residues present in the biopolymer (S3).

As the value of this charge, a literature value calculated for each amino acid, for example, a Weiner value (Weiner, SJ, Koll
man, PACase, DA, Singh, UC, Ghio, C., Alagona, G., Pr
ofeta, S., Jr. and Weiner, P., J. Am. Chem. Soc., 106 (198
4) 765-784) can be used.

Next, a hydrogen bond number is added to the heteroatom constituting the hydrogen bondable functional group present in the biopolymer (S4).

Hydrogen bond numbers can be added according to Table 1 below.

Next, a ligand binding region is designated (S5).

As the above-mentioned region, a region including an arbitrary site of the biopolymer can be selected. Depending on the purpose, a region containing a ligand binding pocket, a site of a biopolymer surrounding the site, and, if necessary, a site of a biopolymer to which another molecule such as an effector is bound may be, for example, as shown in FIG. As shown in (a), it can be designated as a rectangular parallelepiped. The ligand binding pocket refers to an inner hole on a concave molecular surface of a biopolymer, to which a ligand molecule such as a substrate or an inhibitor binds.

To specify the range of the ligand binding region, use the program GREE
N (Journal of Computer Ages Molecular Design, vol.
1, pp. 197-210 (1987), Nobuo Tomioka, Akiko Itai, and
You can use some of the functions of Yoichi Iitaka).

Three-dimensional grid points are generated in the area specified in S5, and numbering and grid point information are calculated for each three-dimensional grid point (S6).

The three-dimensional lattice points are points on the three-dimensional lattice generated at regular intervals in the ligand binding region of the biopolymer. Lattice point information refers to van der Waals and electrostatic interaction energies acting between the biomolecule and the probe atoms, calculated assuming that there is a probe atom on each three-dimensional lattice point, and hydrogen bonding properties. This is a concept that includes local physicochemical information of the ligand binding region.

By using the grid point information of the three-dimensional grid points, it becomes possible to speed up the approximate calculation of the interaction energy between the biopolymer and the ligand molecule performed in the subsequent steps, and to set in the subsequent steps. The position of the dummy atom can be rationally determined. As a result, a docking model of a biopolymer and a ligand molecule can be comprehensively searched in a short time.

The three-dimensional lattice points may be generated at a fixed interval of 0.3 to 2.0 angstroms, preferably 0.3 to 0.5 angstroms, in the area designated by S5.

As the probe atom, it is preferable to adopt all the types of atoms contained in the compound which may be a ligand.

The van der Waals interaction energy acting between each probe atom arranged at each three-dimensional lattice point and the biopolymer can be calculated by an ordinary atom pair calculation using an experimental potential function. As an experimental potential function, Lennard-Jon as expressed by the following equation
The es type function can be used.

(In the formula, i represents a number indicating the position of the probe atom, j represents the atomic number of the biopolymer, and A and B represent parameters for determining the position and magnitude of the minimal potential.
_ij is the i-th probe atom and j of the biopolymer
Represents the distance between the atoms. The van der Waals interaction energy between the ligand molecule and the biopolymer when the ligand molecule is placed on this three-dimensional lattice point can be calculated from the following equation.

(Where m is also the number of atoms in the ligand molecule, N is the number of atoms in the ligand molecule, and m ₀ is the number of the three-dimensional lattice point closest to the mth atom.) As parameters of A and B, values of Weiner et al. (Weiner, SJ, Kollman, PA, Case, DA, Singh, UC,
Ghio, C., Alagona, G., Profeta, S., Jr.and Weiner, P., JA
m. Chem. Soc., 106 (1984) 765-784) can be used.

The electrostatic interaction energy acting between each probe atom arranged at each three-dimensional lattice point and the biopolymer can be calculated from the following equation.

( _Where i, j, and r _ij are as defined above, q _j is the charge of the jth atom of the biopolymer, K is a constant for converting energy units, and ε is A constant value may be used as the above-mentioned dielectric constant.
rshel et al. suggested r _ij dependent values (Warshe
1, A., J. Phys. Chem., 83 (1979) 1640-1652).

The electrostatic interaction energy between the ligand molecule and the biopolymer when the ligand molecule is placed on this three-dimensional lattice point can be calculated from the following equation.

(Wherein, m, N, and m ₀ are as defined above.) Hydrogen bonding refers to the arrangement of either a hydrogen-donating atom or a hydrogen-accepting atom at its three-dimensional lattice point. If it is, it can form a hydrogen bond with the hydrogen bonding functional group of the biopolymer, or it can not form a hydrogen bond with the hydrogen bonding functional group of the biopolymer regardless of which atom is arranged Information about.

The hydrogen bonding property of a three-dimensional lattice point can be determined as follows. A distance DP between a certain three-dimensional lattice point P and a hydrogen donor atom D present in the biopolymer is a distance (for example, 2.5 to 3.1 angstroms) at which a hydrogen bond can be formed; Angle ∠DH formed by hydrogen H and D
If P is an angle capable of forming a hydrogen bond (for example, 30 ° or more), the three-dimensional lattice point may be considered to be hydrogen-accepting. Similarly, a distance DA between a certain three-dimensional lattice point P and a certain hydrogen acceptor atom A present in a biopolymer is a distance capable of forming a hydrogen bond, and P, a lone electron pair L, If the angle ∠ALP formed by A and A is an angle at which a hydrogen bond can be formed, this three-dimensional lattice point may be regarded as hydrogen-donating. If a certain three-dimensional lattice point is neither hydrogen-accepting nor hydrogen-donating, it may be considered that there is no hydrogen bonding. By adding 1 to the three-dimensional lattice points that are hydrogen-accepting, 2 to the three-dimensional lattice points that are hydrogen-donating, and 0 to the three-dimensional lattice points that do not have hydrogen bonding properties, the hydrogen bonding properties of the three-dimensional lattice points are added. Can be represented.

A hydrogen bonding functional group present in a biopolymer,
From those existing in the region specified in S5, those that are expected to form hydrogen bonds with the ligand molecule are selected (S7).

When a plurality of hydrogen-bonding functional groups are present, selection is made according to their importance.

Based on the grid point information of the three-dimensional grid points calculated in S6,
A dummy atom is set for each hydrogen-bonding functional group selected in S7 (S8).

In this step, first, based on the hydrogen bonding properties of the three-dimensional lattice points calculated in S6, a region capable of forming a hydrogen bond with the highly hydrogen bonding functional group selected in S7 (hereinafter referred to as a “hydrogen bonding region”) ), And then construct a hydrogen bonding region in which there is an appropriate number, for example, 5 to 20 three-dimensional lattice points, within the hydrogen bonding region and outside the van der Waals radius of the other atoms. This can be performed by arranging a dummy atom at the center of the three-dimensional lattice point. The hydrogen bonding region is a region having the same hydrogen bonding property and configured from a group of three-dimensional lattice points adjacent to each other.

The dummy atom is given the same hydrogen bonding property as the three-dimensional lattice point where the dummy atom is arranged at the center.

Two or more dummy atoms are set from one hydrogen bonding functional group, and one dummy atom is not set.

Next, for the atoms constituting the ligand molecule,
The atom name, atom type number, atomic coordinates, charge, and if the atom is a hydrogen-bonding heteroatom, enter the hydrogen bond number (S9).

The atomic type numbers of the atoms that make up the ligand molecule are We
iner et al. (Weiner, SJ, Kollman, PA, Cese, DA, Singh,
UC, Ghio, C., Alagona, G., Profeta, S., Jr. andg Weiner,
P., J. Am. Chem. Soc., 106 (1984) 765-7840).

The atomic coordinates of the ligand molecule can be obtained by analyzing the structure of a single crystal,
It can be calculated from a three-dimensional structure or the like obtained from a model database based on a crystal database or energy calculation. The three-dimensional structure represented by the atomic coordinates of the ligand molecule is
Basically, the conformation does not matter at all, as long as the geometric quantities such as the bond distance between atoms and bond angles are accurate. However, in order to calculate the charge of each atom constituting the ligand molecule as accurately as possible, it is preferable to input the atomic coordinates that give a stable three-dimensional structure independently to some extent.

The charge of each atom constituting the ligand molecule can be calculated by molecular orbital calculation using the MNDO method or the AM1 method in the MOPAC program.

The hydrogen bond number of the hydrogen bondable hetero atom present in the ligand molecule can be added according to Table 1 described above.

By changing the type of the ligand molecule and repeating steps S9 to S31, it is possible to search for the binding mode of various ligand molecules to the same biopolymer.

The minimum value 1 _min and the maximum value 1 _max of the number of hydrogen bonds formed between the biopolymer and the ligand molecule are specified (S10).

Although any value may be set as these values, it is preferable to specify values before and after the value predicted from the number of hydrogen bonds formed between the biopolymer and the known ligand molecule.

For the ligand molecule, a bond capable of torsional rotation and a mode of the rotation are specified (S11).

10 ゜ -120 if the torsionally rotatable connection is a single connection
It is preferable to designate a rotation mode for systematically rotating at a constant rotation angle of ゜. If there is a torsionally rotatable bond in the ring structure, it is preferable to specify a method of sequentially inputting various possible ring structures from a file.

A correspondence is established between the dummy atom and the hydrogen-bonding heteroatom in the ligand molecule (S12).

Assuming that the number of dummy atoms is m and the number of hydrogen-bonding heteroatoms present in the ligand molecule is n, the number N (i) of correspondences (combinations) forming i hydrogen bonds is _m P _i
× _n C _i (where P represents a sequence and C represents a combination). It is preferred to select all possible combinations of dummy atoms and hydrogen bonding heteroatoms in the ligand molecule.

For all i satisfying the relationship 1 _min ≤ i ≤ 1 _max , S1
Repeat steps 2 to S30. As a result, the difference between the dummy atom and the hydrogen-bonding heteroatom existing in the ligand molecule Are all selected. As a result, all combinations of hydrogen bonds formed between the biopolymer and the ligand molecule are selected, and the mode of binding between the biopolymer and the ligand molecule can be systematically and efficiently searched.

One of the combinations of the dummy atom and the hydrogen-bonding heteroatom in the ligand molecule, which has been associated in S12, is selected (S13).

At this time, a combination in which the hydrogen bonding property of the dummy atom and the hydrogen bonding property of the hydrogen bonding hetero atom in the ligand molecule do not match is not selected.

For all combinations that have been matched in S12,
Steps S30 to S30 are repeated.

Next, for the dummy atoms included in the combination selected in S13, the distance between each dummy atom is calculated (S14).

If both the number of dummy atoms and the number of hydrogen-bonding heteroatoms in the ligand molecule are one, skip steps S14 to S21, skip to step S22, and then skip steps S23 to S25. Then, jump to step S26. When one of the number of dummy atoms or the number of hydrogen-bonding heteroatoms present in the ligand molecule is one, S140 to
Skip to step S22 without performing step S21.

Next, the ligand molecule is divided into a hydrogen bonding part and a non-hydrogen bonding part (S15).

For example, as shown in FIG. 1 (b), the ligand molecule is divided into a hydrogen bonding part and a non-hydrogen bonding part.

For the hydrogen bonding portion of the ligand molecule divided in S15, a bond capable of torsional rotation and a mode of the rotation are selected (S16).

The bonds selected in S16 are rotated according to the rotation mode specified in S11, thereby sequentially generating the conformation of the ligand molecules (S17).

The following steps S18 to S30 are repeated for each generated conformation.

For the conformation generated in S17, the interatomic distance of each hydrogen-bonding heteroatom in the ligand molecule included in the combination selected in S13 is calculated (S18).

The value of F, which is the sum of the square of the difference between the interatomic distance of the dummy atom obtained in S14 and the interatomic distance of the hydrogen bonding heteroatom in the corresponding ligand molecule obtained in S18, is constant. The conformation of the ligand molecule exceeding the range is removed (S19).

By this step, it is possible to efficiently cover the hydrogen bonding mode between the biopolymer and the ligand molecule and the possibility of conformation of the ligand molecule.

Where k _a ≦ F ≦ k _a ′ (where n is the number of hydrogen bonds, r _di is the distance between the i-th dummy atoms, and r _li is the i-th hydrogen-bonding heteroatom in the ligand molecule. , K _a is the lower limit of F, k _a ′
Represents the upper limit of F. ) Remove conformations other than those that are. Examples of the above-mentioned k _a,
0.6 to 0.84 and k _a ′ are preferably 1.2 to 1.4.

With respect to the conformation of the ligand molecule not removed in S19, the conformation of the hydrogen bonding portion of the ligand molecule is optimized so that the value of F is minimized (S20).

This step modifies the conformation of the ligand molecule so that the biopolymer and the ligand molecule form stable hydrogen bonds.

The value of the above function F is calculated by the Fletcher-Powell method (R. Fletcher, MJ DPowell, Comp.
uter J., 6 , 163. (1968)), it is possible to optimize the conformation of the hydrogen bonding portion of the ligand molecule.

Next, the intramolecular energy of the hydrogen bonding portion of the ligand molecule optimized in step S20 is calculated, and the conformation in which the value of the intramolecular energy is equal to or more than a certain value is removed (S2
1).

For example, when the intramolecular energy of the hydrogen bonding portion of the ligand molecule is calculated using the molecular force field of AMBER 4.0, the conformation in which the intramolecular energy becomes 100 Kcal / mol or more may be removed.

The atomic coordinates of the ligand molecule are changed so that the coordinates of the hydrogen-bonding heteroatom of the ligand molecule having the conformation obtained in the step up to S21 and the corresponding dummy atom are identical. Replace with the coordinate system (S22).

To do that, Kabsh's least squares method (W. Kabsh, Acta Crys
t., A321 , 922 (19764), W. Kabsh, Acta Cryst., A34 , 827
(1978)) can be used. In this way, possible hydrogen bonding modes and conformations of the hydrogen bonding portion of the ligand molecule can be roughly estimated at the same time.

Next, the intermolecular interaction energy (sum of van der Waals interaction energy and electrostatic interaction energy) between the biopolymer and the hydrogen bonding portion of the ligand molecule,
Then, the intramolecular energy of the hydrogen bonding portion of the ligand molecule is calculated (S23).

The intermolecular interaction energy E _inter between the biopolymer and the hydrogen bonding portion of the ligand molecule is represented by the van der Waals interaction energy G _vdW (lattice point information of the three-dimensional lattice point closest to each atom k in the ligand molecule). Using k) and the electrostatic interaction energy G _elc (k), it can be calculated from the following equation.

(In the formula, q _k represents the charge of the atom k.) Further, the intramolecular energy E _intra of the hydrogen bonding portion of the ligand molecule can be calculated using a known force field published by a known method. it can. For example, E _intra is AMBER
It can be calculated by the following equation using a known known force field such as 4.0.

(Where V _n is a constant given to the arrangement of the atom types of the four atoms constituting the torsion angle, ng is a constant representing the symmetry of the potential with respect to the torsion angle rotation, and Φ is the torsion angle. , Γ represent the phase of the potential with respect to the torsion angle rotation (given for the arrangement of the atom types of the four atoms constituting the torsion angle), and A _ij and B _ij represent the i-th and j-th positions in the ligand molecule. The constant set for the atom type pair of atoms, R _ij is the distance between the i-th and j-th atom, q _i is the charge of the i-th atom in the ligand molecule, q _i
Represents the charge of the j-th atom in the ligand molecule, and ε represents the dielectric constant. The conformation of the ligand molecule in which the sum of the intermolecular interaction energy and the intramolecular energy calculated in S23 has reached a certain value or more is removed (S24).

For example, it is preferable to remove a conformation in which the sum of the above energies is 1500 kcal / mol or more.

The structure of the hydrogen bonding portion of the ligand molecule having a conformation not removed by S24 is optimized (S25).

By optimizing the torsion of the hydrogen bonding portion of the ligand molecule and the position and orientation of the ligand molecule by the structure optimization calculation, the structure of the hydrogen bonding portion of the ligand molecule can be optimized. For example, to minimize the _total energy E _total calculated from the following equation,
The structure of the hydrogen bonding part of the ligand molecule can be optimized by the mplex method (HITACH. Library program, SOFSM) or the like.

EB _total = E _inter + E _intra + Whb · Nhb · Chb (where E _inters and E _intra are as described above, Whb is the weight, Nhb is the number of hydrogen bonds, and Chb is the stabilization energy by one hydrogen bond. (For example, 2.5 kcal / mol). By rotating the bond capable of torsional rotation of the non-hydrogen bonding part of the ligand molecule according to the rotation mode specified in S11, the conformation of the ligand molecule is sequentially generated. (S26).

The following steps S27 to S30 are repeated for each generated constellation.

Next, the intermolecular interaction energy between the biopolymer and the ligand molecule and the intramolecular energy of the ligand molecule are calculated (S27).

The above-mentioned intermolecular interaction energy and intramolecular energy can be calculated in the same manner as in S23. However, in S23,
Although the calculation is performed only for the hydrogen bonding portion of the ligand molecule, in S27, the calculation is performed for the entire ligand molecule.

The conformation of the ligand molecule in which the sum of the intermolecular interaction, energy, and intramolecular energy calculated in S27 is equal to or more than a certain value is removed (S28).

The upper limit of the sum of the above energies is 1500 kcal / m
ol is preferred.

The structure of the entire ligand molecule having the conformation not removed in S28 is optimized (S29).

This step results in a stabilized complex of the biopolymer and the ligand molecule and an active conformation of the ligand molecule.

In the same manner as in the step of S25, the structure of the ligand molecule can be optimized. However, in S25, only the structure of the hydrogen bonding portion of the ligand molecule is optimized, but in S29
Then, the structure is optimized for the entire ligand molecule.
The structure of the complex obtained in the steps up to S29 is referred to as an initial complex model.

The structure of the biopolymer-ligand molecule is optimized by changing the conformation of the biopolymer (S30).

Until the step of S29, the biopolymer has been treated as a rigid one, but in S30, the conformation of the biopolymer is changed to minimize the energy of the initial complex model. By minimizing this energy using the publicly known energy minimization calculation program AMBER, the structure of the initial complex model can be optimized. In this step, it is preferable that all the water molecules and other molecules found in the crystal structure of the complex are arranged in the complex structure to optimize the complex structure.
The structure of the complex obtained by the step of S30 is referred to as a final complex model.

The energy of the biopolymer-ligand molecule complex structure obtained in the steps up to S30 is calculated, and the complex structures are ranked in ascending order of the energy (S31).

Of the structures of the final complex model, the structure giving the minimum energy often matches well with the complex structure already obtained experimentally. Therefore, it is reasonable to consider the structure with the lowest energy as the correct composite structure. However, the order of the magnitude of the complex energy varies slightly depending on the precision of the force field used, the precision of the atomic coordinates of the biopolymer, and other factors. It is preferable to consider

Also, in a more preferred embodiment of the invention, where it is expected that a very large number of dummy atoms will be created from the biopolymer, and if the ligand molecule has a considerable degree of conformational freedom or In the case of having a hetero atom of, after inputting information about the ligand molecule in the step of S9, specify the partial structure of the ligand molecule, perform the steps from S10 to S24 for the specified partial structure, the partial structure Combination of a dummy atom and a hydrogen-bonding heteroatom in a ligand molecule that gives an impossible hydrogen-bonding mode in a molecule, and information on a hydrogen-bonding heteroatom in a ligand molecule that cannot form a hydrogen bond with a dummy atom .
The partial structure of the ligand molecule can be arbitrarily determined without any structural restriction, but a structural portion containing three or more hydrogen-bonding functional groups is preferable, and the degree of conformational freedom of the structural portion is not limited. It doesn't matter. Next, steps S10 to S31 are performed on the entire structure of the ligand molecule, and the combination including the conserved hydrogen bonding heteroatom and the combination of the conserved dummy atom and the hydrogen bonding heteroatom are S12. Is removed from the correspondence combination created in the step. As a result, the number of combinations and conformations to be tested is reduced, and the time required for searching for the structure of a stable biopolymer-ligand molecule complex is greatly reduced. Such a method is referred to as a Pre-Pruning method (hereinafter, referred to as a “PP method”).
Shall be called.

The PP method is based on the presumption that hydrogen bonding modes and ligand conformations that are impossible in the partial structure of the ligand molecule are also impossible in the overall structure of the ligand molecule. The PP method has made it possible to significantly reduce the search time without affecting the accuracy and reliability of the resulting biopolymer-ligand complex structure at all.

[Example]

In order to demonstrate the effectiveness of the method of the present invention, the method of the present invention was applied to a structure search for a complex of dihydrofolate reductase and its inhibitor or its substrate. Dihydrofolate reductase is one of the enzymes widely studied as targets of clinical drugs, and the crystal structures of complexes with various inhibitor molecules have been analyzed. For example, the crystal structure of the dihydrofolate reductase-inhibitor methotrexate complex has been experimentally analyzed, and is therefore often used worldwide as a model for methodology research for the analysis of stable complex structures. .

The atomic coordinates of dihydrofolate reductase (hereinafter referred to as “DHFR”) are as shown in the Protein Data Bank entry.
The values obtained from the atomic coordinates of the binary complex of DHFR-methotrexate of E. coli available from 4DFR were used. DHFR-
Based on the crystal structure of the folate-NADP + ternary complex,
The atomic coordinates of the NADP + molecule were added to the atomic coordinates of the binary complex. In addition, all water molecules were removed from the atomic coordinates of the binary complex except for the two water molecules that bind very strongly to DHFR and cannot be removed chemically. 2 not removed
One of the water molecules has a hydrogen bond to tryptophan at position 21 and aspartic acid at position 26 of DHFR, and the other has a leucine at position 114 and a threonine at position 116 of DHFR (see FIG. 8). thing).

The atomic coordinates of the hydrogen atoms that make up the amino acid residues that are present in DHFR are represented by P, one of the programs that make up GREEN.
Calculated by DBFIL. In addition, the results of ab initio calculations by Weiner et al. (Weiner, SJ, Kollman, PA, Case, DA, Singh,
UC, Ghio, C., Alagona, G., Profeta, S., Jr, and Weiner,
Based on P., J. Am. Chem. Soc., 106, (1984) 765-784), a charge was imparted to an atom constituting an amino acid residue present in DHFR. However, it was assumed that the carboxylic acid side chain of aspartic acid at position 27 was ionized to the carboxylate anion. A hydrogen bond number was added to the hetero atom constituting the hydrogen bondable functional group present in DHFR by PDBFIL.

Dihydrofolate reductase - in the crystal structure of methotrexate binary complex, the size rectangular area of 12.80 × 16.8 × 11.6Å ³ holes near inner containing the inhibitor is a molecule methotrexate specified as the ligand binding region .

Three-dimensional lattice points were generated at intervals of 0.4 ° in the above rectangular parallelepiped region, and lattice point information was calculated for each three-dimensional lattice point. Hydrogen, carbon, nitrogen, and oxygen atoms present in the methotrexate molecule were used as probe atoms.

In addition to selecting the amino acid side chains and the main chain hydrogen bonding functional groups of DHFR protruding into the rectangular parallelepiped region, a total of 10 hydrogen bonding functional groups were selected. Thirteen dummy atoms were set for the hydrogen bonding functional group.

As the ligand molecule, ぱ, inhibition example methotrexate (hereinafter referred to as “MTX”) and substrate dihydrofolate (hereinafter “D
HF ". ).

Example 1 Methotrexate Molecule The following treatment was performed except for the terminal acetic acid group of the MTX molecule.

For the atomic coordinates of MTX, see the Cambridge Crystal Graphics Database (Cambrzidge Crystallographic Da
tabose), the atomic coordinates of the single crystal structure available. The charge of each atom of MTX is
It was calculated using the MNDO method. In addition, 1 of the pteridine ring of MTX
Since it is known that the nitrogen at the position is susceptible to protonation, the charge was calculated as being protonated.
In the structural formula of FIG. 3, the heteroatoms encircled are designated as hydrogen-bonding heteroatoms, and hydrogen bond numbers are added.
For the bonds a to d with an arrow in the MTX molecular structure shown in FIG. 3, the torsion angle is rotated by 60 ° in the range of 0 to 360 ° for the bond a, and 0 to 180 ° for the bond b. Rotate by 60 ° in the range of ゜, set 0 ° or 180 ° for the coupling c, and 180 ° for the coupling d
The conformation of the MTX molecule was generated.

In this method, the number of hydrogen bond combinations N (i) when i hydrogen bonds are formed becomes an enormous number as shown below.

_{N (7) = 11 P 7} × 10 C 7 = 199,584,000 ways _{N (6) = 11 P 6} × 10 C 6 = 69,854,400 types _{N (5) = 11 P 5} × 10 C 5 = 13,970,880 types N (4) = ₁₁ P ₄ × ₁₀ C ₄ = 207,900 ways It is possible to dock the entire structure of the MTX molecule to DHFR at once (Iris4D / 25 workstation, MTX
Approximately 30 hours for processing from the input of the atomic coordinates of the molecule to the optimization of the structure of the complex), and after docking the partial structure of the MTX molecule from the pteridine ring to the benzene ring to DHFR (the required time is slightly more than 4 minutes) ), The time required to search for a stable complex structure in consideration of the entire structure of the MTX molecule based on the results (PP method) was reduced (from inputting the atomic coordinates of the MTX molecule to optimizing the complex structure) Processing time about 30 minutes). Also, when the stable complex structure is searched by docking the entire structure of the MTX molecule to DHFR at a stretch, the partial structure of the MTX molecule is docked to DHFR, and based on the result, the entire structure of the MTX molecule is considered and stable. When searching for complex structures,
The same stable composite structure was obtained.

From the obtained 11 initial complex models, 5 models were selected in ascending order of energy, and AM
Energy was minimized using the BER program to obtain the final composite model. The energy minimization of the AMBER program improved both the energy value and the number of hydrogen bonds in each model, and changed the ranking of some models. These facts suggest that the effects of local modification of the biopolymer structure are not negligible. Table 2 shows the number of hydrogen bonds and the total energy of the higher model among the finally obtained stable complex models. FIG. 4 shows a most stable complex model showing the minimum energy (−176.71 kcal / mol). In Fig. 4, the thick line is MTX.
The dotted line shows a part of the molecular structure of DHFR. The binding mode of the DHFR-MDX complex in this complex model well reproduced the binding mode of the DHFR-MTX complex in the crystal structure of the complex obtained by X-ray crystallography. The conformation of the MTX molecule in the most stable complex model (FIG. 5 (b)) was
The conformation of the molecule (FIG. 5 (a)) and the conformation of the MTX molecule in the crystal structure of the DHFRG-MTX complex are shown in FIG. 5 together with the conformation of the MTX molecule (FIG. 5 (c). In contrast, the twist angle of C5-C6-C9-C10 of the MTX molecule is -139.1 ° in the crystal structure alone, whereas 25.3 ° in the most stable complex model, the crystal of the DHFR-MTX complex 28.4 構造 in the structure, MT in the most stable complex model
The conformation of the X molecule was significantly different from the initial structure and showed very good agreement with that in the crystal structure of the complex.

[Example 2] The binding mode between dihydrofolate molecule DHFR and its substrate DHF has not been elucidated by X-ray analysis, but the binding mode differs from MTX due to the stereospecificity of tetrahydrofolate, a product of oxygen reaction. Is estimated. That is, the hydrogen at the C6 position of tetrahydrofolate generated by the reducing action of DHFR is the coenzyme NADP
Although it is known to come from H as a hydride ion, the same binding mode as that of MTX to DHFR in the crystal structure of the ternary complex of DHER-MTX-NADPH analyzed by X-ray is DHFR-DHF Assuming the complex takes
Tetrahydrofolate having the opposite steric structure will be produced, and the stereospecificity of the DHFRG enzyme reaction cannot be explained (see FIG. 6).

Here, the structure of the stable DHFR-DHF complex was searched for in the same manner as in Example 1 without any such prejudice. However, as the atomic coordinates of DHF, the atomic coordinates created by performing atom substitution based on the atomic coordinates of the crystal structure of MTX alone were used. The charge of each atom of DHF was calculated by the same method as in Example 1. In addition, in the structural formula of FIG. 3, the heteroatoms surrounded by circles are designated as hydrogen-bonding heteroatoms, and hydrogen bonding numbers are added.

Table 2 shows the number of hydrogen bonds and the total energy of the higher-order stable composite model finally obtained. The most stable complex model showing the minimum energy (-173.46 Kcal / mol) is shown in FIG. In FIG. 7, the thick line shows the molecular structure of DHF, and the dotted line shows a part of the molecular structure of DHFR. The conformation of DHF in the most stable complex model was different from that of MTX in the crystal structure of the DHFR-MTX complex, and the face of the pteridine ring was inverted. DH in the most stable complex model
The binding mode of F to DHFR was in good agreement with the binding mode of DHF to DHFR (FIG. 8 (b)) estimated from the stereospecificity of the DHFR enzyme reaction product. This mode of binding of DHF to DHFR depends on the stereospecificity of the DHFR enzyme reaction (see FIG. 6).
Can be explained well. The conformation of DHF in the most stable complex model was very close to the conformation of folate in the crystal structure of the complex between DHFR and the substrate analog, folic acid.

These results demonstrate the effectiveness of the method of the present invention.

According to the method of the present invention, it has become possible to search for a stable complex structure of a biopolymer-ligand molecule in a short time.

Further, the method of the present invention has made it possible to reliably select a small number of complex structures including the most stable structure from all possible biopolymer-ligand molecule complex structures.

Furthermore, the method of the present invention enables an operator to search for a stable complex structure of a biopolymer-ligand molecule without subjective intervention.

Further, the method of the present invention makes it possible to search for a stable and highly reproducible stable structure of a biopolymer-ligand molecule.

Furthermore, by the method of the present invention, a search for a stable complex structure of a biopolymer-ligand molecule was conducted in which the degree of freedom of the conformation of the ligand molecule was considered and the quantitative evaluation of the intermolecular interaction was sufficiently performed. It is now possible to do.

Furthermore, the method of the present invention makes it possible to search for a stable complex structure even when the number of hydrogen bonds between the biopolymer and the ligand molecule is small, for example, when there is only one hydrogen bond.

INDUSTRIAL APPLICABILITY Using the method of the present invention, it is possible to predict what kind of binding mode a newly designed ligand molecule is likely to form with a target biopolymer and its stability. This allows efficient design and research of drugs and active compounds.

In addition, a novel lead compound can be searched for by applying the method of the present invention to a large number of ligand candidate compound groups for a certain target biopolymer.

Claims (3)

(57) [Claims] (1) Combining the association between a dummy atom set at the position of a heteroatom that can be a hydrogen bond partner of a hydrogen bondable functional group in a biopolymer and a hydrogen bondable heteroatom in a ligand molecule A first step of covering the hydrogen bonding mode between the biopolymer and the ligand molecule by comprehensively covering, (2) comparing the distance between the dummy atoms and the distance between the hydrogen-bonding hetero atoms By the biopolymer-
A second step of simultaneously estimating the hydrogen bonding mode between the ligand molecules and the conformation of the hydrogen bonding portion of the ligand molecule; and (3) for each hydrogen bonding mode and conformation obtained in the second step, A third step of obtaining a complex structure of a biopolymer-ligand molecule by replacing the coordinates of all atoms of the ligand molecule with a biopolymer coordinate system based on the correspondence between the hydrogen-bonding hetero atom and the dummy atom. A method for searching for a structure of a stable complex of a biopolymer and a ligand molecule. (1) Correspondence between a dummy atom set at a position of a heteroatom which can be a hydrogen bond partner of a hydrogen bondable functional group in a biopolymer and a hydrogen bondable heteroatom in a partial structure of a ligand molecule. A first step of covering the hydrogen bonding mode between the biopolymer and the ligand molecule by covering the attachments in combination; (2) reducing the distance between the hydrogen bonding heteroatoms of the distance between the dummy atoms; By comparison, the biopolymer-
A second step of simultaneously estimating the hydrogen bonding mode between the ligand molecules and the conformation of the partial structure of the ligand molecule; (3) the part of the ligand molecule based on the hydrogen bonding mode and the conformation obtained in the second step A third step of preserving a combination of a dummy atom and a hydrogen-bonding heteroatom that gives a hydrogen bonding mode impossible in the structure, and a hydrogen-bonding heteroatom that cannot form a hydrogen bond with the dummy atom; (4) third step Except for the combination containing the hydrogen-bonding heteroatom conserved in the step and the combination of the dummy atom and the hydrogen-bonding heteroatom conserved in the third step, the combination of the dummy atom and the hydrogen-bonding heteroatom in the ligand molecule A fourth step of covering the hydrogen bonding mode between the biopolymer and the ligand molecule by covering the association in combination; (5) the distance between the dummy atoms and the distance between the hydrogen bonding heteroatoms; The fifth step of simultaneously estimating the hydrogen bonding mode between the biopolymer and the ligand molecule and the conformation of the hydrogen bonding portion of the ligand molecule by comparing the distances, and (6) the hydrogen bond obtained in the fifth step By replacing the coordinates of all atoms of the ligand molecule with the biopolymer coordinate system based on the correspondence between the hydrogen-bonding heteroatom and the dummy atom in the ligand molecule for each mode and conformation, the biopolymer-ligand molecule A method for searching for a structure of a stable complex of a biopolymer-ligand molecule, comprising the sixth step of obtaining a complex structure of the above. 3. Combination of association between a dummy atom set at a position of a heteroatom which can be a hydrogen bond partner of a hydrogen bondable functional group in a biopolymer and a hydrogen bondable heteroatom in a ligand molecule. A first step of covering the hydrogen bonding mode between the biopolymer and the ligand molecule by comprehensively covering, (2) comparing the distance between the dummy atoms and the distance between the hydrogen-bonding hetero atoms By the biopolymer-
A second step of simultaneously estimating the hydrogen bonding mode between the ligand molecules and the conformation of the hydrogen bonding portion of the ligand molecule, (3) while maintaining the hydrogen bonding mode obtained in the second step,
Optimize the conformation of the ligand molecule so that the position of the hydrogen-bonding heteroatom in the dummy atom and the ligand molecule coincide,
Next, a third step of removing the conformation of the ligand molecule having a high intramolecular energy, (4) for each conformation not removed in the third step, the correspondence between the hydrogen-bonding hetero atom and the dummy atom in the ligand molecule A fourth step of obtaining the complex structure of the biopolymer-ligand molecule by replacing the coordinates of all atoms of the ligand molecule with the coordinate system of the biopolymer based on the relationship; (5) the complex obtained in the fourth step Removing from the structure a complex structure with a high intramolecular energy of the hydrogen-bonding moiety in the ligand molecule and a high intermolecular interaction energy of the biopolymer-hydrogen-bonding moiety in the ligand molecule; A fifth step of optimizing the structure of the structure; (6) for each complex structure obtained in the fifth step, a new complex structure is obtained by generating the conformation of the dehydrogenation binding portion of the ligand molecule 6th construction And (7) removing, from the complex structure obtained in the sixth step, a complex structure having a high intramolecular energy of the entire ligand molecule and a high intermolecular interaction energy of the biopolymer-ligand molecule, and then removing the remaining complex structure. A method for searching for a structure of a stable biopolymer-ligand molecule complex, comprising a seventh step of optimizing the structure of the complex structure.