KR20020028925A - New method for desiging molaculars of new medicines - Google Patents

Abstract

PURPOSE: A molecular design method is provided to construct chromosome with weighting values of probes as genes, to optimize an average weighing value of the probes by applying a least square method and a genetic algorithm alternatively or repeatedly, and analyze three dimensional quantitative structure activity relation so that it can estimate a physiological activity of unknown chemical compound. CONSTITUTION: The method comprises steps of generating probes for calculating probe interaction energy(100), generating initial objects by expressing weighting values of the probes as genes(200), obtaining a linear coefficient of chromosome by using a least square method for expressing a relation between the weighted probe interaction energy and the physiological activity(300), obtaining the average weighting value of the probes, and then obtaining the linear coefficient based the average weighting value(400), obtaining better weighting value of the probe by using a genetic algorithm(600), and obtaining the physiological activity by using spatial coordinates, a partial charge, and a final weighting coefficient of the probes(700).

Description

New method for desiring molaculars of new medicines

The present invention relates to a new molecular design method for drug development by analyzing the quantitative structure-activity relationship of the three-wave source, which is a molecular design method that improves the efficiency in the drug development process.

In general, developing a new drug requires an astronomical investment because it involves a complex process of synthesizing numerous substances to measure their physiological activity and verifying their toxicity.

Therefore, there is a global trend in approaching rational drug design, which understands pharmacological action at the molecular level and selects synthetic substances based on it.

Pharmacological action is a phenomenon caused by three-dimensional interaction between a target enzyme in the living body and a protein such as a receptor and a ligand, a molecule that has a medicinal effect. Structure-based drug design is used.

However, when not knowing the three-dimensional structure of the receptor, only the ligand molecules that already know the drug can be used to find out the quantitative relationship between the three-dimensional molecular structure and the drug, and from this relationship can predict the drug's efficacy by calculation. There is also a three-dimensional quantitative structure-activity relationship (3D-QSAR), and many universities and pharmaceutical research institutes around the world are actively developing and applying various three-dimensional quantitative structure-activity relationships. Is going on.

The most widely used software that uses three-dimensional quantitative structure-activity relations in the drug development process is the method invented by Cramer in 1988 and Tripos, Inc., St. Louis, USA. CoMFA (abbreviated for Comparative Molecular Field Analysis).

The algorithms used in this software are registered in US Pat. Nos. 5,025,388 and 5,307,287.

In this method, based on the general chemical theory that the logarithm of the ligand's potency or activity is proportional to the binding energy with the receptors, you overlap and lattice the space around it. Calculate the interaction energy between atomic probes and ligand molecules located at each point by selecting the points and choosing them as the values of the lattice points and contributing to the efficacy of each ligand molecule. coefficients) and then use these constants to predict the activity of the unknown compound.

At this time, the number of grid points, that is, the number of variables, is much greater than the number of molecules, that is, the number of measurements, so the partial least squares method, which is difficult to understand because of multiple regression, which is a common statistical method, cannot be obtained. A special algorithm called partial least squrres is used.

However, at least two types of probes must be used, and there are many mathematical problems such as variable selection in the partial least squares method and results vary depending on various environmental variables that the user must decide. The interpretation of the results is also indirect, rather than intuitive.

However, despite the problems of the above methodology and the disadvantage of being available only on expensive workstations, the comparator is widely used all over the world because of the effectiveness of the three-dimensional quantitative structure-activity method in the development of new drugs and the lack of alternatives. have.

In the present invention, it is possible to more accurately predict the activity of an unknown compound by using a different concept approach and a different mathematical solution than the conventional methodology including the compa, minimizing the user's decision making process and intuitive interpretation of the result. In addition, it is to provide a new molecular design method for new drug development by using a new three-dimensional quantitative structure-activity relationship analysis method that can be made with software that is short in computation time and can be used in general personal computers.

1 is a three-dimensional quantitative structure-activity relationship analysis system using the weighted probe interaction energy according to the present invention.

2 is a non-zero weighted probe and one compound state diagram obtained by training molecules of 21 steroid compounds that bind to corticosteroid binding protein.

3 is a graph showing experimental values of physiological activities of unknown molecules not included in a training molecule and calculations obtained using the probes of FIG. 2.

The three-dimensional quantitative structure-activity relationship analysis method using the weighted probe interaction energy of the present invention is based on the following theoretical background.

In the process of drug development, many compounds are synthesized to measure drug efficacy, and they exhibit different drug effects (IC ₅₀ ) depending on the molecular structure.

