Quantitative Pharmacophore Models with Inductive Logic Programming

Quantitative Pharmacophore Models with Inductive Logic Programming Ashwin Srinivasan1, David Page2, Rui Camacho3, and Ross D. King4 1IBM India Research Lab, Block 1, Indian Institute of Technology, New Delhi. 2Department of Biostatistics, University of Madison, Wisconsin. 3LIACC - CIUP, R. Campo Alegre, 4150 Porto, Portugal. 4Department of Computer Science, University of Wales, Aberystwyth. Abstract. Three-dimensional models, or pharmacophores, describing Euclidean constraints on the location on small molecules of functional groups (like hydrophobic groups, hydrogen acceptors and donors, etc.), are often used in drug design to describe the medicinal activity of potential drugs (or ‘ligands’). This medicinal activity is produced by interaction of the functional groups on the ligand with a binding site on a target protein. In identifying structure-activity relations of this kind there are three principal issues: (1) It is often difficult to “align” the ligands in order to identify common structural properties that may be responsible for activity; (2) Ligands in solution can adopt different shapes (or ‘conformations’) arising from torsional rotations about bonds. The 3-D molecular substructure is typically sought on one or more low-energy conformers; and (3) Pharmacophore models must, ideally, predict medicinal activity on some quantitative scale. It has been shown that the logical representation adopted by Inductive Logic Programming (ILP) naturally resolves many of the difficulties associated with the alignment and multiconformation issues. However, the predictions of models constructed by ILP have hitherto only been nominal, predicting medicinal activity to be present or absent. In this paper, we investigate the construction of two kinds of quantitative pharmacophoric models with ILP: (a) Models that predict the probability that a ligand is “active”; and (b) Models that predict the actual medicinal activity of a ligand. Quantitative predictions are obtained by the utilising the following statistical procedures as background knowledge: logistic regression and naive Bayes, for probability prediction; linear and kernel regression, for activity prediction. The multi-conformation issue and, more generally, the relational representation used by ILP results in some special difficulties in the use of any statistical procedure. We present the principal issues and some solutions. Specifically, using data on the inhibition of the protease Thermolysin, we demonstrate that it is possible for an ILP program to construct good quantitative structure-activity models. We also comment on the relationship of this work to other recent developments in statistical relational learning. 1 Introduction The primary goal of the pharmaceutical industry is to find, develop and market new drugs for previously untreatable diseases or which have better properties than existing drugs. The development of a new drug is a time-consuming, expensive process—it can take anywhere from 12 to 16 years from the start of a research program to final approval; and the cost can be up to US$700 million. It is clearly of significant humanitarian and commercial interest to investigate techniques and tools that can assist in making the process more efficient. Most drugs molecules work by binding to “target sites”—commonly proteins— within the body. By interaction with these targets, drugs can modulate their actions. The process of drug development begins with the selecting an appropriate target with which the drug could interact to modulate disease. In some cases, the detailed structure of the target binding site is known (by the use of X-ray crystallographic techniques) or can be guessed (by structure prediction techniques or from the structure of similar molecules); but this situation is still rare and the identification of potential drugs has to proceed without this knowledge. Chemistry for drug development (see Fig. 1) is concerned with the identification of “ligands”: small molecules that are potential drugs. The search for ligands commences with some compounds that are known to interact with the target. These compounds, or “leads”, may be identified in a number of ways, for example from large scale empirical testing of available chemicals. They will have some activity, i.e. ability to interact with the target, but may do so only weakly, and may possess other undesirable properties, for example metabolic instability. Once a set of leads have been isolated, their chemical structure can be refined to improve biological activity and reduce side effects. A computational activity that assists this is concerned with constructing models that relate activity to molecular structure. These “structure-activity relationships”, or SARs, describe how the structural differences of a set of leads affects their activity, and can be used to suggest new molecules to make which should have enhanced activity. Chemical synthesis and testing of such molecules leads to the discovery of more active molecules until, ultimately, ligands of the desired activity and properties are discovered. These compounds now begin the long process of development including biological trials for efficacy, safety trials, formulation and advanced testing, patent development and finally, application for clinical trials. The whole process is critically dependent on effective initial chemistry: proposal of a large number of ligands for development and testing, most of which turn out to be useless is clearly to be avoided. In this, good SAR models can play an important role. “Traditional” SAR methods (“1-D” and “2-D” methods in Fig. 1), generate many features for each lead (effectively, a form of domainspecific propositionalisation) and use statistical modelling tools on the resulting, table containing the feature-values for a set of leads and their activity. The features, or molecular descriptors, can be properties of the whole molecule, such as partition coefficient, molecular volume, number of rings, etc., and, if the series of leads share a common core, properties of the substructures at the positions of variation The modelling task is to construct a predictive model that relates Chemistry Pharmacology, Safety etc.Research Compounds Synthesise Get Leads Potential Activity Predict Protein to inhibit Choose Target Biological Trials ‘‘Grid’’−based Structure−Activity Relationships (SARs) 1−D 2−D 3−D Target−Ligand Models Ligand−only Models Models Pharmacophore Models Qualitative Quantitive Fig. 1. A simplified view of the chemistry involved in the drug-design process (adapted from a figure kindly provided by Dr. Stuart Green, University of Leeds). The focus of this paper is shown in bold-face. feature-values to activity. Such models have been widely applied, and there are many examples of successful SAR analyses [12]. There are four principal difficulties with this approach. First, much depends on the ability to identify molecular descriptors that retain all the information necessary for obtaining a good model. Second, constructing reliable models is not easy. The molecular descriptors can number into several millions, and special care must be taken avoid chance correlations. Third, the models are not always easy to use in refining leads. It is easy to calculate the value of a molecular descriptor for a given lead (for example, an electronic partial charge), but quite a different matter to design a molecule that will possess that value. Finally, by dealing with bulk molecular properties, rather than a more explicit representation of the molecules, no account is taken for the three-dimensional shape of the leads. Leads, being small molecules are flexible and can adopt different shapes (conformations) by torsional rotations about bonds within the molecule. Molecules rapidly convert from one conformation to another, so that no single conformation can be isolated and tested for a given activity. Hence in general one does not know a priori which conformation of a molecule is an active conformation for a given form of activity. Computational chemists usually employ one of two approaches to this problem. In the first, the shape and electrostatic interactions are dealt with via a calculated interaction with a “probe” atom or group at points on a three-dimensional grid which surrounds the molecules. Statistical methods are then used to identify those parts of the molecule which are responsible for activity. The analysis can be displayed graphically to aid in the design of new molecules. The major disadvantage with this approach is that, in order to compare values at the calculated points between molecules, the conformation of the molecules, and a common coordinate frame, or alignment, must be chosen in advance of the analysis. This is equivalent to deciding on the manner in which the molecules interact with the target. If the molecules contain a large common structural element this alignment may be straightforward, but this is often not the case. The second approach uses a representation of biological activity called a pharmacophore. This is an abstraction of the molecular structure to the, usually, small number of key features which contribute the majority of the activity, together with their geometric arrangement represented by pairwise distances (see Fig. 2, from [8]). An advantage of the pharmacophore representation is that it expresses biological activity in a language that is familiar to chemists within the pharmaceutical industry. These representations are also readily convertible into search queries of compound databases, which is an effective means of identifying additional active molecules [9]. That Inductive Logic Programming (ILP) is particularly well-suited to represent and discover SARs has been argued persuasively elsewhere [8, 20]. The most important reasons given are: (i) the comprehensibility of the models constructed allow an easy translation of the models into new chemical structures; (ii) ILP programs do not need to align the leads to a common spatial reference frame; and (iii) the first-order representation used by ILP naturally deals with d1 d2 d3 Target Lead Fig. 2. A schematic diagram of how a “pharmacophore” definition (right) is related to a lead-target interaction (from [8]). The pharmacophore contains the key functional interactions and the geometric relationships between them expressed as distances d1, d2, and d3. For example, the functional groups may be a hydrogen donor and two hydrogen acceptors; d1 might be 4.5 Angstroms, d2 might be 5.0 Angstroms, and d3 might be 3.75 Angstroms. difficulties arising from the fact that the leads may assume several shapes (or conformations), not all of which may be responsible for the biological activity. Non-relational learning algorithms—those that need a feature vector represention of data—require modification to handle point iii. Furthermore, given even just two shapes for each molecule, these approaches either require worst-case exponential time (in the number of molecules) to perform the alignment mentioned in point ii, or must at times work with an incorrect alignment—one in which the parts of the molecules responsible for interaction are not aligned with each other. Regarding point i, comprehensibility, the most complex SARs found by ILP using only the structures of leads and their activities are pharmacophore models in the form of rules like the following: Molecule M is “active” if: M has a hydrogen donor at position P1 and M has a hydrogen acceptor at position P2 and M has a hydrogen acceptor at position P3 and the distance between P1 and P2 is 4 ±1˚ A and the distance between P1 and P3 is 3 ±1˚ A and the distance between P2 and P3 is 5 ±1˚ A. Here, positions P1, P2, and P3 are points in 3-dimensional space and the rule represents a “3-point pharmacophore”. The rule, evidently simple to understand, nevertheless highlights an important shortcoming of SARs constructed by ILP: they are classificatory in nature and unable to quantify their predictions of ac- tivity1. In this paper, we investigate identification by ILP of two kinds of quantitative SARs, namely, class-probability models of the form: The probability of molecule M being “active” is P if: M has a hydrogen donor at position P1 and M has a hydrogen acceptor at position P2 and M has a hydrogen acceptor at position P3 and the distance between P1 and P2 is X1 and the distance between P1 and P3 is X2 and the distance between P2 and P3 is X3 and P is the probability of being “active” given X1, X2, X3. and regression models of the form: The affinity of molecule M is A if: M has a hydrogen donor at position P1 and M has a hydrogen acceptor at position P2 and M has a hydrogen acceptor at position P3 and the distance between P1 and P2 is X1 and the distance between P1 and P3 is X2 and the distance between P2 and P3 is X3 and A is the expected affinity given X1, X2, X3. We will call these “quantitative pharmacophore models”: specifically, the rules above are 3-point quantitative pharmacophore models. The reason for interest in the regression model is evident: when accurate quantitative information of lead-activity is available, SARs that relate lead-structure to activity are most useful. Even when we only have categoric information about the lead-activity (for example, leads are “active” or “inactive”), class-probability models can still be extremely useful for the reasons below: – Inherent uncertainties in the domain may make categoric classification difficult. For example, data may be laboratory measurements obtained from an imprecise assay; – Applications involving decision-making often require probability estimates for use in cost/benefit calculations. For example, a synthesis of a particular molecule may proceed only if there was a very high probability of it being biologically active; – It may be important to rank alternatives with a given class value. For example, laboratory constraints may require synthesis of molecules to proceed in small batches. An ordering on molecules predicted to be active is then needed. 1These classificatory rules have been used to construct boolean features, which have then be used by regression techniques to build quantitative models for activity [21, 37]. We are concerned here with constructing quantitative models solely with the use of ILP. In this paper, both the geometric constraints and the “numeric” constraint in the rules above (that is, computing P or A given X1, X2 and X3) are to be obtained during hypothesis construction by an ILP program. It is our intention to investigate the use of statistical procedures as background knowledge for obtaining the numeric constraint. The following questions—novel for both ILP and statistical methods—arise for rules above: Model Identification. Molecule M may have several hydrogen donors and acceptors, at different positions, resulting in different sets of values for X1, X2 and X32. The novelty for ILP is that the domain dictates only one of these sets is relevant for identifying the the appropriate numeric constraint. Which one? Prediction. Each set of values for X1, X2 and X3 results in a prediction for P or A. Which of these values should be returned3? This question has not arisen previously in either ILP or statistical methods. These two tasks are unlike those confronted by a statistician during normal discourse. There it is common, due to random variation, to obtain different affinity values for repeated measurements of the same values for X1, X2 and X3. Machine learning researchers will recognise the questions here as consequences of a multiple-instance representation [6]. Novel here, however, is the real-valued prediction task in a multiple-instance setting. While this combination has been addressed before, in two papers [32, 1], as discussed in the next paragraph an adequate solution to model identification and prediction in this setting has not been found. And the combination of real-valued prediction in a multiple-instance setting has not been addressed explicitly before in ILP, although the paper by Ray and Page [32] pointed to its importance. The approaches proposed in the two papers already mentioned performed well on model identification from synthetic data. Note that model identification is easy to assess on synthetic data, where the data generator can record the model it used to generate the real-valued response for each data point. Model identification is difficult or impossible to assess on real-world data, where one sees only the response value and not the correct model, or set of bindings for variables. Hence one must rely on prediction to assess performance on most realworld tasks. The approaches in the two papers have not been demonstrated to perform well on real-world tasks, such as predicting quantitative drug activity values. The present paper presents procedures to address both model identification and prediction, and it tests these procedures on real-world data for quantitative drug activity prediction. For illustrative purposes, we will use linear regression as 2The logic programmer will recognise this as different substitution-sets arising from a non-determinate definition for the donor and acceptor predicates. 3We could consider using some representative set of values for distances (for example, the average) for both model identification and prediction. However, what this representative should be may not be apparent and in any case, will require a different rule to the ones shown. a background predicate to perform quantitative prediction, although the procedures proposed are not restricted to that particular modelling technique (as will be demonstrated in Section 6.1). The paper is organised as follows. In Section 2 we present a simple example to illustrate the principal issues. Section 3 introduces relevant terminology from the literature on multiple-instance learning. “Solutions” to the problems that arise during model identification and prediction are in Sections 4 and 5 respectively. Application of the procedures developed to quantitative structure-activity relations is in Section 6. Section 7 concludes the paper. The paper is accompanied by two appendices. Appendix A describes the statistical procedures used in the empirical study. Appendix B relates the work in Section 6.1 to other, more general, work on incorporating probabilities in ILP and relational learning generally. 2 An Example Consider learning the following simple 1-point quantitative pharmacophore model using an ILP system: The affinity of molecule M is Y if: M has a hydrogen donor at position X and Y = mX + c. Here “Y = mX + c” is a statistical model in which mand care parameters to be estimated and Suppose data available are as follows: (1) affinity values of 5 different molecules; and (2) records of hydrogen donor locations on each molecule (a molecule can have more than one donor). The data are tabulated in Fig. 3. Logically speaking, the rule above states that the affinity of the molecule is a linear function of one of the hydrogen donor positions (without specifying which one: the logician will recognise this as a consequence of the variable X being existentially quantified within the rule body). Since for each N, it is not apparent which of the X values are to be paired with the corresponding Y, a combinatorial problem arises4. For this simple problem, the combinatorics are manageable: there are only 3×1×2×1×2 = 12 different tables containing exactly 1 entry for each of the 5 values of X. Two such tables are shown in Fig 4. Parameter estimates can now be obtained using each such table. The estimates returned are those that result in the best fitting model. For the data in Fig 4, the best fitting model Y = 2.572 X - 33.155 is obtained with the table in Fig 5. The resulting rule is therefore: The affinity of molecule M is Y if: M has a hydrogen donor at position X and Y = 2.572 X - 33.155. 4Some theoretical results known for multiple-instance learning can be found in [7]. Mol (M) Aff. (Y) Donor (X) 25 1 50 33 29 2 100 50 3 150 73 75 4 200 90 5 250 120 110 Fig. 3. “Training” data for the problem of predicting affinity Y using hydrogen donor location X. For a given value of N, it is known that the value of Y depends on one of X values—but it is not known which one. Mol (N) Aff. (Y) Donor (X) Mol (N) Aff. (Y) Donor (X) 1 50 25 1 50 33 2 100 50 2 100 50 3 150 73 3 150 75 4 200 90 4 200 90 5 250 120 5 250 110 Fig. 4. Example tables obtained by combinatorial enumeration of the data in Fig. 3. The table on the left yields the model Y = 2.161 X - 4.739, with sample correlation coefficient 0.997. The one on the right yields Y = 2.565 X - 33.677, with sample correlation coefficient 0.998. Section 4 describes procedures that equips an ILP system to reach such a model (or at least, an approximation to it) given the data in Fig 3. Now consider using this rule on the data shown in Fig. 6. It is evident that simply executing the rule above will yield two predictions: Y = 188.037 (using X = 86) and Y = 203.469 (using X = 92). Which of these predictions should be returned is the problem posed by (Q2). In Section 5, we examine some solutions to this problem. 3 Terminology It is convenient at this point to introduce some terminology from the machine learning literature on multiple instance learning. The following statement of the multi-instance learning task is largely from [32]. Data consists of a set of nbags (1 ≤n < ∞). The ith bag consists of miinstances (1 ≤mi<∞) and a label yi(which may be nominal or real-valued). Instance jof bag iis described by a d-dimensional vector of values (1 ≤d < ∞). The task is to construct a model 6 Quantitative Pharmacophore Models for Thermolysin Inhibition Our test-bed for the construction of quantitative pharmacophore is the inhibition of Thermolysin. Thermolysin is a zinc-containing protease that consists of two spherical domains separated by a deep cleft that constitutes the active site. Zinc-containing proteases like Thermolysin play an important role in physiological processes like digestion and blood pressure regulation. Data on a number of inhibitors of Thermolysin are readily available in the literature: we use the molecules studied by [21]. The data consist of crystal structures and corresponding activity values (pKi = -logKi) of 31 inhibitors. All experiments use the ILP system Aleph [35] (specifically, Aleph Version 4). The experiments were performed on machine equipped with two 512 Mhz Pentium III processors with 128 megabytes of random access memory. We follow the work of [8] in providing the ILP system with background knowledge of the following: –Compound-specific knowledge. This is in the form of the atom and bond structure of each compound, as well as its 3-dimensional conformation (for each of the three lowest energy conformers identified by [21]). This information is represented by first-order atomic formulae: we refer the reader to [8] for examples of this encoding. –General chemical and geometric knowledge. Generic chemical knowledge provided is in the form of a library of elementary chemical concepts (for example definitions of hydrogen donors, hydrogen acceptors, hydrophobic groups, esters, ethers, etc.). We have provided the same concepts as those in [21] (the authors there have used 39 such concepts). The only geometric knowledge needed for this task is a procedure for calculating the Euclidean distance between two points. Additional background knowledge required for the construction of class-probability and regression models will be described below. 6.1 Class-probability models It is our goal in this section to demonstrate that an ILP system capable of using the procedure in Fig. 7 and statistical procedures for conditional probability estimation can construct class-probability models for structure-activity prediction. Recall that these models were of the form: The probability of molecule M being active is P if: M has a hydrogen donor at position P1 and M has a hydrogen acceptor at position P2 and M has a hydrogen acceptor at position P3 and the distance between P1 and P2 is X1 and the distance between P1 and P3 is X2 and the distance between P2 and P3 is X3 and P is the probability of being active given X1, X2, X3 Candidate statistical procedures for computing P considered are logistic regression and “naive” Bayes (see Section A.1). The reasons for selecting these techniques are: –They are simple; –There is abundant support for their use in class probability estimation, both in the statistical and machine learning literature (see for example, [29]). In addition, it can be shown that the logistic function is the appropriate Bayesian choice under some fairly general conditions (see [17]); and –They can be seen as special cases of more general work on combining Bayesian networks and ILP (see Appendix B). Data Of the 31 Thermolysin inhibitors available, we have designated the top 15 inhibitors to be be “active” and the remaining 16 inhibitors to be “inactive”. In order to build class-probability models, the actual examples are of the form “probability of molecule mbeing active is p”, where pis one of 1.0 or 0.0 (depending on whether mis active or inactive). Additional background knowledge In addition to the background information described earlier we also include the following: –Constraints on legitimate models. Legitimate models are required to contain functional groups, pairwise distances between the groups and a numeric constraint that predicts the probability of being active. –Model identification. This definition implements the procedure described in Fig. 7. The actual model construction within this procedure is done by one of logistic regression or naive Bayes. Method Our method is straightforward. Using background definitions of logistic regression or naive Bayes: 1. Construct the “best” class-probability model for the largest pharmacophore possible for the inhibitors. 2. Estimate the performance of the model obtained in the previous step. The following details are relevant: (a) Looking for the largest pharmacophore model for a set of active molecules is a characteristic of this kind of domain (see [8] for more details). Models can be seen as containing two sorts of constraints: first, those in the “prefix”, consisting of the functional groups and their pairwise distances (the pharmacophore per se); and second, the numeric constraint that predicts probabilities using the distances between these functional groups. In Step (1) we restrict models to be a single clauses. This is similar to [8] and offers the most comprehensible models for activity). Further, the search for this model is done in 2 steps:(i) the largest prefixes common to all inhibitors are found; and (ii) each prefix is extended to include the numeric constraint. For the Thermolysin data, the largest pharmacophore models contain 4-point pharmacophores (that is, there is no single clause 5-point pharmacophore model that can be used to explain the activity of all of the 31 inhibitors). The search procedure in Fig. 7 searches for the “best” numeric constraint by minimising a loss-function. Here we will use a standard quadratic loss function (sum of squared differences between actual and predicted probabilities) summed over all training instances. (b) In Step (2) performance estimates will be estimated using a leave-one-out procedure. That is, each inhibitor is, in turn, categorised as a “test” instance. A class-probability model is constructed using the remaining instances and then used to predict the probability of the test instance being active. The non-determinate nature of the domain will result in a set of predictions for the probability. We have elected to select the highest of these as the final prediction (as described in Section 5.2). This is a reasonable choice for the domain considered, as we are interested in the best chance of a molecule being able to inhibit Thermolysin (if the set of predictions is empty, then the test instance is not included in estimating the performance). The performance is summarised by an ROC curve, obtained by changing the threshold probability above which the test instance is to be classified as being active (for example, if this threshold is 0.3, and the model predicts the probability of the test instance being active as being 0.2, then the test instance is classified as being inactive). Results Figure 10 shows ROC curves generated by changing threshold probabilities in the manner described above. The curves are, in fact, the convex hull of the points denoting classificatory models obtained from the corresponding class-probability models (it has been shown elsewhere [31] that the convex hull contains the optimal classifiers under some very general decision-theoretic conditions). The curves clearly show how we can obtain higher positive prediction rates at the expense of an increased false-positive rate: obtaining a set of models with such properties is straightforward once we have class-probability models. For false-positive rates above 0.2, the ROC curve obtained with naive Bayes dominates that with the logistic regression. That is, as long as a false alarm rate of at least 20% is tolerable, the use of naive Bayes as a background predicate results in better predictive performance. Otherwise, the use of logistic regression yields a better model. The position is more clear-cut on the comprehensibility front: the use of logistic regression results in a simple equation for predicting probability. In contrast, the when using the naive Bayes procedure, probabilities are kernel-density estimates obtained from the training instances. This makes the corresponding rule returned by the ILP program significantly harder to un- 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 True Positive Rate False Positive Rate NBayes Logistic Random Fig. 10. ROC curves for Thermolysin Inhibition. “NBayes” is the curve obtained using a naive Bayesian procedure for probability estimation. “Logistic” is the corresponding curve using logistic regression as background knowledge. “Random” is the performance of a classifier that assigns class-values using a coin toss. derstand. Example rules obtained with logistic regression and naive Bayes are shown in Figures 11 and 12. 6.2 Regression models The structure-activity relations constructed so far do not include predictions of actual binding affinities. This is usually harder and in this section we demonstrate that an ILP system capable of using the procedure in Fig. 7 and appropriate statistical procedures for conditional mean estimation can construct a regression model that predicts binding affinities. Recall that an example is a rule of the form: The affinity of molecule M is A if: M has a hydrogen donor at position P1 and M has a hydrogen acceptor at position P2 and The probability of molecule M being “active” is P if: M has a hydrogen donor at position P1 and M has a hydrogen donor at position P2 and M has a hydrogen donor at position P3 and M has a hydrogen donor at position P4 and the distance between P1 and P2 is X1 and the distance between P1 and P3 is X2 and the distance between P1 and P4 is X3 and the distance between P2 and P3 is X4 and the distance between P2 and P4 is X5 and the distance between P3 and P4 is X6 and P = 1 1+e−87.5+16.1X1−6.4X2+0.1X3−36.3X4+43.8X5−80.1X6 Fig. 11. An example of a class-probability model obtained with the use of logistic regression. The regression procedure computes the numbers in the equation for P. M has a hydrogen acceptor at position P3 and the distance between P1 and P2 is X1 and the distance between P1 and P3 is X2 and the distance between P2 and P3 is X3 and A is the expected affinity given X1, X2, X3 Once again our choice of statistical procedures—linear and kernel regression (see Section A.2)—is primarly based on simplicity and prevalance of use. Data Data consists of the 31 Thermolysin inhibitors along with their activity (pKi) values. Examples are thus of the form “the affinity of molecule mis y” where yis some floating-point number. Additional background knowledge In addition to the background information described earlier we also include the following: –Constraints on legitimate models. Legitimate models are required to contain functional groups, pairwise distances between the groups and a numeric constraint that predicts the affinity. –Model identification. This definition implements the procedure described in Fig. 7. The actual model construction within this procedure is done by one of linear or kernel regression. Method Our method is similar to that used for the construction of classprobability models. That is, using background definitions of linear or kernel regression: 1. Construct the “best” regression model for the largest pharmacophore possible for the inhibitors. The probability of molecule M being “active” is P if: M has a hydrogen donor at position P1 and M has a hydrogen donor at position P2 and M has a negatively charged atom at position P3 and M has a negatively charged atom at position P4 and the distance between P1 and P2 is X1 and the distance between P1 and P3 is X2 and the distance between P1 and P4 is X3 and the distance between P2 and P3 is X4 and the distance between P2 and P4 is X5 and the distance between P3 and P4 is X6 and P is the kernel density estimate using instances I Where I is: P X1 X2 X3 X4 X5 X6 0.0 5.6 2.2 4.3 4.5 2.2 4.1 0.0 5.3 2.2 4.8 5.7 2.2 5.3 ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... Fig. 12. An example of a class-probability model obtained with the use of naive Bayes. The entries in the tabulation are the training instances returned by the search procedure in Fig. 7. 2. Estimate the performance of the model obtained in the previous step. The following details are relevant: (a) Models will again contain two sorts of constraints: first, those in the “prefix”, consisting of the functional groups and their pairwise distances (the pharmacophore per se); and second, the numeric constraint that predicts affinity using the distances between these functional groups. As before, we restrict models to be a single clause and the search is done in 2 steps:(i) the largest prefixes common to all inhibitors are found; and (ii) each prefix is extended to include the numeric constraint. We use the standard quadratic loss function (sum of squared differences between actual and predicted affinities) summed over all training instances. (b) For kernel regression, we use the Epanechnikov quadratic kernel [14] whose window size is determined by the kth-nearest neighbour (with k= 3). These choices are arbitrary, although the Epanechnikov kernel has some optimality properties that are described in [13]. (c) In Step (2) performance estimates will be estimated using a leave-one-out procedure. That is, each inhibitor is, in turn, categorised as a “test” instance. A regression model is constructed using the remaining instances and then used to predict the affinity of the test instance. The non-determinate nature of the domain will result in a set of predictions for the affinity. As with class-probability models, we have elected to select the highest of these as the final prediction (and as before, if the set of predictions is empty, then the test instance is not included in estimating the performance). Spearman’s rank correlation coefficient (rCV ) between actual and predicted values will be taken to be representative of the performance of the model. Results Figure 13 shows the correlation coefficients obtained by the leave-oneout procedure described. Also tabulated are the best results obtained with the standard QSAR 3-D technique CoMFA (comparative molecular field analysis) and a two-stage approach by King and colleagues ([21]). Method rCV Linear 0.34 Kernel 0.68 CoMFA 0.78 King 0.86 Fig. 13. Leave-one-out estimates of the rank correlation between actual and predicted values of affinity. “Linear” and “Kernel” stand for linear and kernel regression respectively. “CoMFA” represents the best QSAR model obtained with comparative molecular field analysis. “King” represents the approach described in [21]. The models obtained with linear and kernel regression are not as good as those with CoMFA or King. This is not surprising, given the models constructed here—single clauses with a regression constraint—are very simple (in contrast, King uses a kind of voting with multiple pharmacophore models). Nevertheless, the performance using linear regression is particularly poor, suggesting that the assumption of a linear model may be inappropriate. There are also good biological reasons to believe that this may be the case: assuming an ideal binding geometry, affinity values should decrease as distances deviate (on either side) from the ideal distances between functional groups. Quadratic regression should yield a better model under these circumstances: this is confirmed by an improved rCV value of 0.55. This is, of course, still well short of the mark achieved by King: some portion of the blame appears to lie in the selection rule used to obtain a final prediction (that is, the highest value of all predictions obtained). Better correlation values are obtainable, but the selection rule is not evident a priori. Example rules obtained with linear and kernel regression are shown in Figures 14 and 15. 7 Concluding Remarks The process of developing a new drug is long, laborious and expensive. A large part of the time and effort goes into the testing and assesment of compounds that The affinity of molecule M is A if: M has a hydrogen donor at position P1 and M has a hydrogen acceptor at position P2 and M has a hydrogen donor at position P3 and M has a negatively charged atom at position P4 and the distance between P1 and P2 is X1 and the distance between P1 and P3 is X2 and the distance between P1 and P4 is X3 and the distance between P2 and P3 is X4 and the distance between P2 and P4 is X5 and the distance between P3 and P4 is X6 and A = 1.06 + 0.55 X1 - 0.63 X2 + 0.29 X3 - 0.26 X4 + 0.06 X5 + 0.91 X6. Fig. 14. An example of a regression model obtained with the use of linear regression. The regression procedure computes the numbers in the equation for A. ultimately prove unsuitable as medications. Given a biological target to modulate, the role of chemistry is to identify the initial set of chemical compounds that can be taken forward for development. In this, structure-activity relationships (SARs) can play an important role in ensuring that a large proportion of the compounds proposed are also effective. In the past, the relational representation used by ILP has been repeatedly shown to be particularly well-suited to the task of constructing good SARs, with one caveat: the representation of activity has been a simple categorisation (usually “active” or “inactive”). In this paper, we have sought to extend ILP-constructed SARs to true quantitative models by utilising some standard statistical procedures as background knowledge. In doing so, we have had to confront some specific issues that arise during model construction and prediction when using statistical procedures within a first-order logic setting. To address these, we have proposed an ILP system that: (a) identifies statistical models in the SAR by employing the procedure in Section 4; and (b) uses the resulting SAR within domain-independent or domain-specific procedures to ensure deterministic prediction (as described in Section 5). The obvious limitation of the procedures proposed here is the lack of provable properties, in particular, about the optimality of models constructed. In the absence of such properties, we have attempted to demonstrate the practical utility of the procedures using using data on the inhibition of Thermolysin. To the best of our knowledge, the results represent the first examples of a constructing truly quantitative 3-dimensional SARs with ILP. The literature on attempting to incorporate statistical models within ILP hypotheses is relatively sparse. Both [18] and [36] are concerned with the use of linear regression by an ILP program (as a built-in definition in [18] and as a background predicate in [36]). In the former, the multiple instance problem is avoided by restricting background predicates that introduce the independent variables to be strictly deterministic (that is, functional). The latter simply ig- The affinity of molecule M is A if: M has a hydrogen donor at position P1 and M has a hydrogen donor at position P2 and M has a negatively charged atom at position P3 and M has a negatively charged atom at position P4 and the distance between P1 and P2 is X1 and the distance between P1 and P3 is X2 and the distance between P1 and P4 is X3 and the distance between P2 and P3 is X4 and the distance between P2 and P4 is X5 and the distance between P3 and P4 is X6 and A is the kernel regression estimate using instances I Where I is: A X1 X2 X3 X4 X5 X6 3.3 3.7 2.2 5.5 2.9 2.3 4.2 7.5 6.6 6.8 4.7 12.0 10.0 4.6 ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... Fig. 15. An example of a regression model obtained with the use of kernel regression. The entries in the tabulation are the training instances returned by the search procedure in Fig. 7. nores the multi-instance problem and treats all instances as independent data points. Adopting the same approach, Slattery and Craven [4] examine the use of a built-in naive Bayesian classifier. To this extent, the work here appears to the first to confront squarely the problems of using statistical procedures with non-determinate background knowledge. More generally, there is an increasing awareness within machine learning of the importance of learning statistical models from relational data. In Appendix B, we place this work within the context of broader efforts in the emerging field of statistical relational learning. The techniques presented are not confined to the particular statistical procedures used, a particular ILP system, or to the construction of SARs. This suggests three interesting ways in which the work here could be extended. First, the use of other conditional estimation procedures may yield better SAR models than those obtained here. Second, the same techniques could be used to extend the capabilities of other quantitative ILP approaches like those that construct first-order regression trees. Third, and most interestingly, the same procedures should provide any ILP system with the tools necessary to construct complex theories that combine first-order and ‘propositional’ aspects. An example of this are the theories described in [34]. A.D. Shapiro’s work on structured induction requires expert assistance to decompose hierarchically a complex induction task into sub-problems that can solved inductively (in Shapiro’s case, using a tree learner). While the approach was shown on difficult chess endgames to yield novel and comprehensible theories, the need to provide a complete decomposi- tion of the task has remained a principal difficulty with the technique. It is of interest to see the extent to which an ILP system can construct such theories automatically (equipped, for example, with a tree learner as background knowledge and appropriate domain-specific constraints). Acknowledgements Much of this work was done when the first author was at the Computing Laboratory, Oxford. Thanks are also due to Steve Moyle for several interesting discussions on the techniques described here and his generous computational support; to Nathalie Marchand-Geneste for the data on the Thermolysin inhibitors; and to Ravi Kothari for suggesting the use of kernel regression. References 1. R.A. Amar, D.R. Dooly, S.A. Goldman, and Q. Zhang. Multiple-Instance Learning of Real-Valued Data. In C.E. Brodley and A.P. Danyluk, editors, Proceedings of the Eighteenth International Conference on Machine Learning (ICML 2001), pages 3–10, San Francisco, 2001. Morgan Kaufmann Publishers. 2. W. Buntine. Operations for learning with graphical models. Journal of Artificial Intelligence Research, 2:159–225, 1994. 3. V. S. Costa, D. Page, M. Qazi, and J. Cussens. CLP(BN): Constraint Logic Programming for Probabilistic Knowledge. In Proceedings of the Nineteenth International Conference on Uncertainty in AI (UAI-2003), San Francisco, CA, 2003. Morgan Kaufmann. (to appear). 4. M. Craven and S. Slattery. Relational Learning with Statistical Predicate Invention: Better Models for Hypertext. Machine Learning, 43(1-2):97–119, 2001. 5. J. Cussens. Stochastic logic programs: Sampling, inference and applications. In Proceedings of the Sixteenth Annual Conference on Uncertainty in Artificial Intelligence (UAI-2000), pages 115–122, San Francisco, CA, 2000. Morgan Kaufmann. 6. T. Dietterich, R. Lathrop, and T. Lorano-Perez. Solving the multiple-instance problem with axis-parallel rectanges. Artificial Intelligence, 89:31–71, 1997. 7. D.R. Dooly, S.A. Goldman, and S.S. Kwek. Real-Valued Multiple-Instance Learning with Queries. In N. Abe, R. Khardon, and T. Zeugmann, editors, Proceedings of the Twelfth Conference on Algorithmic Learning Theory (ALT 2001), volume 2225 of LNAI, pages 167–180, Berlin, 2001. Springer-Verlag. 8. P. Finn, S. Muggleton, D. Page, and A. Srinivasan. Pharmacophore Discovery using the Inductive Logic Programming system Progol. Machine Learning, 30:241–270, 1998. 9. P.W. Finn. Computer-based screening of compound databases for the identification of novel leads. Drug Design Today, 1(9):363–370, 1996. 10. N. Friedman, L. Getoor, D. Koller, and A. Pfeffer. Learning probabilistic relational models. In Proceedings of the 16th International Joint Conference on Artificial Intelligence. Stockholm, Sweden, 1999. 11. L. Getoor, N. Friedman, D. Koller, and A. Pfeffer. Learning probabilistic relational models. In S. Dzeroski and N. Lavrac, editors, Relational Data Mining, pages 307– 335. Springer, Berlin, 2001.