Method of identifying candidate molecules

RELATED APPLICATIONS

The present application is a continuation-in-part of U.S. Ser. No. 09/422,303 filed Oct. 21, 1999, now abandoned.

The present invention relates to a method of identifying candidate molecules expected to be biologically active.

Evaluating receptor or target suitability of molecules is an important task in pharmaceutical drug research. With the increasing employment of automation techniques over the last years within Drug Discovery processes, methods like High-Throughput-Screening (HTS) and High-Throughput-Synthesis have become industry standards in pharmaceutical research. Nowadays, it is possible to test more than 20,000 molecules per day for their biological activities in certain disease targets. Also in the area of chemical synthesis, combinatorial chemistry in combination with automation processes, hundreds of molecules per day can be made physically available. Since based on today's chemical knowledge, more than 10

100

molecules could theoretically be synthesized and tested and several hundreds of thousands molecules are commercially available, computer assisted methods have been developed to select subsets of molecules which are actually supposed to be tested based on their predicted potential of biological activity for certain disease targets.

Two categories of computer assisted methods serve the purpose of discovering (selecting and/or prioritizing) molecules from data sets of theoretically available molecules for biological activity testing. The first category comprises diversity or similarity based discovery methods, whereas the second category comprises structure based discovery methods. Among the second category, there are database search techniques, as well as (Q)SAR methods and Docking methods.

Only the (Q)SAR methods and the Docking methods implicitly consider information related to specific targets, either common structural patterns of a series of active molecules ((Q)SAR) or the 3-dimensional structure of a target protein (Docking) and therefore deliver the most specific results. In practice, methods based on (Q)SAR or Docking are applied to smaller data sets (up to 50,000 sets), since they need relatively high computing power. However, although parallel computing techniques can be used to gain speed, still data sets consisting of more than 10

6

molecules are not predictable with respect to their biological activity in a reasonable time frame.

The term biological activity is hereinafter used to comprise in particular pharmaceutical as well as agrochemical activity with respect to a certain receptor or target.

The search for candidate molecules also comprises the search for lead compounds.

SUMMARY OF THE INVENTION

It is therefore an object of the present invention to provide a method of and a system for finding candidate molecules expected to be biologically active, which method and system can be applied on molecule libraries comprising high amounts of data and yields results in a reasonable time.

According to the invention, the method of identifying candidate molecules expected to be biologically active, comprises the following steps:

a) creating a set consisting of different molecules;

b) to each of said molecules of said set, assigning a descriptor representing a predetermined number of molecular properties;

c) mapping said set of molecules onto points of a two-dimensional grid with regard to a predetermined similarity relation of the respective assigned descriptors such that the grid distance between grid points of two molecules is a measure for the similarity of said two molecule descriptors;

d) forming a three-dimensional surface over said grid of molecules, said surface representing the distribution of biological activity of the molecules on the grid approximatively according to a predetermined quality criterion;

e) selecting from said three-dimensional surface candidate molecules satisfying a predetermined criterion with respect to their biological activity.

According to the invention, said three-dimensional surface in step d) may be formed by applying the following steps:

da) taking the whole two-dimensional grid as initial region of approximation;

db) selecting molecules on predetermined grid points of said region, and calculating their respective values of biological activity;

dc) approximating the surface over said region using said previously determined values of biological activity of said molecules on said predetermined grid points;

dd) determining whether said approximated surface satisfies a predetermined quality criterion; if so, goto step e); if not so, refining the approximation of the surface by selecting molecules on further grid points, calculating their respective values of biological activity, and repeating step dc) and this step dd).

The method used for forming said three-dimensional surface in step d) is preferably an approximation method of the Delauney Triangulation type.

Thus, the inventive method consists in performing two major steps. In the first step, the molecules are sorted and mapped onto a 2-D grid according to their similarity of their descriptors. In the second step, the biological activity of the mapped molecules is approximated by modelling the distribution as surface over the molecule map. From the surface, suitable candidate molecules for further evaluation can be determined. According to the invention, only a very small amount of molecules within the data set have to be really calculated. This results in a considerable gain of performance. The recursive proceeding allows to study the data base based on customizable quality criteria. Error and quality criteria for the analysis can be tailored according to a given problem. Docking simulations of collections of molecules can easily be run parallel, which additionally leads to a performance gain.

Thus, the method according to the present invention translates predicted/measured biological activity into topographical information of a three-dimensional surface, which is analyzed iteratively, using approximation algorithms. Only those regions of the surface are thoroughly analyzed, that represent regions of high biological activity while those of low binding energies for a given protein binding site are approximated by few data points. Thus, as examples have shown, active molecules can be identified from data sets by explicitly calculating/measuring just 4-6% of the molecules within the data set.

By using the method according to the invention, drug lead candidates can be identified without the need of making large molecule sets physically available and testing them.