The reason is that the physiological activity A _i of the i-th compound can be expressed as the first-order function of E _i as shown in Equation (1) due to the difference in binding energy E between each ligand molecule and specific protein receptor such as enzyme. have.

Therefore, E _i is a calculable amount, so if the linear coefficients a and b can be determined, the biological activity A _i or IC ₅₀ can be calculated.

(1) physiological activity (A _i ) = -log (IC ₅₀ ) = a + b E _i

In the present invention, the drug varies depending on the interaction between ligand molecules and the surrounding environment. Therefore, when the expression of the environment is replaced by a set of probes, some probes are related to the difference in drug efficacy and some are not related. The hypothesis is that three-dimensional quantitative structure-activity relationships can be identified by determining the contributions of the probes to the interaction energy that may explain the difference in efficacy.

In other words, replace E _i in Equation (1) with the binding energy (or interaction energy) between the probes and the i-th ligand molecule (known as training molecules), and the weight of the j-th probe ( E _i can be represented by equation (2) when weight) is expressed as wj having a value from 0 to 1.

This concept is simple but is one of the key points of the present invention.

(2) E _i = ∑ _j ^probes w _j E _ij

Substituting Equation (2) into Equation (1) again yields the physiological activity A _i of the i-th compound (3).

(3) A _i ^cal = a + b∑ _j ^probes w _j E _ij

If we set a, b, and w _j to minimize the difference between the calculated value of A _i using the equation (3) and the value obtained from the experiment, it becomes the three-dimensional quantitative structure-activity relationship.

However, while there are only a few dozen experiments, the number of probes is at least hundreds, so setting a, b, and w _j is an underdetermined problem that is difficult to solve mathematically, as in a comparator.

In the present invention, this problem is solved repeatedly by using a linear least squares method and a genetic algorithm known as an optimization algorithm.

In other words, the weights of the probes w _j are expressed as genes with n binary bits per probe, and then a chromosome composed of these genes is generated. The common genetic algorithm is used to compare the calculated values of A _i with experimental values. We developed a new algorithm to find the w _j of each probe to minimize the difference.

Since the linear coefficients a and b can be determined by using the least-square method for each evolutionary step, the optimized w _j can be obtained by repeatedly applying the genetic algorithm and the least-squares method.

Since w _j is a real number value between 0 and 1, each weight w _j may be a chromosome composed of a gene represented by a real number, but if w _j is expressed as an n-bit integer instead of a real value, real Compared to constructing a chromosome, the program efficiency is increased and the calculation speed is much faster.

In other words, replace w _j with a two-bit integer that can represent 0,1,2,3, and divide it by 3 to express it as a number with intervals of 1/3 (that is, 0.3333) .If you substitute a 4-bit integer and divide by 15, It can be expressed as a real number of 1/15 (ie 0.0667) interval.

In general, constructing a gene with n-bit integers is equivalent to representing w _j by a real number of 1 / (2 ⁿ -1) intervals.

As a result of testing several conditions, it is best to express the weight w _j as a 4-bit integer in terms of the balance between the precision of the result and the speed of the calculation.

Once the optimized a, b, w _j are obtained, superposition of the three-dimensional structures of the undifferentiated molecules that have not been synthesized and tested in the superposition in the space with the training molecules is then performed using equations (1) and (2). Drug efficacy can be predicted in advance by calculation.

Hereinafter, the detailed calculation methodology of the present invention will be described in detail with reference to FIG. 1.

In the present invention, as in most three-dimensional quantitative structure-activity relationship analysis method, three-dimensional conformation of several training molecules known to be effective is used after superimposing them.

The probe generation and energy calculation unit 100 of FIG. 1 performs the following operations.

Two different solvent accessible surfaces are created, spaced apart from the overlapping molecular surface, and the atomic probes are set to have an interval of about 2 microns (angstroms) on the surface.

On the inner solvent contact surface, a probe to calculate the hydrogen bond energy is placed on the outer surface, and a probe to calculate the van der Waals or steric energy and electrostatic energy.

The sum of these three energies is the interaction energy between the ligand molecule and the probe.

The generation of probes is described in the Marching cube algorithm (Heiden, W .; Goetze, T .; Brickmann, J. Fast Generation of Molecular Surfaces from 3D Date Fields with an Enhanced "Marching Cube" Algorithm. J. Comput. Chem. 1993, 25, 246-250).

Hydrogen bond energy probes preselect only those capable of hydrogen bonding with ligand molecules and remove the rest from the probe list.