The selected candidate molecules are suitable for chemical synthesis.

Preferably, the molecular properties represented by said descriptors are at least two of:

molecular weight,

number of rotatable bonds,

number of hydrophobic groups,

number of hydrophilic groups,

number of acid groups,

number of basic groups,

number of neutral groups,

number of zwitter groups,

number of heavy atoms,

number of H-bond donors,

number of H-bond acceptors,

number of 1-2 dipoles,

number of 1-3 dipoles,

number of 1-4 dipoles.

Preferably, said molecule-mapping is performed using self-organizing maps within a neural network or statistical methods like linear vector quantization.

BRIEF DESCRIPTION OF THE DRAWINGS

The invention comprises also a computer system having means for performing the identifying method, means for inputting commands to the system, and means for outputting the result of performing the method.

The invention and examples thereof are described in detail with reference to the accompanying figures, in which

FIG. 1

illustrates the type of descriptor used herein,

FIG. 2

illustrates the distribution resulting from a molecule sorting step,

FIG. 3

illustrates the de-clustering step,

FIG. 4

illustrates the distribution resulting from de-clustering,

FIG. 5

illustrates four intermediate steps of approximation of a surface by Delauney triangulation,

FIG. 5

a

illustrates the first intermediate step with 15 molecules of approximation of a surface by Delauney triangulation.

FIG. 5

b

illustrates the second intermediate step with 25 molecules of approximation of a surface by Delauney triangulation.

FIG. 5

c

illustrates the third intermediate step with 50 molecules of approximation of a surface by Delauney triangulation.

FIG. 5

d

illustrates the first intermediate step with 80 molecules of approximation of a surface by Delauney triangulation.

FIG. 6

illustrates the final approximation of the surface by Delauney triangulation,

FIG. 6

a

illustrates the final approximation of the surface of Delauney triangulation.

FIG. 6

b

illustrates the final approximation of the surface by Delauney triangulation.

FIG. 7

shows an a sorted data set on the map for an example of a data set of 40000 molecules,

FIG. 8

displays the amount of active molecules for a first target according the example of

FIG. 7

,

FIG. 9

displays the amount of active molecules for a second target according the example of

FIG. 7

,

FIGS. 10 and 11

show the percentage of active molecules related to the total amount of active molecules contained in the dataset for the results according to

FIGS. 8 and 9

, respectively.

DETAILED DESCRIPTION

The inventive method is illustrated using the example of finding inhibitors of Dihydrofolatreductase (PDB Code: 4 dfr).

According to the invention, prior to evaluation of particular molecules, the so-called virtual library is created, which comprises all possible molecules. That means that the virtual molecule library contains such molecules which can be purchased or produced with reasonable costs, that are commercially available molecules or molecules which can be produced using combinatorial synthesis approaches. Not comprised are molecules which are a priori not suitable for drug synthesis, in particular such molecules which contain toxic groups, or which have a molecular weight greater then 500 u or more than 5 donors, or molecules having a log P value of greater than 5. The library is held in a computer database.

The database in this example comprises 2,000 molecules. Each of the molecules is represented by 2-D structural data in a machine-readable form.

Upon storing the molecules in the library, a descriptor is assigned to each molecule within the library, which descriptor correlates with the biological activity of the respective molecule. The descriptor is a vector of several scalar molecular properties. This vector comprises the following values:

molecular weight,

number of rotatable bonds,

number of hydrophobic groups,

number of heavy atoms,

number of H-bond donors,

number of H-bond acceptors.

In order to perform a pre-selection of molecules, it is possible to use values covering economical or technical aspects, such as availability and production costs of molecules.

FIG. 1

displays four data streams (denoting four molecules) of the descriptor used herein. The first lines specifies the dimension of the descriptor, the second to fifth lines specify the molecules, the last column contains the ID of the corresponding molecules.

The descriptors are adapted for further processing the molecule library in order to find out the molecule candidates for drug synthesis. In order to allow further processing, the descriptors chosen for the molecules of the database are all of the same dimension.

In the next step, the molecules of that library are sorted according to their similarity with regard to potential ligand-receptor-activity. The sorting is performed by mapping the molecules onto a 2-dimensional grid, whereby the similarity of the molecules to each other are represented by the distances of their positions on the grid. The similarity is determined by evaluating the descriptors of the molecules. The sorting is based on a neural network of the Kohonen Self Organizing Map type. Self organizing maps are well-suited since they have been found out to be able to project higher dimensional problems to lower dimension even if there are linear dependencies in the data.

Thus, the sorting is performed in the following way: Initially, all the molecules are put in an arbitrary way onto the grid. The grid must be large enough to comprise all the molecules of the data base. During the learning phases of the neural network, the molecules are propagated over the grid with the aim to bring them into an order with respect to their similarity. Three learning phases are applied. In the first learning phase, the amount of propagation is limited to a grid radius of about 1/10 of the maximum number of grid points in one direction, whereby the convergence is 0.3. In the second learning phase, the radius is one grid point, and the convergence is 0.2; in the third learning phase the radius is 0.5, and the convergence is 0.2. It is to be understood that the grid radius and the convergence may be adapted to the concrete problem.

The result of the sorting is a molecule distribution over the grid, wherein a grid point can be populated by several molecules, whereas other grid points may remain free.

FIG. 2

illustrates the resulting map after three learning phases.

Although such an irregular population of the molecule map does actually not hinder the further processing of the molecule map, the molecule map is de-clustered, i.e., molecules which are placed together with others on the same grid point, are placed onto free grid points in their neighbourhood. This is accomplished by a search for a free grid point according to the following scheme (see FIG.

3

): The search begins at one of the four direct adjacent grid points (upper, lower, left, right grid point). The first one of these points is randomly selected. The molecule is placed thereon, whether that point is free or not. If that point has been free before the search is successfully terminated. Otherwise, if that point has been occupied before, the search for a free point continues with the next three adjacent grid points. If there is no free point found, the four diagonal neighbours (at the comers) are searched. If there is no free point neither, the respective molecule remains at the place where it has been placed, i.e. the first one of the adjacent grid points. For the point which has been initially at that place, the same search is performed. If after a predetermined number n of iterations there is one molecule left, that molecule is discarded. However, it has been empirically found that already for n=4, there is no molecule to be discarded.

Note that the border of the grid is taken as occupied points. If a border point is reached, the search is continued until the search is back within the grid area.

The result of this step is a distribution of the molecules over the grid wherein no grid point is occupied by more than one molecule.

FIG. 4

shows the resulting distribution. It is to be noted that this de-clustering method conserves the neighbourhood relations between the molecules.

Up to now, the inventive method is not specific for a certain receptor or certain target. The following steps, however, are specific for the receptor or target. Remember that the aim is to find out the candidate molecules with the highest biological activity with respect to the receptor or target. The biological activity of a molecule can be determined by many different approaches, namely calculations, simulations, or experiments.

Since the determination of the biological activity of all the molecules in the grid would be too much time-consuming, a strategy must be applied to reduce the number of possible candidates for further evaluation. Key idea is that adjacent molecules on the grid will have similar values of biological activity due to their similarity in properties (as expressed by the descriptors). However, note that concrete values of the biological activity are not directly known from the descriptors. These values must be determined by using the above-mentioned (very time consuming) methods.

The most straightforward approach to search those molecules having the highest values of biological activity over the molecule distribution would consist in computing all values of the molecule grid. However, such an exhaustive approach would be too much time consuming. Thus, an approximative modelling of the surface has to be performed.

An approximative modelling of the surface by determining activity values of molecules on regularly spaced grid points all over the grid would exhibit the danger of stepping over some interesting minimum/maximum values.

Therefore, an adaptive search is performed in which interesting regions of the grid are evaluated more precisely than other regions of the grid.

The distribution of the biological activity over the molecule grid is modelled as a 3-D surface representing the biological activity of the molecules over the grid. The problem of finding out suitable molecules is transformed into the problem of finding the grid points (i.e., the molecules) having the highest values of biological activity.

As most suitable approach for approximating the biological activity surface, the Delauney Triangulation method has been found. The Delauney Triangulation method consists in approximating a surface by iteratively decomposing it into triangles, whereby the number of iterations is increased at regions of the surface which are more shaped, i.e., regions with great differences in activity values.

Initially, the activity values of the four corner points of the grid are calculated. Based on these four points, the surface is approximated (i.e., decomposed) by two triangles. For each of the resulting triangles, a quality criterion is computed which defines the quality of approximation of the surface by the respective triangle. As quality criterion, the following expression has been found to be suitable:

DF

⁢

⁢

(

i

)

=

∑

j

⁢

⁢

ε

⁢

⁢

n

⁢

⁢

b

i

⁢

⁢

&LeftDoubleBracketingBar;

∇

i

⁢

-

∇

j

&RightDoubleBracketingBar;

·

∑

j

⁢

⁢

ε

⁢

⁢

n

⁢

⁢

b

i

⁢

⁢

&LeftDoubleBracketingBar;

f

⁡

(

j

)

&RightDoubleBracketingBar;

#

⁢

n

⁢

⁢

b

i

⁢

λ

⁢

⁢

(

i

)

·

κ

The quality criterion DF(i) is proportional to the sum of all gradient values ∇

i

, ∇

j

at the corner points i, j and to the average of the functional values f(j) over the number of corner points #nb

i