The partial charge of the probe required to calculate the electrostatic energy is equal in magnitude but equal in magnitude to the electrostatic potential of the ligand molecules, based on the general fact that there is a complementary electrostatic bond between the ligand and the receptor. Give it the opposite value.

Although any atom may be used as the probe, it was found that a simple result using only one type of atom can be obtained by various test results. According to the present invention, the oxygen probe is used at a distance of 1.6Å and 2.0Å away from the surface of the molecule. Cited.

In addition, experiments with various drugs have shown that the inclusion of hydrogen bond energy probes as well as van der Waals energy probes plays a significant role in improving the accuracy of the results.

The initial object generation unit 200 of FIG. 1 creates chromosomes randomly arranged 0 or 1 for each bit to represent the weight w _j of the probes by the population size, which is input data specified by the user. .

The initial individual evaluation unit 300 of FIG. 1 determines a good chromosome.

In this case, a fitness function for determining whether a chromosome including w _j causes a difference between the calculated value of A _i and the experimental value to be small is defined by Equation (4).

The smaller X _k ² with respect to chromosome k is the better chromosome.

(4) X _k ² = ∑ _i ^training [A _i ^obs -A _i ^cal ] ²

= ∑ _i ^training [A _i ^obs- (a _k + b _k ∑ _j ^probes w _jk E _ij )] ²

To obtain the X _k ² from the given w _j for each evolution Σ _j w _jk ^probes calculate E _ij and A _iobs and Σ _j w _jk ^probes E using a least square method to calculate a, b each chromosome from _ij.

In the mean weight determination unit 400 and the convergence check process 500 of FIG. 1, the average value of w _j is obtained by using the Boltzmann equation from only chromosomes, elites, of which X _k ² is smaller than one chromosome. To obtain the final w _j and obtain the a, b of the whole system from them. Also, if X _k ² is no longer smaller, stop the evolution and change w _j to reduce the evolution if X _k ² can be further reduced. Perform convergence checks to continue.

The genetically engineered entity generator 600 of FIG. 1 performs an evolution process of manipulating a gene to change w _j of each probe in a direction in which X _k ² decreases.

Evolution may be any genetic algorithm, but the present invention uses the following simple method.

That tournament selection method to select the two chromosomes at random using the (tournament selection) performed twice through the process of selecting that X _k ² smaller as parents of the next generation, and even crossed the two parents became chose (uniform crossover) to the next generation The process of making two children of (generation) was repeatedly performed to make the desired number of children, and then some of them were mutated.

By repeating the processes of the individual evaluation unit 300 and the genetically engineered individual generation unit 600 of FIG. 1 for hundreds of generations until convergence is achieved, the optimized values can be given closest to the experimentally measured physiological activities. The average weights w _j of the probes and the linear coefficients a, b can be obtained.

In the physiological activity predictor 700 of FIG. 1, the training molecules are used to train the three-dimensional structure of the unknown molecules using the spatial coordinates, partial charges, weights (w _j ), and linear coefficients a and b of probes whose weight is not zero. After superpostion and space, the physiological activity of unknown molecules is calculated using Eq. (3).

Experiments with several drugs resulted in weights of probes of about 70-80% of the original probe number being zero.

The remaining set of non-zero probes may be called a hypothetical receptor.

If these probes are shown as molecules using three-dimensional computer graphics, one can intuitively see how changing any part of the molecule reduces the interaction energy and improves its efficacy, making it an important tool for molecular design.

All of the calculations in the methodology were made in computer programs written in C and can be used on personal computers.

Using a program called the WeP (weighted probe intraction energy method), a file that specifies a three-dimensional atomic coordinate file and a minimum user variable (initial population size) can be found. You can get it automatically.

Actually, the test results using the WeP program are shown in FIGS. 2 and 3.

Test samples are a variety of steroid compounds that bind to corticosteroid binding globulin (CBG).

These CBG steroids are described by Cramer in a Comparator study (Cramer III, RD; Patterson, DE; Bunce, JD Comparative Molecular Field Analysis (CoMFA). 1.Effect of Shape on Binding of Steroids to Carrier Proteins.J. Am. Chem Soc. 1988, 110, 5959-5967) has been the benchmark sample of almost all 3D-QSAR methodology studies and has been used to compare the results.

The physiological activity, that is, 21 steroid compounds of known efficacy, is used as a training molecule to obtain a structure-activity relationship, and the same is used to calculate the effects of 10 steroid compounds, which are also known as drugs but not used as training molecules.

The smaller the difference between the calculated value and the experimental value, the better the method.

The standard deviation of the error of predictions (SDEP) is often used as a parameter to express this difference quantitatively. The smaller the SDEP, the better the result.

As a result of WeP calculation using the structure and pharmacological data of 21 steroid compounds in Cramer's research paper, the initial total number of probes is generated 240, but as a result of optimizing the weight of probes by genetic algorithm, non-weight probes Was 40.

The 40 probes have different weights and partial charges, which are shown in FIG. 2 using a computer graphics program with one of the training molecules in a sphere proportional to the weights.

The results of calculating the physiological activities of 10 steroidal compounds using known 40 probes but not included in the training molecule are shown in FIG. 3.

If the calculated physiological activity (physiological activity) is exactly the same as the experimental value, all the points in this plot are on a straight line of y = x plotted on the plot.

The SDEP of the WeP result is 0.547.

Recently published papers (Robinson, DD; Winn, PJ; Lyne, PD; Richards, WG Self-Organizing Molecular Field Analysis: A Tool for Structure-Activity Studies.J. Med. Chem. 1999, 42, 573-583) We compared the results obtained by using a primitive comparator with six new three-dimensional quantitative structure-activity relationship methods that improved or modified it. The SDEP of the compa was 0.837 and the SDEP of the other methods was 0.640 to 0.716. Only the SDEP of the method called SOMFA (Robin's paper, above) was 0.583.

Therefore, it can be seen that the prediction of the drug efficacy using the WeP method giving a SDEP of 0.547 has the best accuracy among the existing methods.

According to the present invention, the WeP method is more accurate than the conventional three-dimensional structure-activity quantitative relationship identification method, and does not use a complex statistical analysis method such as partial least squares method that requires cross validation. It is easy to understand, easy to interpret the results, requires minimal user decision making process and input data, and the results can also be used intuitively in molecular design using computer graphics.

Therefore, this WeP program can help increase competitiveness by saving time and cost for synthesis and drug efficacy test in new drug development process.

In addition, since it uses a concept approach and a mathematical solution that are different from the three-dimensional structure-activity quantitative relationship known methods, the development of commercial software to cope with the monopolization of the world software market in this field is also possible. It will be possible.

Claims (8)

Probe generation and energy calculation unit 100 for calculating the probe interaction energy; Initial entity generation unit 200 expressing the weight of the probe in the gene of the chromosome, entity evaluation unit 300 for determining the linear coefficient of the chromosome representing the relationship between the weighted probe interaction energy and the physiological activity by the least square method, and the average of the probes An average weight determiner 400 for determining the weight and linear coefficients from the least square method, a genetic manipulation entity generator 600 for obtaining better probe weights using a genetic algorithm, spatial coordinates, partial charges, New molecular design method for drug development, characterized in that consisting of the physiological activity predictor 700 to obtain a calculation of physiological activity using the final weight and the linear coefficient. The method of claim 1, wherein the probe generation and energy calculation unit 100 has two different solvent contact surfaces for placing the probe to calculate the hydrogen bonding energy on the inner surface and the probe to calculate the van der Waals energy and electrostatic energy on the outer surface A new molecular design method for drug discovery, which consists of making a drug. The method of claim 1, wherein the initial individual generation unit 200 consists of a real number obtained by dividing a gene represented by an integer using a binary n-bit by (2 ⁿ -1) by the weight of the probe and obtaining a chromosome from these genes. New molecular design method for drug development characterized in that. The method according to claim 1, wherein the individual evaluation unit 300 is configured to determine the linear coefficient for each chromosome by the least square method under the assumption that the weighted probe interaction energy has a linear relationship with the physiological activity. New Molecular Design Method. The method of claim 1, wherein the mean weight determination unit (400) is configured to obtain the average weight of the probe and to determine the linear coefficient of the entire system from the least square method from the new molecular design method for drug development. The method of claim 1, wherein the genetically engineered object generation unit 600 genetically manipulates chromosomes consisting of genes expressed in n-bit to obtain a new molecule design for drug development, characterized in that to obtain chromosomes consisting of a better probe weight Way. The physiological activity predictor 700 uses the spatial coordinates, partial charges, average weights, and linear coefficients a, b of probes whose weight is not zero to determine physiological activities of unknown molecules in which three-dimensional atomic coordinates overlap. A new molecular design method for drug discovery, characterized by computing. The method of claim 1, wherein the individual evaluation unit 300, the average weight determination unit 400, and the genetic manipulation individual generation unit 600 alternately and repeatedly use a least square method and a genetic algorithm to obtain an optimized probe weight. New molecular design method for drug discovery, characterized in that consisting of.