of the respective triangle i. Further, the quality criterion is inversely proportional to the surface area &lgr;(i) of that triangle. The factor &kgr; is a scaling factor, which is determined according to the needs for displaying the surface. Thus, using this quality criterion leads to a better decomposition of regions which comprise molecules having great differences in activity. Since the surface area of the triangles is taken into account, large triangles will be treated even if they do not comprise molecules having particular great differences in activity.

As further grid point whose activity is calculated, the center of gravity of the triangle having the worst quality factor (corresponding to the grid region under the surface having the greatest changes in activity) is selected. Based on this calculated point and all precedingly calculated points, a new decomposition of the whole grid surface into triangles is performed. Again, for each of the triangles, the respective quality factor is calculated in order to select a further point of calculation for refining the decomposition.

This step is iteratively performed.

FIG. 5

shows the respective decomposition for four iteration steps, namely

15

,

25

,

50

,

80

calculated molecule grid points. The iteration is stopped when the decomposition has become exact enough, i.e., when a predetermined stop criterion is reached, e.g., 5% of all grids points are calculated, or the General error function (sum of all quality factors) reaches a predetermined threshold. The grid surface can then be considered to be optimally approximated with respect to the interesting regions of the grid, i.e., the regions with grid points having the highest values of biological activity, as can be seen in FIG.

6

.

Then, in the last step, all the molecules within a predetermined radius around the grid points having the highest biological activity are explicitely evaluated.

All these points are the resulting molecules identified as candidates for further pharmaceutical research.

By applying this method, at least 94% of computing time with respect to an exhaustive evaluation can be saved. Further, due to the pre-sorting step, the optimum candidate molecules are nearly in every case found.

The identified molecules may be tested in suitable biological assays as described for instance by R. Bolger, “High-throughput screening: new frontiers for the 21

st

century”, published in DDT, Vol. 4, No 6, pp. 251-253, June 1999, or by J. S. Major, “Challenges of high throughput screening against cell surface receptors”, J. of Receptor and Signal Transduction Research, 15(1-4), pp. 595-607, 1995.

The method according to the invention was evaluated on a data set consisting of 40000 molecules from the World Drug Index (WDI). Some modifications compared to the basic method have been made in order to accelerate the processing.

The molecular descriptor was calculated for all the molecules. The sorting step took place on a rectangular 250×250 neuron SOM (i.e., on 62500 grid points). Instead of putting the molecules in an arbitrary way onto the grid (FIG.

1

), the weight vectors of the 62500 neurons of the self organizing map (SOM) were initialized according to the first eigenvectors of the complete molecule descriptor matrix, i.e., the 40000×6 matrix. Therefore, the de-clustering step as described above could be stepped over.

Then the training of the SOM was performed with a convergence of 0.2, and an amount of propagation of 15. A Gaussian neighborhood function was used. The learning was performed over 5 phases while the amount of propagation was decreased stepwise to 1 proportional to the elapsed learning time.

FIG. 7

shows the sorted data set on the map.

Then, the distribution of the biological activity was modeled. Prior to the triangulation steps, a convex hull around all molecules on the grid was calculated, by using 8 to 12 points of the grid. The resulting points were taken as a part of the initial data points. The other points used are the mean and the means of the four quadrants.

At these points, the biological activity of the underlying molecules was determined, and the initial triangulation was calculated. The docking calculations were performed for both thrombin and dihydrofolatreduktase (DHFR) receptors. For each of the resulting triangles, the quality criterion DF(i) was computed. The number of iterations for the triangulation was 2000.

As a last step, candidate molecules with a binding energy below −30 kJ/mol were selected and the biological activity of all molecules lying in the first and second neighborhood circle around these points was determined. The total amount of found molecules with an activity below −30 kJ/mol was calculated at this point and used as a measure for the performance of the algorithm.

The complete algorithm was applied to both receptors thrombin and DHFR. As explained earlier, the SOM needed to be calculated only once. For both receptors, the calculation was terminated after 2000 iterations, i.e. 5% of molecules were calculated.

The whole procedure including the subsequent refinement took about 4000 minutes for the Docking (ca. 2 minutes per molecule) and 12 hours for the SOM calculation, which had to be performed only once. This means an acceleration of about factor 8 for one investigation and an acceleration of 20 for all further runs for further targets based on the same SOM.

FIG. 8

displays the amount of active molecules (hits) found during the iterations for the thrombin receptor. The algorithm was finished after 2000 iterations. The horizontal gray line shows the total amount of hits present in the WDI for that specific receptor.

FIG. 9

display the amount of hits during the 2000 iterations for the DHFR receptor. The horizontal gray line shows the total amount of hits present in the WDI for that specific receptor.

FIGS. 10 and 11

show the percentage of hits related to the total amount of hits as a function of the number of iterations for thrombin, an DHFR (Dihydrofolatreductase), respectively.