Interpretative Neural Networks for QSAR

James H. Wikel*, Ernst R. Dow, and Michael Heathman

Abstract

A neural network program AUTONET was used to develop multiple neural network models trained to local minima. The network models were easily interpreted by visualization of the hidden weight coefficients.

Introduction

The role of the medicinal chemist has remained essentially unchanged for the past 50 years. Their role is dictated by their quest for rapid and efficient methods that optimize biological activity through structural variations. This need is historically driven by the fact that, on average, approximately 10,000 compounds are prepared and evaluated for every 1 that becomes a marketable drug. According to Lehman Brothers, financial analysts estimate that the current cost of drug development is nearly $600 million [1]. Additionally, increased development time has shortened the useful patent life in which companies can recover their costs; one-third of which are estimated to occur in the lead generation, discovery, and optimization phase. The dramatic increase in the cost of discovery resources is highlighted by the fact that a single, traditionally synthesized compound is estimated to cost $6,000 per research size sample [2]. These factors have contributed to a paradigm shift in the way pharmaceutical research is being conducted; companies are adopting approaches which reduce costs early in the drug development process.

Combinatorial chemistry and high through-put screening (HTS) are exciting techniques that are being adopted by the pharmaceutical and agrochemical industries in an effort to reduce costs and shorten discovery and optimization time. Computational scientists are contributing to this effort through combinatorial chemistry library analysis, diversity analysis, and quantitative structure activity relationship (QSAR) studies. QSAR studies rely heavily upon statistics to derive mathematical models which relate the biological activity of a series of compounds to one or more properties of the molecules. These properties, or descriptors, may be derived from numerous sources including refractive index, octanol/water partition coefficient or spectral data. In cases where experimental values for these properties are not available, several programs, including the popular CLOGP program, can be used for the computation of octanol/water partition coefficients [3]. Alternatively, theoretical properties also may be obtained from computational programs such as MOPAC [4]. A plethora of graph theory-based topological descriptors are available from programs such as MOLCONN-X [5]. Extensive lists of substituent parameters describing electronic (sigma), lipophilic (pi), and steric properties (cMR and Taft coefficients) are also available [6]. The initial phase of a QSAR study requires the collection of many of these descriptors prior to model building.

The seminal work in the field of QSAR was report by Hansch et. al. who demonstrated the use of regression analysis for model building [7]. In the intervening years since Hansch introduced regression analysis to chemistry, other methods have been developed and explored to circumvent some of the problems associated with this technique.

The success of regression analysis in QSAR model building depends upon an assumed linear relationship between the biological activity and one or more descriptors. As the number of descriptors increase however, regression analysis becomes problematic. One problem likely to occur in large descriptor sets, for example, is redundancy in information when descriptors are correlated. Latent variable techniques have become accepted methods of addressing this issue [8]. These techniques include the use of principal components in regression analysis and the method of partial least squares. A second problem encountered in using regression analysis is the a priori assumption of a model form (i.e. quadratic, cubic, use of cross terms, etc.). In order to address this issue, variable selection techniques such as stepwise forward and stepwise backward multiple linear regression analysis (MLR) were introduced. One recurrent problem in all of these methods is the fact that by using computational methods to generate descriptors, a modern dataset may contain more descriptors than compounds, that is, more columns of parameters than rows of compounds. This results in the introduction of the insidious problem described by Topliss and Edwards [9] - that correlations observed may be chance correlations.

Intriguing approaches using machine learning methods have been under study in the field of chemistry for the past decade. The first description of a simple neural network was provided in 1943 [10]. Interest in neural networks was slow until the 1980s when new computer architecture and learning algorithms began to appear. The use of artificial neural networks (ANN) in all fields has since grown substantially. In 1988, Hoskins et. al. reported the first use of process control in chemistry [11]. This was followed by two reports using ANN for prediction of protein secondary structure [12]. The use of ANN in chemistry has further expanded into the analysis of spectral data, pharmaceutical product development, classification of anticancer compounds, prediction of chemical reactivity, physical properties, electrostatic potential, ionization potentials as well as QSARs [13].

Neural networks are part of a new era of evolving computer technology in which a computer system has been designed to learn from data in a manner emulating the learning pattern in the brain [14]. Neural networks are typically used when there are a large number of observations and when the problem is not understood well enough to write a procedural program or expert system. Using neural networks, the solution to the problem is sought as follows:

1. an answer is calculated by multiplying each input by the connection weight;
2. products are summed at each hidden unit where a non-linear transfer function is applied; and
3. the output of each hidden unit is then multiplied by the connection weight and summed and interpreted.

The neural network "learns" by repeatedly passing through the data and adjusting its connection weights to minimize the error; in this case, the predicted, versus the actual biological activity. A neural network is thus a mathematical model to describe a non-linear hypersurface. The increasing interest and availability of neural network software has prompted several groups to apply this technology in QSAR studies. Aoyama reported the application of neural networks as a substitute for discriminant analysis [15]. Aoyama et. al. and subsequently Andrea et. al. applied neural networks in QSAR in a manner similar to multiple regression analysis [16]. So and Richards reported a comparative study of neural networks and regression analysis using a set of dihydrofolate reductase inhibitors [17]. Their results indicated neural networks were superior to regression analysis in providing accurate predictions, but that the design of the neural net was critical to obtaining these results. Reports by Manallack et. al. furthers solidifies the importance of the number of hidden units in the neural network on the predictive value of the derived model [18]. Song et. al. described an empirical based rule for dynamically adjusting the learning rate parameter during training [19]. This rule was developed following an extensive study for the optimization of numerous neural network parameters.

All of these reports may be viewed as successes for the application of neural networks in chemistry. However, their successes also highlighted four important factors which must be considered when using neural networks. First the design of the network is critical with respect to the number of hidden units involved. The network will overfit or memorize the data if too many hidden units are used. Conversely, the network will fail to generalize and become unstable if too few hidden units are used. The second factor which must be considered is the length of the training time. It is possible that networks may be overtrained, and thus destabilized, through the use of excessive training periods. Third, the selection of an appropriate test set and training set are important. The test set should adequately represent the entire dataset and be sufficiently large in order to properly train the neural network. The test set should also be contained within the neural network model. In addition it should be large enough to provide for an assessment of the model. Finally, the results obtained from neural networks can be difficult to interpret and apply to the drug design problem. This issue is especially troublesome for the medicinal chemist who is not an expert in the use and interpretation of neural network technology.

We have reported our initial results in the use of neural networks to identify the descriptors most relevant to biological activity [20]. The present report describes our continued work in this area and the enhancements we have made in the methodology. Our objective is to provide additional tools for the medicinal chemist in the area of molecular design. Our focus is to apply neural network technology early in the development of the SAR in a manner that, for the medicinal chemist, is easy to use. A further goal is to provide a technique whose results are both relevant and interpretable. Methodology has been developed and incorporated within a program, named AUTONET, that represents a self training neural network. Results from the neural network are presented visually in order to rapidly and easily convey to the medicinal chemist the important features derived by the neural network.

Methods

A detailed description of the theory behind a neural network has been adequately described elsewhere (see for example, Andrea [16b], So [17] or Manallack [21]). We will provide only a brief description of points most relevant here.

An artificial neural network consists of a number of "neurons" or "hidden units" that receive data from the outside, process the data, and output a signal. A "neuron" is essentially a regression equation with a non-linear output. When more than one of these neurons is used, non-linear models can be fitted. These networks have been shown to work well for modeling a number of different problems, including QSAR. Neural networks are known for their ability to model a wide set of functions without knowing the model a priori. The back propagation network receives a set of inputs which are multiplied by each neuron's weights (Figure 1). These products are summed for each neuron and a non-linear transfer function is applied. The bias has the effect of shifting the transfer function to the left or right. The transformed sums are then multiplied by the output weights where they are summed a final time, transformed, and interpreted. Since a back-propagation network is a supervised method, the desired output must be known for each input vector so an error (the difference between the desired output and the network's predicted output) can be calculated. This error is propagated backwards through the network (thus the name), adjusting the weights so that the next time the network sees the same input pattern, it will come closer to the desired output. The patterns are shown many times until the network either learns the relation or determines that there is none.

Figure 1. Sample neural network showing connections

for 2 inputs and 3 hidden units.

For our purposes, the input vector and output values were normalized between 0.1 and 0.9 by column. This ensures that no exceptionally large valued descriptors will have an undue effect on the network. All of the connection weights are initialized to very small random numbers (+/- 0.0005). This is necessary so each hidden unit will respond to a slightly different feature in the input vector. Each hidden unit outputs the hyperbolic tangent of the sum of the products of the inputs and the weights (Equation 1).

The hyperbolic tangent function has a range of -1 to 1, with the highest gain near 0. This compresses the output of the unit which defines a maximum contribution for each hidden unit. The output unit takes the sum of the products of the hidden units and the weights (Equation 2) and applies the transfer function (Equation 3).

x is the value of the output_unit. This function has minimum and maximum values of 0 and 1 respectively. Since the output values were normalized between 0.1 and 0.9, this allows the network to slightly exceed the minimum and maximum values that were given in the original data file.

Once the output is calculated, it is compared to the desired output value for that particular vector (the biological activity). An error, , is calculated according to Equation 4 and is used in a gradient descent algorithm to adjust the weights of the network (Equation 5).

= experimental - predicted (4)

where is the learning rate which controls the step size of the gradient descent algorithm. The learning rate is typically between 0 and 1 and is decreased during training as the solution is reached. The term (out (1-out)) is the derivative of the transfer function (Equation 3). The hidden unit weights are adjusted in a similar manner (Equation 6). The term (1 + out)(1 - out) is the derivative of the hyperbolic tangent transfer function.

Each input vector and desired output pair for the entire training set was presented to the network and the weights were adjusted. The training set was generated by sorting all the data based on the output, biological activity, and then every fourth compound was placed in a testing set and the remaining compounds were used for the training set. Sorting the small datasets that are typical of QSAR studies ensured that the test set was as representative as possible. One complete cycle through the data is called an epoch. During each epoch, the order in which the compounds were presented was randomized. This procedure improved the overall performance of the neural network. The training and testing errors were calculated (Equation 7) for the testing set every 10 epochs and this value was saved with the network weights.

If the testing error had not decreased in 250 epochs, the network was returned to that set of weights. Often this technique will produce a large number of networks of different sizes that have similar training and testing errors. The best network is the one with the smallest testing error. The r² value is also checked since it is possible to have a small error and a poor r²'d value. If there are several networks of similar errors, the smallest network is often the easiest to interpret.

The model knowledge in the neural network can be discerned by examining the weights. As the weight from an input descriptor to a hidden unit approaches zero, then the effect that the descriptor can have on the model approaches zero. However, it was not obvious which descriptors were contributing to the model and which only had chance effects since every descriptor had a weight coefficient. In order to make a cutoff criteria, three random descriptors were included as input vectors. The networks were trained with these random descriptors along with the other descriptors. Now it was possible to compare the various chemical descriptors against the random descriptors to determine which descriptors were more significant than random noise.

A second set of networks was trained with a reduced set of descriptors. Only those descriptors whose absolute value of the weight coefficient were larger than the largest of the random descriptors were used. In addition the descriptors from the best (lowest testing error) network were automatically included. Often, the networks using this "reduced" set of descriptors will outperform the original set. This increases the likelihood that a difficult model can be solved and also that an easy-to-interpret network will be constructed. Non-linear effects are also determined. An examination of the weights for a descriptor where the largest weight and the second largest are of opposite sign and are at least half the magnitude of the largest weight in the network suggests the presence of a non-linear effect. This tends to identify non-linear effects with a fewer number of compounds than otherwise possible.

In order to present the chemist with useful information from the neural network, certain data are visualized. The hidden unit weights for each descriptor for each network are displayed in a color map. A green color indicates that a weight value is near zero, blue is a negative weight, and red is positive. If an output weight is negative, all of the weights entering that particular hidden unit are multiplied by (-1). This can be done since the hyperbolic tangent function is symmetric about the y axis. If a descriptor has all red weights, then increasing the value of that descriptor will have a positive effect on the output of the network. The chemist thus can quickly see which descriptors are consistently having an effect on the different models. The numerical value of the hidden unit weights are also available. The chemist can also examine the weights of the best network by testing error, or if there are several that are close, the smallest hidden unit size. Each network is given a brightness which is proportional to its testing error. Random descriptors are given colors as well. The chemist may opt to use these descriptors in a multiple regression study.

All networks were of the back-propagation type and trained on a Silicon Graphics Workstation. The program AUTONET is written in the language C. The networks were trained to predict activity. A hyperbolic tangent transfer function was used with a user definable learning rate coefficient between 0.1 and 1.0. All inputs were normalized between 0.1 and 0.9. Multiple networks, at least 3, were developed for each level of hidden units specified. The number of hidden units were predetermined at discrete levels of 1,3,5,7,9,11 and 21 depending on the number of compounds in the dataset. Presentation of inputs to the ANN was randomized after every epoch. Each network was also started from a randomized input in order to start each one at a different point on the response surface. Results from all of the networks were compared. The total number of networks built is two times the product of the number of passes and the number of hidden units. This is the combined total of networks constructed with the complete set of descriptors and the number of networks built with the reduced set of descriptors. For example, 30 networks are built for each training/test set with 3 passes and 5 hidden units.

Results and Discussion

In a typical neural network application, the dataset is randomly divided into two subsets. One group, the larger of the two, is used to train the network while the smaller subset is used to evaluate the predictive power of the network. QSAR datasets are typically small in the early stages of the project and thus it becomes impractical to reduce them substantially. In the present study, the datasets were sorted by activity and exemplars were removed for training purposes by one of two methods. The first method requires the removal of compounds from the sorted dataset at predetermined intervals. The second method, used for smaller datasets, is a leave-one-out procedure requiring the removal of each compound, one at a time, to serve as the test case. In each case the observation(s) that was removed served as the test case for the network. If the datasets were not sorted before these techniques were applied, the learning power of the resulting neural network was compromised.

Published data were used for this study. Dataset 1, reported by Selwood et. al., consisted of 31 compounds and 53 descriptors [22]. The second dataset, reported by Dunn et. al., consisted of 13 compounds and 5 descriptors [23]. We added 58 additional descriptors to this dataset in order to compound the problem for the neural network[24]. These two datasets were selected to represent a QSAR study in the early phase (13 compounds) and in a more developed phase (31 compounds) of progress. Since it is not possible to represent ortho- or meta-substituted compounds appearing in the dataset from Dunn as words to the network, a binary representation was used: ortho 0 1, meta 1 0. Missing values for individual atom descriptors, which were due to differing numbers of atoms, were assigned a value of zero. The third dataset, reported by Howbert et. al. consisted of 47 compounds and 7 descriptors [25]. In this dataset, compounds were classified as active or inactive based upon their in vivo potency as antitumor agents.

In order to determine if our method would allow the neural networks to merely memorize data, even random noise, we conducted the following experiment. Datasets were created containing 20 "compounds" and 10, 20, 40, or 80 descriptors of random numbers. The output representing biological activity was a random number. A network with 13 descriptors (10 descriptors plus 3 random noise) and 1 hidden unit has 14 adjustable parameters. The same dataset with 21 hidden units has 274 adjustable parameters and thus a definite potential to overfit the data exists. With the current methodology using the leave-one-out protocol, all the networks that were created consistently memorized the data as evidenced by three characteristics: 1) a low training error on the learning set and a high testing error on the test set; 2) generally (>75%) less than 100 epochs in each network; and 3) the hidden weight coefficients, when examined in the color maps, were similar to known random noise descriptors (Color Plate 1).

Dataset 1 (Selwood Dataset). A critical step prior to the construction of a neural network is the selection of the appropriate size of the training set and the test set. Care must be taken to ensure that each set is representative of the other. The training set should be as large as possible in order to provide the neural network with the best opportunity to learn. The test set should be sufficiently large to provide new cases in order to fairly evaluate the neural network. The Selwood dataset was sorted according to the output response and every fourth record was removed. This created a test set of 31/4 = 7 records (records 4, 8, 12, 16, 20, 24, 28) and a training set of 24 records. Training sets were randomized before the first epoch and before each subsequent epoch. Three descriptors of random numbers were added, the learning rate was set at 0.6, and a total of 15 networks were generated with all of the descriptors. Three networks were generated for each hidden unit level of 1, 3, 5, 7, and 9. The results are depicted in Color Plate 2 using 3 colors (positive weight coefficients in red and negative weight coefficients in blue). The left side of the panel represents the networks that include all of the descriptors (one descriptor per row and one network per column). Each network is given a brightness which is proportional to its testing error and thus the darker appearing columns are those networks with the largest testing error. A few networks failed to find any descriptor more important than another; as evidenced by a green vertical bar. However, the remaining networks had low testing errors and found one or more descriptors to be consistently important to the learning behavior. The right side of the panel represents 15 new networks using only those descriptors that frequently provided weights greater than any of the three random number descriptors. The results indicate that the learning behavior of the ANNs were most often related to several descriptors. The descriptors ATCH4, ESDL3 and CLOGP were identified as important to the learning of the ANN and had positive weight coefficients while ATCH6, DIPV_X, DIPV_Z, NSDL1, and NSDL7 were important but had negative weight coefficients.

Regardless of the learning behavior of the individual neural network, useful information is readily conveyed by the color panels. Notably the display illustrates that not all network configurations behave similarly. Some networks failed to learn; as evidenced by the vertical green columns. The networks with the smaller number of hidden units, 1 and 3 located at the left edge of the panel, appear to have the best learning results. The r² training and r² testing in the output file suggests that these networks were predictive; although this is not required of the AUTONET derived neural network nor is it the focus of our interest. Due to our deliberate undertraining of the neural network, it is probable that this method found local minima which may be responsible for the results. Since each network starts from a different point on the response surface, these minima may or may not be the same. The actual predictive power of ANN derived from local minima may be suspect. The critical observation is that the solutions to these local minima are derived from a similar set of descriptors. The interpretation of the learning behavior of the network is possible due to the commonality of the solutions (e.g., the important descriptors for learning) to these local minima from the entire collection of multiple networks.

As reported previously multiple regression analysis may then applied to the dataset using the most important descriptors identified by the neural network[20]. The Selwood dataset has been the subject of analysis by numerous approaches[26]. These studies illustrate a general point in model building that many models exist to explain a dataset. The purpose of the exercise is to find reasonable models upon which to base additional experiments.

Dataset 2 (Dunn Dataset). This was a small dataset of only 13 compounds and 5 descriptors as originally reported. We added 58 descriptors in order to better represent a realistic situation in which no descriptor bias is assumed. Since this dataset was small we used a true cross validation technique to train the neural networks in which every compound is removed once to serve as the training set. For each network 3 descriptors of random data were added and the learning rate set equal to 0.6. The results for these networks are presented in Color Plates 3. The color scheme for these plates is the same as described for Dataset 1. The results are more complex than the Selwood results as more network are produced. The results following the removal of one compound to serve as the test case are depicted in a set of 15 columns representing the 15 networks (3 passes x 5 sets of hidden units). The results in the panel are those from a total of 195 networks. The interpretation is further complicated by the fact that the first and last columns are black indicating the network was not able to train properly when the observations at either extreme are removed to serve as the test cases. This behavior is indicative of a model from which extrapolations are not possible. Similarly, darker columns within the body of the color plate indicate those compounds whose removal produced an less satisfactory network training. However useful information becomes apparent from the color plates. The sigma descriptor and it component Swain-Lupton R descriptor appear important and have a negative weight coefficient. Thus one may conclude that the electronic effects of the substituent may be important to the biological activity of these compounds. This in agreement with the results from regression analysis reported[23] in previous studies.

Dataset 3 (Howbert Dataset). The largest dataset studied contained 47 compounds. Unlike the previous examples, the dataset represents a classification problem (active/inactive) based upon the in vivo biological potency. Using the same methodology as described, the neural networks were trained to identify features for the correct classification of the compound. Active compounds were designated with the number 1 and inactive compounds were designated with the number 0 according to the definition provided by Howbert et.al[25]. Three neural networks were trained at each of the hidden unit levels (1,3,5,7,9,11,21) following the addition of 3 descriptors of random numbers with a learning rate equal to 0.6.

It is apparent from the results depicted in Color Plate 4 the neural networks strongly identified the VDWVOL, a negative weight coefficient, and pi, a positive weight coefficient, as descriptors important in the correct classification of the compounds into active or inactive groups. Both of these descriptors were identified by Howbert et. al. using cluster significance analysis of their data. Multiple regression analysis did not yield a statistically valid model. Examination of the training and testing errors indicate, as expected, the neural network was not predictive.

Conclusion

The applications of neural networks have required large datasets and have shown extensive training periods in order to achieve a predictive solution. We have show using the techniques with the neural network program AUTONET that it is possible to extract information using a neural network from relative small datasets and from networks that are not statistically predictive. The AUTONET method uses a series of multiple, short training neural networks to provide local minima as solutions. The information content is extracted from the coefficients of the hidden weights associated with the input descriptors with the overall solution provided by a consensus of solutions to the local minima. We found that in spite of short training periods the neural network memorized random data very quickly. However a characteristic profile of an overtrained network was identified as: 1) a low training error on the learning set and a high testing error on the test set; 2) generally (>75%) less than 100 epochs in each network; and 3) the hidden weight coefficients similar to each other and similar to the coefficients of know random noise descriptors. The addition of three descriptors of random numbers allow for the establishment of a level from which to judge background noise and chance correlation. The focus of the AUTONET approach is not to achieve a predictive solution. We have attempted to broaden the scope of the utility of neural networks in QSAR by gaining information in the absence of a predictive solution. An often encountered difficulty with neural networks is their lack of interpretation. The AUTONET approach addresses this through the visual display of the hidden unit weights and thus rapidly conveys useful and informative results to the user.

References

1. SCRIP World Pharmaceutical News, January 9, 1996, PJB Publications, Ltd.

2. Huston, S. Integrated Strategies in Drug Discovery. 1995, 23, 19-21.

3. Daylight Chemical Information Systems, Inc. Sante Fe, NM USA.

4. Stewart, J. J. P. MOPAC, QCPE Bull., 1989, 9, 10. QCPE Program 255, MOPAC

5. MOLCONN-X, Hall Associates Consulting, Quincy, MA. USA 02170

6. a) Hansch, C.; Leo, A. Substituents Constants for Correlation Analysis in Chemistry and Biology. John Wiley & Sons, New York, 1979. b) Boyd, D.; Seward, C. M. QSAR: Rational Approaches to the Design of Bioactive Compounds.; Silipo, C.; Vittoria, A., Ed.; Elsevier Science Publishers B. V.: Amsterdam, 1991; 167-170.

7. a) Hansch, C. ; Muir, R. M.; Fujita, T.; Maloney, P.P.; Geiger, F.; Streich, M. The Correlation of Biological Activity of Plant Growth Regulators and Chloromycetin Derivatives with Hammett Constants and Partition Coefficients. J. Am. Chem. Soc. 1963, 85, 2817-2824.

8. a) Rawlings, J. O. Applied Regression Analysis: A Research Tool. Wadsworth & Brooks, 1988. b) Dunn, W. J.; Wold, S.; Edlund, U.; Hellberg, S.; Gasteiger, J. Quant. Struct.-Act. Relat., 1984, 3, 131-137.

9. Topliss, J. G.; Edwards, R. P. Chance Factors in Studies of Quantitative Structure -Activity Relationships. J. Med. Chem., 1979, 22, 1238- 1244.

10. McCulloch, W. S.; Pitts, W. A Logical Calculus of the Ideas Immanent in Nervous Activity. Bull. of Math. Bio., 1943, 5, 115-133.

11. Hoskins, J. C.; Himmelbau, D. M. Artificial Neural Network Models of Knowledge Representation in Chemical Engineering. Comput. Chem. Eng., 1988, 12, 881-890.

12. a) Qian, N.; Sejnowski, t. J. Predicting the Secondary Structure of Globular Proteins Using Neural Network Models. J. Mol. Biol., 1988, 202, 865-884. b) Bohr, H.; Bohr, J.; Brunak, S.; Cotterill, R.; Lautrup, B.; Norskov, L.; Olsen, O.; Petersen, S. Protein Secondary Structure and Homology by Neural Networks FEBS Lett., 1988,241, 223-228.

13. See references cited in Sigman, M. E.; Rives, S. S. Prediction of Atomic Ionization Potentials I-III Using an Artificial Neural Network. J. Chem. Inf. Comput. Sci., 1994, 34, 617-620.

14. Rumelhart, D. B. Parallel Distributed Processing, Feldman, J. A.; Hayes, P. J.; Rumelhart, D. B., Ed.; The MIT Press, London, 1982, 1, 318-363.

15. Aoyama, T.; Suzuki, Y.; Ichikawa, H. Neural Networks Applied to Structure- Activity Relationships. J. Med. Chem., 1990, 33, 905-908.

16. a) Aoyama, T.; Suzuki, Y.; Ichikawa, H. Neural Networks Applied to Quantitative Structure-Activity Relationship Analysis. J. Med. Chem., 1990, 33, 2583-2590. b) Andrea, T. A.; Kalayeh, H. Application of Neural Networks J. Med. Chem., 1991, 34, 2824-2836.

17. So, S.; Richards, W. G. Application of Neural Networks: quantitative Structure-Activity Relationships of the derivatives of 2,4-Diamino-5- (substituted-benzyl)pyrimidines as DHFR Inhibitors. J. Med. Chem., 1992, 35, 3201-3207.

18. Manallack, D. T..; Ellis, D. D.; Livingstone, D. J. Analysis of Linear and Nonlinear QSAR Data Using Neural Networks. J. Med. Chem., 1994, 34, 3758-3767.

19. Song, X.; Chen, Z.; Yu, R. Artificial Neural Networks Applied to Odor Classification for Chemical Compounds. Computers Chem., 1993, 17, 303-308. b) Gakh, A. A.; Gakh, E. R.; Sumpter, B. G.; Nord, D. W. Neural Network-Graph Theory Approach to the Prediction of the Physical Properties of Organic Compounds. J. Chem. Inf. Comput. Sci., 1994, 34, 832- 839.

20. Wikel, J. H.; Dow, E. R. The Use of Neural Networks for Variable Selection in QSAR. Bioorg. Med. Chem. Lett., 1993, 3, 645-651.

21. Manallack, D. T.; Livingstone, D. J. Neural Networks and Expert Systems in Molecular Design. Methods Princ. Med. Chem. 1995, 3(Advanced Computer- Assisted Techniques in Drug Discovery), 293-318

22. Selwood, D. L.; Livingstone, D. J.; Comley, J. C.; O'Dowd, A. B.; Hudson, A. T.; Jackson, P.; Jandu, K. S.; Rose, V. S.; Stables, J. N. Structure-Activity Relationships of Antifilarial Antimycin Analogues: A Multivariate Pattern Recognition Study J. Med. Chem., 1990, 33, 136-142.

23. Dunn, W. J.; Greenberg, M. J.; Callejas, S. S. Use of Cluster Analysis in the Development of Structure-Activity Relations for Antitumor Triazenes. J. Med. Chem., 1976, 19, 1299-1301.

24. The descriptors added to the original dataset were the sterimol term, L1, B1 through B4, (Verloop, A.; Hoogenstraaten, W.; Tipker, J. Drug Design; Ariens, E. J., Ed.; Academic Press: New York, 1976; Vol 7, pp. 165-207.) and molecular connectivity indices ( Kier, L. B.; Hall, L. H. Molecular Connectivity Analysis; Bawden, D., Ed.; Research Studies Press Ltd.: England, 1986) obtained from the program MOLCONN-X, Hall Associates Consulting, Quincy, MA 02170.

25. Howbert, J. J.;Grossman, C. S.; Crowell, T. A.;Rieder, B. J.; Harper, R. W.; Kramer, K. E.; Tao, E. V.; Aikins, J.; Poore, G. A.; Rinzel, S. M.; Grindey, G. B.; Shaw, W. N.; Todd, G. C. Novel Agents Effective against Solid Tumors: The Diarylsulfonylureas. Synthesis, Activities, and Analysis of Quantitative Structure-Activity Relationships. J. Med. Chem., 1990, 33, 2393- 2407.

26. a) McFarland, J. W.; Gans, D. J. On Identifying Likely Determinants of Biological Activity in High Dimensional QSAR Problems. Quant. Struct.-Act. Relat., 1994, 13, 11-17. b) Rogers, D.; Hopfinger, A. J. Application of Genetic Function Approximation to Quantitative Structure-Activity Relationships and Quantitative Structure-Property Relationships J. Chem. Inf. Comput. Sci., 1994, 34, 854-866. c) Kubinyi, H. Variable Selection in QSAR Studies. II. A Highly Efficient Combination of Systematic Search and EvolutionQuant. Struct.-Act. Relat., 1994, 13, 393- 401.

An Introduction to QSAR Methodology

Allen B. Richon

Introduction

Drug design is an iterative process which begins with a compound that displays an interesting biological profile and ends with optimizing both the activity profile for the molecule and its chemical synthesis. The process is initiated when the chemist conceives a hypothesis which relates the chemical features of the molecule (or series of molecules) to the biological activity. Without a detailed understanding of the biochemical process(es) responsible for activity, the hypothesis generally is refined by examining structural similarities and differences for active and inactive molecules. Compounds are selected for synthesis which maximize the presence of functional groups or features believed to be responsible for activity.

The combinatorial possibilities of this strategy for even simple systems can be explosive. As an example, the number of compounds required for synthesis in order to place 10 substituents on the four open positions of an asymmetrically disubstituted benzene ring system is approximately 10,000. The alternative to this labor intensive approach to compound optimization is to develop a theory that quantitatively relates variations in biological activity to changes in molecular descriptors which can easily be obtained for each compound. A Quantitative Structure Activity Relationship (QSAR) can then be utilized to help guide chemical synthesis. This chapter develops the concepts used to derive a QSAR and reviews the application of these techniques to medicinal research.

Statistical Concepts

Computational chemistry represents molecular structures as a numerical models and simulates their behavior with the equations of quantum and classical physics. Available programs enable scientists to easily generate and present molecular data including geometries, energies and associated properties (electronic, spectroscopic and bulk). The usual paradigm for displaying and manipulating these data is a table in which compounds are defined by individual rows and molecular properties (or descriptors) are defined by the associated columns. A QSAR attempts to find consistent relationships between the variations in the values of molecular properties and the biological activity for a series of compounds so that these "rules" can be used to evaluate new chemical entities.

A QSAR generally takes the form of a linear equation

Biological Activity = Const + (C₁ P₁) + (C₂ P₂) + (C₃ P₃) + ...

where the parameters P₁ through P_n are computed for each molecule in the series and the coefficients C₁through C_n are calculated by fitting variations in the parameters and the biological activity. Since these relationships are generally discovered through the application of statistical techniques, a brief introduction to the principles behind the derivation of a QSAR follows.

The work reported from The Sandoz Institute for Medical Research on the development of novel analgesic agents¹ can be used as an example of a simple QSAR. In this study, vanillylamides and vanillylthioureas related to capsaicin were prepared and their activity was tested in an in vitro assay which measured ⁴⁵Ca²⁺influx into dorsal root ganglia neurons. The data, which was reported as the EC₅₀ (µM), is shown in Table 1 (note that compound 6f is the most active of the series).

TABLE 1
Capsaicin Analogs Activity Data

Cmpd Number	Cmpd Name	X	EC50(µM)
1	6a	H	11.80 ± 1.90
2	6b	Cl	1.24 ± 0.11
3	6d	NO2	4.58 ± 0.29
4	6e	CN	26.50 ± 5.87
5	6f	C6H5	0.24 ± 0.30
6	6g	N(CH3)2	4.39 ± 0.67
7	6h	I	0.35 ± 0.05
8	6i	NHCHO	???

In the absence of additional information, the only way to derive a best "guess" for the activity of 6i is to calculate the average of the values for the current compounds in the series. The average, 7.24, provides a guess for the value of compound 8 but, how good is this guess? The graphical presentation of the data points is shown in Graph 1.

GRAPH 1
Capsaicin Analogs Activity Data

The standard deviation of the data, s, shows how far the activity values are spread about their average. This value provides an indication of the quality of the guess by showing the amount of variability inherent in the data. The standard deviation is calculated as shown below.

Rather than relying on this limited analysis, one would like to develop an understanding of the factors that influence activity within this series and use this understanding to predict activity for new compounds. In order to accomplish this objective, one needs:

• binding data measured with sufficient precision to distinguish between compounds;
  
  
• a set of parameters which can be easily obtained and which are likely to be related to receptor affinity;
  
  
• a method for detecting a relationship between the parameters and binding data (the QSAR) and
  
  
• a method for validating the QSAR.
  
  

The QSAR equation is a linear model which relates variations in biological activity to variations in the values of computed (or measured) properties for a series of molecules. For the method to work efficiently, the compounds selected to describe the "chemical space" of the experiments (the training set) should be diverse. In many synthesis campaigns, compounds are prepared which are structurally similar to the lead structure. Not surprisingly, the activity values for this series of compounds will frequently span a limited range as well. In these cases, additional compounds must be made and tested to fill out the training set.

The quality of any QSAR will only be as good as the quality of the data which is used to derive the model. Dose-response curves need to be smooth, contain enough points to assure accuracy and should span two or more orders of magnitude. Multiple readings for a given observation should be reproducible and have relatively smaller errors. The issue being addressed is the signal-to-noise ratio. The variation of the readings obtained by repeatedly testing the same compound should be much smaller than the variation over the series. In cases where the data collected from biological experiments do not follow these guidelines, other methods of data analysis should be utilized since the QSAR models derived from the data will be questionable.

Once biological data has been collected, it is often found that the data is expressed in terms which cannot be used in a QSAR analysis. Since QSAR is based on the relationship of free energy to equilibrium constants, the data for a QSAR study must be expressed in terms of the free energy changes that occur during the biological response. When examining the potency of a drug (the dosage required to produce a biological effect), the change in free energy can be calculated to be proportional to the inverse logarithm of the concentration of the compound.

G₀ = - 2.3RTlogK = log 1/[S]

Further, since biological data are generally found to be skewed, the log transformation moves the data to a nearly normal distribution. Thus, when measuring responses under equilibrium conditions, the most frequent transformation used is to express concentration values (such as IC₅₀, EC₅₀, etc.) as log[C] or log 1/[C]. The transformed data for the capsaicin agonists are shown in Table 2.

TABLE 2

Capsaicin Analogs

Transformed Data

Cmpd Number	Cmpd Name	X	EC50	Log EC50	Log 1/EC50
1	6a	H	11.80 ± 1.90	1.07	-1.07
2	6b	Cl	1.24 ± 0.11	0.09	-0.09
3	6d	NO2	4.58 ± 0.29	0.66	-0.66
4	6e	CN	26.50 ± 5.87	1.42	-1.42
5	6f	C6H5	0.24 ± 0.30	- 0.62	0.62
6	6g	N(CH3)2	4.39 ± 0.67	0.64	-0.64
7	6h	I	0.35 ± 0.05	- 0.46	0.46
8	6i	NHCHO	?? ± ??	??	??

The effect of this transformation on the spread of the data relative to the average is shown in Graph 2. Note that the data points, projected onto the Y-axis, have become more uniformly distributed.

GRAPH 2
Capsaicin Analogs
Transformed Data

Given the transformed data, our best guess for the activity of 6i is still the average of the data set (or 0.40). As before, the error associated with this guess is calculated as the square root of the average of the squares of the deviations from the average.

This is an example data set intended to show the general approach; real data sets would have many more compounds and descriptors. Since the purpose of a QSAR is to highlight relationships between activity and structural features, we would like to find one or more structural features which relate these molecules and their associated activity. Additionally, we would like to find a parameter that works consistently for all of the molecules in the series.

There are several potential classes of parameters used in QSAR studies. Substituent constants and other physico-chemical parameters (such as Hammett sigma constants) measure the electronic effects of a group on the molecule. Fragment counts are used to enumerate the presence of specific substructures. Other parameters can include topological descriptors and values derived from quantum chemical calculations.

The selection of parameters is an important first step in any QSAR study. If the association between the parameter(s) selected and activity is strong, then activity predictions will be possible. If there is only weak association, knowing the value of the parameter(s) will not help in predicting activity. Thus, for a given study, parameters should be selected which are relevant to the activity for the series of molecules under investigation and these parameters should have values which are obtained in a consistent manner.

The Sandoz group divided their analysis of capsaicin analogs into three regions: the A-region which was occupied by an aromatic ring; the B-region which was defined by an amide bond; and the C-region which was occupied by a hydrophobic side-chain (See figure in Table 1). The hypothesis for the C-region assumed that a small, hydrophobic substituent would increase activity. Given this assumption, the parameters selected to best define this characteristic were molar refractivity (size) and , the hydrophobic substituent constant. These values are given in Table 3.

TABLE 3
Capsaicin Analogs
Parameter Values

Cmpd Number	Cmpd Name	X	Log EC50		MR
1	6a	H	1.07	0.00	1.03
2	6b	Cl	0.09	0.71	6.03
3	6d	NO2	0.66	- 0.28	7.36
4	6e	CN	1.42	- 0.57	6.33
5	6f	C6H5	- 0.62	1.96	25.36
6	6g	N(CH3)2	0.64	0.18	15.55
7	6h	I	- 0.46	1.12	13.94
8	6i	NHCHO	??	??	??

The data above can be analyzed for relationships by two means: graphically and statistically. The most visual approach to a problem with a limited number of variables is graphical. In this case, a plot of activity versus either molar refractivity or hydrophobicity gives some insight into the relationship between the parameters and activity. The plots derived by the Sandoz group are reproduced in Graph 3.

GRAPH 3
Capsaicin Analogs
Parameter Values

Does the graph provide insight into the the activity for compound 6i? Does knowing the value for either the hydrophobicity or molar refractivity parameters for this compound provide a good estimate for activity?

Since this is a simple example where only two values are examined, the answers to these questions are a qualified yes. In more complex situations however, where multiple parameters are correlated to activity, statistics is used to derive an equation which relates activity to the parameter set. The linear equation which defines the best model for this set of data is

Log EC₅₀ = 0.764 - (0.817)

How much confidence should we place in this model? The first step to answering this question is to determine how well the equation predicts activities for known compounds in the series. The equation above estimates the average value for the EC₅₀ based on the value for ; because assays vary, it is not surprising that individual values will differ from the regression estimate. The difference between the calculated values and the actual (or measured) values for each compound is termed the residual from the model. The calculated values for activity and their residuals (or the errors of the estimate for individual values) are shown in Table 4.

TABLE 4
Capsaicin Analogs
Calculated Values

Cmpd Number	Cmpd Name	X	Log EC50		Calculated Log EC50	Residual
1	6a	H	1.07	0.00	0.79	0.28
2	6b	Cl	0.09	0.71	0.21	- 0.12
3	6d	NO2	0.66	- 0.28	1.02	- 0.36
4	6e	CN	1.42	- 0.57	1.26	0.16
5	6f	C6H5	- 0.62	1.96	- 0.81	0.19
6	6g	N(CH3)2	0.64	0.18	0.65	- 0.01
7	6h	I	- 0.46	1.12	- 0.12	- 0.34
8	6i	NHCHO	??	- 0.98	1.60	??

The residuals are one way to quantify the error in the estimate for individual values calculated by the regression equation for this data set. The standard error for the residuals is calculated by taking the root-mean-square of the residuals (in this calculation, the denominator shown as decremented by two to reflect the estimation of two parameters).

In order to be an improved model, the standard deviation of the residuals calculated from the model should be smaller than the standard deviation of the original data. The standard error about the mean was previously calculated to be 0.76 whereas the standard error from the QSAR model is 0.28. Clearly, the the use of linear regression has improved the accuracy of our analysis. The plot of measured values versus calculated is shown in Graph 4 with a 45° line.

GRAPH 4
Capsaicin Analogs
Predicted Versus Actual EC₅₀ Values

There are several assumptions inherent in deriving a QSAR model for a series of compounds. First, it is assumed that parameters can be calculated (or measured in some cases) more accurately and cheaply than activity can be measured. Second, it is assumed that deviations from the best fit line follow a normal (Gaussian) distribution. Finally, it is assumed that any variation in the line described by the QSAR equation is independent of the magnitude of both the activity and the parameters. Given these assumptions, the quality of the model can be gauged using a variety of techniques.

Variation in the data is quantified by the correlation coefficient, r, which measures how closely the observed data tracks the fitted regression line. Errors in either the model or in the data will lead to a bad fit. This indicator of fit to the regression line is calculated as:

where the Regression Variance is defined as the Original Variance minus the Variance around the regression line. The Original Variance is the sum-of-the-squares distances of the original data from the mean. This can be viewed graphically as shown in Graph 5.

The calculation is carried out as follows:

Original Variance = (1.07 - 0.40)² + (0.09 - 0.40)² + ...

Original Variance = 3.49

Variance around the line = (0.28)² + (- 0.12)² + (- 0.36)² + ...

Variance around the line = 0.40

Regression Variance = Original Variance - Variance around the line

Regression Variance = 3.49 - 0.40 = 3.09

r² = Regression Variance/Original Variance

r² = 3.09/3.49

r² = 0.89

Possible values reported for r² fall between 0 and 1. An r² of 0 means that there is no relationship between activity and the parameter(s) selected for the study. An r² of 1 means there is perfect correlation. The interpretation of the r² value for the capsaicin analogs is that 89% of the variation in the value of the Log EC₅₀ is explained by variation in the value of , the hydrophobicity parameter.

GRAPH 5
Capsaicin Analogs
Derivation of r² values

While the fit of the data to the regression line is excellent, how can one decide if this correlation is based purely on chance? The higher the value for r² the less likely that the relationship is due to chance. If many explanatory variables are used in a regression equation, it is possible to get a good fit to the data due to the flexibility of the fitting process; a line will fit two points perfectly, a quadratic curve will fit three, multiple linear regression will fit the observed data if there are enough explanatory variables². Given the assumption that the data has a Gaussian distribution, the F statistic below assesses the statistical significance of the regression equation.

The F statistic is calculated from r² and the number of data points (or degrees of freedom) in the data set. The F ratio for the capsaicin analogs is calculated as:

This value often appears as standard output from statistical programs or it can be checked in statistical tables to determine the significance of the regression equation. In this case, the probability that there is no relationship between activity and the value is less than 1% (p=0.01).

We have found that hydrophobicity values correlate well with biological activity. Does the addition of a size parameter (MR) improve our model? In order to analyze a relationship which is possibly influenced by several variables (or properties), it is useful to assess the contribution of each variable. and MR appear to be somewhat correlated in this data set so the order of fitting can influence how much the second variable helps the first. Multiple linear regression is used to determine the relative importance of multiple variables to the overall fit of the data.

Multiple linear regression attempts to maximize the fit of the data to a regression equation (minimize the squared deviations from the regression equation) for the biological activity (maximize the r² value) by adjusting each of the available parameters up or down. Regression programs often approach this task in a stepwise fashion. That is, successive regression equations will be derived in which parameters will be either added or removed until the r² and s values are optimized. The magnitude of the coefficients derived in this manner indicate the relative contribution of the associated parameter to biological activity.

There are two important caveats in applying multiple regression analysis. The first is based on the fact that, given enough parameters any data set can be fitted to a regression line. The consequence of this is that regression analysis generally requires significantly more compounds than parameters; a useful rule of thumb is three to six times the number of parameters under consideration. The difficulty is that regression analysis is most effective for interpolation and it is extrapolation that is most useful in a synthesis campaign (i.e., the region of experimental space described by the regression analysis has been explained, but projecting to a new, unanalyzed region can be problematic).

Using multiple regression for the capsaicin analogs, one can derive the following equation which relates hydrophobicity and molar refractivity to biological activity.

Log EC₅₀ = 0.762 - (0.819) + (0.011)MR

s = 0.313, r² = 0.888

To judge the importance of a regression term, three items need to be considered.

1. Statistical significance of the regression coefficient.
   
   
2. The magnitude of the typical effect b_ix_i (in this case, 0.011 25.36).
   
   
3. Any cross-correlation with other terms.
   
   

As more terms are added to multiple linear regression, r² always gets larger. We recompute the previous calculations (r² = 0.89) carrying three significant figures so that rounding does not lead to confusion.

These results of this analysis indicate that, within this series, steric bulk is not an important factor in activity. The influence of the hydrophobicity constant confirms the presence of a hydrophobic binding site. Given the limited number of substituents in this analysis, it is unlikely that more can be learned from further analysis.

This section has developed the fundamental mathematics of QSAR studies. Several authors have published reviews of QSAR and have discussed various aspects of the methods^3-8. Each of the examples to follow uses these techniques to derive information about the chemical factors which are important for activity.

Approaches to Developing a QSAR

Drugs exert their biological effects by participating in a series of events which include transport, binding with the receptor and metabolism to an inactive species. Since the interaction mechanisms between the molecule and the putative receptor are unknown in most cases (i.e., no bound crystal structures), one is reduced to making inferences from properties which can easily be obtained (molecular properties and descriptors) to explain these interactions for known molecules. Once the relationship is defined, it can be used to aid in the prediction of new or unknown molecules.

The first approach to developing quantitative relationships which described activity as a function of chemical structure relied on the principles of thermodynamics. The free-energy terms E, H and S were represented by a series of parameters which could be derived for a given molecule.

Electronic effects such as electron donating and withdrawing tendencies, partial atomic charges and electrostatic field densities were defined by Hammett sigma () values, resonance parameters (R values), inductive parameters (F values) and Taft substituent values (*, *, E_s). Steric effects such as molecular volume and surface area were represented by values calculated for Molar Refractivity and the Taft steric parameter. Enthalpic effects were calculated using partition coefficients (LogP) or the hydrophobic parameter, , which was derived from the partition coefficient. In addition, an assortment of structural indices were used to describe the presence of specific functional groups at positions within the molecule. The linear equation which described the relationship between activity and this parameter set was the Hansch equation

log 1/[C] = A(logP) - B(logP)² + C(E_s) + D() + E + ...

Multiple linear regression analysis was used to derive the values of the coefficients. In general, Hansch type studies were performed on compounds which contained a common template (usually a rigid one such as an aromatic ring) with structural variation limited to functional group changes at specific sites.

Hansch utilized this approach in his analysis of 256 4,6-diamino-1,2-dihydro-2,3-dimethyl-1-(X-phenyl)-s-triazines which were active against tumor dihydrofolate reductase⁹. It was demonstrated that for 244 of the compounds, activity could be correlated to the presence of hydrophobic groups at the three and four positions of the N-phenyl ring. The parameters used to derive this correlation were the hydrophobic constant () and molar refractivity constant (MR) for meta and para substituents on the N-phenyl ring and six indicator variables I₁-I₆ which were used to indicate the presence (a value of 1) or absence (a value of 0) of specific structural features. The equation which was formulated from these data using the method of least squares is shown below.

FIGURE 5

Analysis of the Baker Triazines

log 1/[C] = 0.680(₃) - 0.118(₃)² + 0.230(MR₄) - 0.024(MR₄)² + 0.238(I₁) - 2.530(I₂) - 1.991(I₃) + 0.877(I₄) + 0.686(I₅) + 0.704(I₆) + 6.489

n = 244, r = 0.923, s = 0.377

The optimal values for MR₄ (4.7) and ₃ (2.9) were obtained from the partial derivatives of the equation. Note that the number of compounds in the data set was reduced to 244. Hansch and Silipo reported improvements in the value for r and s by removing 12 compounds which were incorrectly predicted by a factor of 10 or more.

While there are limits to the Hansch approach, it permitted complex biological systems to be modeled successfully using simple parameters. The approach has been used successfully to predict substituent effects in a wide number of biological assays. The main problem with the approach was the large number of compounds which were required to adequately explore all structural combinations. Further, the analysis methods did not lend themselves to the consideration of conformational effects. Several authors have published articles which provide additional background on the Hansch approach^10-11.

Alternative approaches to compound design have been suggested which avoid the combinatorial problem found in Hansch type analyses. Free and Wilson used a series of substituent constants which related biological activity to the presence of a specific functional group at a specific location on the parent molecule¹². The relationship between biological activity and the presence or absence of a substituent was then expressed by the following equation:

Activity = A + _ijG_ijX_ij

where A was defined as the average biological activity for the series, G_ij the contribution to activity of a functional group i in the jth position and X_ij the presence (1.0) or absence (0.0) of the functional group i in the jth position.

The procedure used the equation above to build a matrix for the series and represented this matrix as a series of equations. Substituent constants then were derived for every functional group at every position. Statistical tests were used to test the importance of the constants. If the models were shown to be valid, the model was used to predict activity values for compounds which had not been prepared. In general, while a large number of compounds are required to explore the effects of multiple substitution patterns, the Free-Wilson approach substantially reduces the number of analogs required. However, the method demands that the effects of substituents are additive.

In 1972, John Topliss published a paper which detailed methodology to automate the Hansch approach². The method assumed that the lead compound of interest contained at least one phenyl ring which could serve as the template for functional group modifications. The first modification to the template was preparation of the para-chloro derivative to examine lipophilicity. Additional substitution patterns were then made sequentially in an attempt to explore and optimize the relationship between activity and the hydrophobic and electronic character of the molecule. While the Topliss approach is easy to follow, it has several drawbacks. The primary problems are that the procedure is not applicable to all types of studies and that there is a high degree of risk associated with its use (it essentially ignores the possibility of interactions between substituents as it changes one substituent at a time).

The use of classical QSAR was expanded during the 1960's as a means of correlating observed activity to chemical properties. However, there are many areas where these techniques could not be used or where they failed to provide useful correlations. These included situations in which activity was found to be determined by 3-dimensional geometry, where poor training sets of compounds were used or the set of compounds were too small or insufficiently diverse and cases where biological activity could not be well quantified. Many of these problems were addressed by extensions to the Hansch method and the development of alternative approaches to QSAR.

There are cases where biological activity values cannot be determined accurately for a variety of reasons, e.g. lack of sensitivity of a particular test system. Alternative statistical techniques can be used in these cases; the problem is simplified to a classification scheme in which compounds are labeled as active, partially active, inactive, etc. The resulting data set is then searched for patterns which predict these categories. The methods which have been used for this type of analysis include SIMCA (Soft, Independent Modeling of Class Analogy)¹³, ADAPT (Automated Data Analysis by Pattern recognition Techniques)¹⁴, CASE (Computer Automated Structure Evaluation)¹⁵ and CSA (Cluster Significance Analysis)¹⁶.

Pattern recognition methods¹⁷ attempt to define the set of parameter values which will result in clustering compounds of similar activity into regions of n-dimensional space. The methods used to accomplish this goal can be parametric or nonparametric. Parametric methods search the n-dimensional space for clusters of compounds based on their calculated properties. These methods do not use derived values (e.g., mean vectors and covariance matrices), but instead use the original data to find clustering definitions and apply iterative procedures to find the linear set of parameters which best define the classification scheme.

Where the methods described above develop discriminant functions, SIMCA methods use Principal Component Analysis (PCA) to describe the data set. The objective of PCA is to create a reduced number variables which describe biological activity or chemical properties into a relatively few independent ones. This is accomplished through an analysis of the correlation matrix of biological or chemical properties.

Principle component analysis can be used to create derived variables for each class (e.g., active and inactive) separately by decomposing the correlation matrix; this method is useful to point out redundancies or interrelationships among the variables. PCA seeks to find simplified relationships in data by transforming the original parameters into a new set of uncorrelated variables which are termed principal components. The symmetric correlation matrix is decomposed by an eigenvalue decomposition. The largest eigenvalue and its eigenvector are used to form a linear combination of the original variables with maximum variance. Successively smaller eigenvalues and vectors produce linear combinations of the original variables with diminishing variance. Successive eigenvectors are independent of one another. The simplification is derived by disregarding eigenvectors associated with small eigenvalues. In summary, the procedure finds the set of orthogonal axes for the data which decompose variance in the data.

Another approach to examining the effects of chemical structure on activity was developed by the Jurs' group. Rather than rely on multivariant statistics to highlight these relationships, Jurs used the combination of cluster analysis and pattern recognition techniques as a tool to develop these correlations. The ADAPT program generated a data set of molecular descriptors (topological, geometrical and physicochemical) derived from three dimensional model building, projected these data points onto an n-dimensional surface and analyzed them using pattern recognition methods. The goal of this analysis was to discriminate between active and inactive compounds in a series.

Jurs has reported several applications of the methodology contained in ADAPT. In one study of chemical carcinogens¹⁸, a linear discriminant function was derived from a set of 28 calculated structure features including fragment descriptors, substructure descriptors, environment descriptors, molecular connectivity descriptors and geometric descriptors. Two hundred and nine compounds from twelve structural classes (130 carcinogens, 79 noncarcinogens) were selected for this study. The program was used to identify a training set of 192 compounds which was used to find the best set of descriptors and analyze the entire data set. A predictive success of 90% for carcinogenic compounds and 78% for noncarcinogenic compounds was obtained in randomized testing.

The CASE program extended the techniques in ADAPT by using topological methods to define substructural fragments which were essential for activity. CASE was able to differentiate between positional isomers. Both CASE and ADAPT are limited to analyzing structurally similar data sets.

The analysis methods described to this point have not explicitly incorporated the contribution of three dimensional shape in the analysis of the activity of a molecule. While the use of chemical graph indexes¹⁸, intermolecular binding distances¹⁹, molecular surface areas²⁰ and electrostatic potentials²¹ contain some information about the 3-D shape of molecules, the Hopfinger²² and Marshall²³ groups were the first to exhaustively analyze these effects.

In 1979, Marshall extended the 2-D approach to QSAR by explicitly considering the conformational flexibility of a series as reflected by their 3-D shape²³. The first step of the Active Analog Approach was to exhaustively search the conformations of a compound which was highly active in a particular biological assay. The result of the search was a map of interatomic distances which was used to filter the conformational searches of subsequent molecules in the series. The implicit assumption of the method was that all compounds which display similar activity profiles were able to adopt similar conformations. Once the "active conformation" was determined, molecular volumes for each molecule were calculated and superimposed. Regression analysis of the volumes was used to establish a relationship to biological activity. Marshall and co-workers commercialized the Active Analog Approach and a suite of other drug design techniques in the SYBYL molecular modeling program.

Hopfinger and co-workers also used 3-D shape in QSAR. In molecular shape analysis²⁴ of the Baker Triazines, the common space shared by all molecules of a series and the differences in their potential energy fields were computed. When these calculations were combined with a set of rules for overlapping the series, comparative indicies of the shape of different molecules were obtained. Inclusion of these shape descriptors in standard Hansch analysis schemes lead to improved descriptions relating computed parameters to biological activity such that no compounds in the original data set had to be eliminated from the calculations. The techniques developed by Hopfinger and co-workers were made available in the CAMSEQ, CAMSEQ-II, CHEMLAB and CAMSEQ-M computer programs.

In 1988, Richard Cramer proposed that biological activity could be analyzed by relating the shape-dependent steric and electrostatic fields for molecules to their biological activity²⁵. Additionally, rather than limiting the analysis to fitting data to a regression line, CoMFA (Comparative Molecular Field Analysis) utilized new methods of data analysis, PLS (Partial Least Squares) and cross-validation, to develop models for activity predictions.

The approach used in the CoMFA procedure requires that the scientist define alignment rules for the series which overlap the putative pharmacophore for each molecule; the active conformation and alignment rule must be specified. Once aligned, each molecule is fixed into a three-dimensional grid by the program and the electrostatic and steric components of the molecular mechanics force field, arising from interaction with a probe atom (e.g., an SP³ C atom), are calculated at intersecting lattice points within the 3-D grid. The equations which result from this exercise have the form

Act₁ = Const₁ + a₁(steric_xyz) + b₁(steric_xyz) + ... + a'₁(estatic_xyz) + b'₁(estatic_xyz) + ...

Act₂ = Const₂ + a₂(steric_xyz) + b₂(steric_xyz) + ... + a'₂(estatic_xyz) + b'₂(estatic_xyz) + ...

Act_n = Const_n + a_n(steric_xyz) + b_n(steric_xyz) + ... + a'_n(estatic_xyz) + b'_n(estatic_xyz) + ...

Traditional regression methods require that the number of parameters must be considerably smaller than the number of compounds in the data set (or the number of degrees of freedom in the data). The data tables which result from CoMFA analysis have far more parameters than compounds. PLS, which removes this limitation, is used to derive the coefficients for all of the steric and electrostatic terms. PLS essentially relies upon the fact that the correlations among near parts of a molecule are similar so that the real dimensionality is smaller that the number of grid points. Since these coefficients are position dependant, substituent patterns for the series are elucidated which define regions of steric bulk and electrostatic charge associated with increased or decreased activity. The size of the model (the number of components²⁷ needed for the best model) and the validity of the model as a predictive tool are assessed using cross-validation.

As opposed to traditional regression methods, cross-validation evaluates the validity of a model by how well it predicts data rather than how well it fits data. The analysis uses a "leave-one-out" scheme; a model is built with N-1 compounds and the N^th compound is predicted. Each compound is left out of the model derivation and predicted in turn. An indication of the performance of the model is obtained from the cross-validated (or predictive) r² which is defined as

r² (cross-validated) = (SD - Press)/SD

SD is the Sum-of-Squares deviation for each activity from the mean. Press (or Predictive Sum of Squares) is the sum of the squared differences between the actual and that predicted when the compound is omitted from the fitting process.

As we have discussed, values for conventional r² range from 0 to 1. Values for the cross-validated r² are reported by the method to range from -1 to 1. Negative values indicate that biological activity values are estimated by the mean of the activity values better than they are by the model (i.e, the predictions derived from the model are worse than no model). Once a model is developed which has the highest cross-validated r², this model is used to derive the conventional QSAR equation and conventional r² and s values. The results of the final model are then visualized as contour maps of the coefficients.

The first CoMFA study reported analyzed the binding affinities of 21 steroid structures to human corticosteroid-binding globulins and testosterone-binding globulins. This class of compounds is rigid and was selected to eliminate conformationally dependant effects from the study. The models for each steroid were built from coordinates from the Cambridge Crystallographic Database which were minimized using the Tripos force field. Side chain positioning was accomplished using systematic conformational searching. The Field Fit algorithm was used to align each structure within the fixed lattice (the 3-D grid used to calculate the CoMFA field effects). The fit of the regression line for the predicted versus actual binding values for the corticosteroids showed a cross-validated r² of 0.65 (conventional r² = 0.897, s = 0.397). For the testosterone-binding steroids, the cross-validated r² was 0.555 (conventional r² = 0.873, s = 0.453).

As noted, CoMFA starts with defined pharmacophore and overlap rules and derives a 3-D model which can be used to predict activity for new chemical entities. The Apex and Catalyst (Accelrys, Formerly Molecular Simulations Incorporated) programs are used to identify pharmacophores from databases of chemical structures and biological activity. These models are then used to predict activities for novel compounds.

Apex-3D is an automated pharmacophore identification system which can identify possible pharmacophores from a set of biologically active molecules using statistical techniques and 3-D pattern matching algorithms. The program classifies molecular structures using three methods: the agreement inductive method identifies common structural patterns in compounds having similar activity; the difference inductive method identifies structural patterns which differentiate active and inactive compounds and the concomitant variations inductive method highlights variations in structural features that explain changes in biological activity for sets of compounds.

The methods defined above follow logic similar to that used by a practicing medicinal chemist: What pharmacophoric patterns are present in the active molecule which are not present in the inactive ones? Pharmacophores are defined by different chemical centers (atom centered functional groups) and the distances between these centers. These descriptor centers can include such things as aromatic ring centers, electron donor ability, hydrogen bonding sites, lipophilic regions, and partial atomic charge. The information for each molecule is stored in a knowledge base in the form of rules which can be used to predict the activity of novel structures.

Apex-3D contains an expert system which automatically selects the best conformation and alignment for structures based on identified pharmacophores. When quantitative biological data is available, a 3D-QSAR model can be developed for any of the possible identified pharmacophores. Depending on the type of biological activity available, it is possible to identify pharmacophores for different binding orientations, receptor subtypes, or agonist versus antagonist activity.

To use this approach, the scientist is required to assign the training set of compounds to one or more activity classes Specific descriptor centers also can be defined if desired. The automated pharmacophore identification portion of the program then builds the knowledge base using the following steps:

• identify all possible binding interaction centers for each compound in the data set;
  
  
• generate topological (2D) or topographical (3D) distance matrices based on the set of descriptors;
  
  
• identify possible pharmacophores from all pairs of molecules using clique selection algorithms;
  
  
• classify these pharmacophores based upon their occurance in compounds in each activity class using Bayesian statistics and their nonchance occurance;
  
  
• set thresholds for probability and reliability statistics associated with a pharmacophore so that all training set molecules are properly classified by the pharmacophore rules;
  
  
• align compounds containing high probability pharmacophores on the pharmacophore.
  
  

Once the knowledge base has been constructed, the scientist can use it to predict biological activity for compounds not included in the training set.

The pharmacophores defined above can be used to build 3D QSAR models by correlating indexes calculated for biophore sites, secondary sites or whole molecule properties. 3-D QSAR models in Apex are generated and screened using a modified scheme of multiple linear regression analysis with variable selection. Special randomization checks are made to estimate the chance of fortuitous correlation. The steps involved are:

• Interactively select the pharmacophore(s) to use in the analysis;
  
  
• Interactively choose parameters to include in the pool of possible parameters to be selected by the program and
  
  
• Calculate the best 3D QSAR models for each pharmacophore using stepwise multiple regression and analysis of statistics to assess the validity and predictive power of the model.
  
  

At this point the knowledge base with 3D QSAR models can be used to calculate activities of novel compounds. Given a set of conformations for a novel structure, an activity range will be calculated based on each conformation which contains one of the pharmacophores and fits the 3D QSAR model.

Few applications of Apex have appeared to date in the literature. One example applied to nonpeptide Angiotensin II antagonists will be discussed for its heuristic value. Several structurally diverse compounds in this activity class have been reported in the literature^31-44. Pharmacophore models have been postulated with some disagreement about whether all of the highly active molecules are binding at the same site³⁰. Automated pharmacophore identification can be used to analyze these compounds and to assess the probability that they could be acting at the same site. A set of 55 compounds with specific binding activity (IC₅₀) values ranging over 6 orders of magnitude was used. Multiple conformations were included using 3D structures whose geometries were optimized by AMPAC. Four activity classes were defined; the most active class (< 100 nM) included 27 compounds. Apex-3D was able to generate rules which properly classified all compounds in the most active class without false negatives or positives. The fact that several biophores were required is consistent with the exsistence of multiple binding sites.

Quantitative 3D QSAR models also were developed for some of these pharmacophores. One model which contained 48 compounds had the following statistical parameters: predicted R² = 0.83, predicted RMSE = 0.86. Predicted activities for compounds excluded from the training set were within statistical boundaries.

Like Apex, Catalyst generates structure-activity hypotheses from a set of molecules of various activities. Once molecular connectivity and activity values are specified for all molecules, Catalyst derives a hypotheses which consist of a set of generalized chemical functions (regions of hydrophobic surface, hydrogen bond vectors, charge centers, or other user-defined features) at specified relative positions. Up to ten hypotheses are produced and ranked by estimated statistical significance. The hypotheses can be examined graphically, fit to new molecules, or fed directly to flexible 3D database search.

In the first step of the process, a set of representative conformers is found that covers the low-energy conformational space of each molecule. Representative structures are chosen to maximize the sampling of conformational space. The second step locates a list of candidate hypotheses that are common among active and rare among inactive compounds (the presence of all identified features is not required for inclusion in the active class since it is unlikely that all active members of a training set possess all the important binding features). The cost of a hypothesis is defined as the number of bits needed to describe the hypothesis as well as the errors in activities as estimated by the hypothesis. The theory of minimum complexity estimation indicates that a predictive hypothesis will minimize this cost. Minimization is carried out over the space of hypotheses covered by the candidates identified above. Statistical significance of the results (low probability of having found a chance correlation) on a variety of medicinal training sets is verified by a non-parametric randomization test.

The Genetic Function Approximation (GFA) algorithm is a novel technique for constructing QSAR models⁴⁵. It was specifically developed for use with data sets containing many more variables than samples, or data sets which contain nonlinear relationships between the variables and the activity.

GFA begins with a population of randomly-constructed QSAR models; these models are rated using an error measure which estimates each model's relative predictiveness. The population is evolved by repeatedly selecting two better-rated models to serve as parents, then creating a new child model by using terms from each of the parent models. The worst-rated model in the population is replaced by this new model. As evolution proceeds, the population becomes enriched with higher and higher quality models.

The different models are, in effect, multiple fits to the data. Scientists can then use their scientific knowledge and intuition to select among the final models. By studying the similarities and differences among the models, the scientist may be provoked to consider alternative mechanisms which explain the data or may plan new experiments to decide between different possible models. Since the models usually use only a subset of the total number of variables in the data set, GFA behaves as a member of the class of variable reduction algorithms and is most useful if the critical information in the data set is concentrated in a few variables.

Experimental results against published data sets demonstrate that GFA discovers models which are comparable to, and in some cases superior to, models discovered using standard techniques such as stepwise regression, linear regression, or partial least-squares regression. As the number of computational and instrumentational sources of experimental data increase, the ability of GFA to perform variable reduction, to discover nonlinear relationships, and to present the user with multiple models representing multiple interpretations of the data set, may become increasingly useful in data analysis.

There are many additional properties beyond binding which contribute to the biological activity of a molecule including transportation, distribution, metabolism and elimination. Since the techniques described above address binding and typically do not account for any of these other processes, they often fail. However, one should note that their utility is one of efficiency and probability. They are generally more successful and efficient than are the ad hoc/intuitive methods. Because no single method works for all cases, many groups are examining alternative approaches to developing SAR equations including electronic based descriptors and topological indicies.

Peter Goodford has reported an energy based grid approach to compound design. The GRID program⁴⁶ is a computational procedure for detecting energetically favorable binding sites on molecules of known structure. It has been used to study arrays of molecules in membranes and crystals or proteins. The energies are calculated as the electrostatic, hydrogen-bond and Lennard-Jones interactions of a specific probe group with the target structure. GRID has also been used to distinguish between selective binding sites for different probes.

In general, topological approaches start with a very different graph theoretic representation of molecular structure in which atoms are represented by vertexes and bonds by edges. Numerical indexes for the structure are then defined which abstract information including atom descriptors (atom types, atomic weight, atomic number, ratio of valance electrons to core electrons, etc.) and sub-group descriptors. In addition, a series of indexes (the indexes) are developed which describe the molecule as a set of fragments of varying size and complexity. Regression analysis of equations which have chi indexes as parameters has been used to correlate chemical structure to physicochemical behavior in applications such as chromatographic retention times, molar refractivity, ionization potential and heats of atomization. While the application of topological indicies has been widely reported in QSAR^47-49, the utility of these methods generally has been limited to predicting structure-property relationships for polymers and hydrocarbons.

SUMMARY

Developing a quantitative structure activity relationship is difficult. Molecules are typically flexible and it is possible to compute many possibly useful properties that might relate to activity. Early in a research program there are typically few compounds to model. Thus we have a few compounds in a very high dimensional descriptor space. Which are the important variables and how do we optimize them? It is clear that many training compounds need to span through the space and model fitting techniques need to address not only deriving a fit, but the predictive quality of the fit. While these methods have not discovered a new compound, they have aided scientists in examining the volumes of data generated in a research program. As the methods evolve, they will find broader application in areas such as combinatorial chemistry.

REFERENCES

1. Christopher S.J. Walpole, Roger Wrigglerworth, Stuart Bevan, Elizabeth A. Campbell, Andy Dray, Iain F. James, Kay J. Masdin, Martin N. Perkins and Janet Winter, J. Med. Chem., 36, 2381 (1993).
   
   
2. John G. Topliss, Utilization of Operational Schemes for Analog Synthesis in Drug Design, J. Med. Chem., 15, 1006 (1972).
   
   
3. Yvonne C. Martin, A Practitioner's Perspective of the Role of Quantitative Structure Activity Analysis in Medicinal Chemistry, J. Med. Chem., 24, 229 (1981).
   
   
4. John G. Topliss, Quantitative Structure-Activity Relationships of Drugs, Academic Press, New York, (1983).
   
   
5. Franke, R., Theoretical Drug Design Methods, Elsevier, Amsterdam, (1984).
   
   
6. Seydel, J.K., QSAR and Strategies in the Design of Bioactive Compounds, VCH, Weinheim, (1985).
   
   
7. Yvonne C. Martin, Accounts of Chem. Res., 19, 392 (1986).
   
   
8. Corwin Hansch, Accounts of Chem. Res., 2, 232 (1969).
   
   
9. Corwin Hansch and Carlo Silipo, J. Amer. Chem. Soc., 97, 6849 (1975).
   
   
10. Robert F. Gould (ed.), Biological Correlations -- The Hansch Approach, Advances in Chemistry Series, No. 114, American Chemical Society, Washington, D.C., 1972.
    
    
11. Yvonne C. Martin, Quantitative Drug Design, Marcel Dekker, New York, (1978).
    
    
12. Spencer M. Free and James W. Wilson, A Mathematical Contribution to Structure-Activity Studies, J. Med. Chem., 7, 395 (1964).
    
    
13. Svante Wold, Pattern Recognition by Means of Disjoint Principal Components Models, Pattern Recognition, 8, 127 (1976).
    
    
14. A.J. Stuper, W.E. Brugger and P.C. Jurs, Chemometrics: Theory and Application, B.R. Kowalski (Ed.), American Chemical Society, Washington, D.C., 1977.
    
    
15. G. Klopman and M.I. Dimayuga, J. Comput.-Aided Mol. Design, 4, 117 (1990).
    
    
16. J. McFarland, J. Med. Chem.,
    
    
17. Peter C. Jurs, Chemometrics and Multivariate Analysis in Analytical Chemistry, in Reviews in Computational Chemistry, Volume 1, K.B. Lipkowitz and B.B. Boyd (Eds.), VCH Publishers, Inc., New York, 1990.
    
    
18. Peter C. Jurs, J. T. Chou and M. Yuan, J. Med. Chem., 22, 476 (1979).
    
    
19. Lemont B. Kier, Molecular Orbital Theory in Drug Research, Academic Press, New York, NY, 1971.
    
    
20. Gordon M. Crippen, J. Med. Chem., 22, 988 (1979).
    
    
21. M. Mabilia, R.A. Pearlstein and A.J. Hopfinger, Computer Graphics in Molecular Shape Analysis, in Molecular Graphics and Drug Design, A.S.V. Burgen, G.C.K. Roberts and M.S. Tute (Eds.), Elsevier Science Publishers, Amsterdam, 1986.
    
    
22. H. Weinstein, in Chemical Applications of Molecular Electrostatic Potentials, Peter Politzer and Donald G. Truhlar (Eds.), Plenum Press, New York, NY, 1981.
    
    
23. R. Potenzone, E. Cavicchi, H.J.R. Weintraub and A.J. Hopfinger, Comput. Chem., 1, 187 (1977).
    
    
24. G.R. Marshall, C.D. Barry, H.E. Bosshard, R.A.. Dammkoehler and D.A. Dunn, The Conformational Parameter in Drug Design: The Active Analog Approach, in Computer Assisted Drug Design, ACS Symposia, 112, E.C. Olson and R.E. Christofferson (Eds.), American Chemical Society, Washington D.C., 1979
    
    
25. Anton J. Hopfinger, J. Amer. Chem. Soc., 102, 7196 (1980).
    
    
26. Richard D. Cramer III, David E. Patterson and Jeffrey D. Bunce, J. Amer. Chem. Soc., 110, 5959 (1988).
    
    
27. Components are defined as the linear combination of all parameters (the independent data) which are used to fit the shape of the parameters to the shape of the activity data (the dependent data).
    
    
28. V. E. Golender and A. B. Rozenblit, Logical Structural Approach to Computer Assisted Drug Design, in Drug Design, Vol IX, Academic Press (1980).
    
    
29. V. E. Golender and A. B. Rozenblit, Logical and Combinatorial Algorithms for Drug Design, Research Studies Press, UK (1983).
    
    
30. V. E. Golender and E. R. Vorpagel, Computer Assisted Pharmacophore Identification, in 3D QSAR in Drug Design: Theory, Methods and Applications, ESCOM Science Publishers, Netherlands (1993).
    
    
31. A. P. Thomas, C. P. Allott, et. al., "New Nonpeptide Angiotensin II Receptor Antagonists. 1. Synthesis, Biological Properties, and Structure-Activity Relationships of 2-Alkyl Benzimidazole Derivatives", J. Med. Chem., 35, 877-885 (1992).
    
    
32. R. H. Bradbury, C. P. Allott, et. al., "New Nonpeptide Angiotensin II Receptor Antagonists. 2. Synthesis, Biological Properties, and Structure-Activity Relationships of 2-Alkyl-4-(biphenylylmethoxy)quinoline Derivatives", J. Med. Chem., 35, 4027-4038 (1992).
    
    
33. R. H. Bradbury, C. P. Allott, et. al., "New Nonpeptide Angiotensin II Receptor Antagonists. 3. Synthesis, Biological Properties, and Structure-Activity Relationships of 2-Alkyl-4-(biphenylylmethoxy)pyridine Derivatives", J. Med. Chem., 36, 1245-1254 (1993).
    
    
34. P. R. Bovy, D. B. Reitz, et. al., "Nonpeptide Angiotensin II Antagonists: N-Phenyl-1H-pyrrole Derivatives Are Angiotensin II Receptor Antagonists", J. Med. Chem., 36, 101-110 (1993).
    
    
35. N. B. Mantlo, P. K. Chakravarty, et. al., "Potent, Orally Active Imidazo[4,5-b] pyridine-Based Angiotensin II Receptor Antagonists", J. Med. Chem., 34, 2919-2922 (1991).
    
    
36. W. T. Ashton, C. L. Cantone, et. al., "Nonpeptide Angiotensin II Antagonists Derived from 4H-1,2,4-Triazoles and 3H-Imidazo[1,2-b][1,2,4]triazoles", J. Med. Chem., 36, 591-609 (1993).
    
    
37. K. S. Atwal, S. Z. Ahmed, et. al., "Dihydro pyrimidine Angiotensin II Receptor Antagonists", J. Med. Chem., 35, 4751-4763 (1992).
    
    
38. J. V. Duncia, A. T. Chiu, et. al., "The Discovery of Potent Nonpeptide Angiotensin II Receptor Antagonists: A New Class of Potent Antihypertensives", J. Med. Chem., 33, 1312-1329 (1990).
    
    
39. B. De, M. Winn, et. al., "Discovery of a Novel Class of Orally Active, Non-Peptide Angiotensin II Antagonists", J. Med. Chem., 35, 3714-3717 (1992).
    
    
40. K. Kybo, Y. Kohara, et. al., "Nonpeptide Angiotensin II Receptor Antagonists. Synthesis and Biological Activity of Potential Prodrugs of Benzimidazole-7-carboxylic Acids", J. Med. Chem., 36, 2343-2349 (1993).
    
    
41. I. Sircar, J. C. Hodges, et. al., "Nonpeptide Angiotensin II Receptor Antagonists. 2. Design, Synthesis, and Structure-Activity Relationships of 2-Alkyl-4-(1H-pyrrol-1-yl)-1H-imidazole Derivatives: Profile of 2-Propyl-1-[[2'-(1H-tetrazol-5-yl)-[1,1'-biphenyl]-4-yl]-methyl]-4-[2-(trifluoroacetyl)-1H-pyrrol-1-yl]-1H-imidazole-5-carboxylic Acid (CI-996)", J. Med. Chem., 36, 2253-2265 (1993).
    
    
42. J. Weinstock, R. M. Keenan, et. al., "1-(Carboxybenzyl_imidazole-5-acrylic Acids: Potent and Selective Angiotensin II Receptor Antagonists", J. Med. Chem., 34, 1514-1517 (1991).
    
    
43. R. M. Keenan, J. Weinstock, et. al., "Imidazole-5-acrylic Acids: Potent Nonpeptide Angiotensin II Receptor Antagonists Designed Using a Novel Peptide Pharmacophore Model", J. Med. Chem., 35, 3858-3872 (1992).
    
    
44. R. M. Keenan, J. Weinstock, et. al., "Potent Nonpeptide Angiotensin II Receptor Antagonists. 2. 1-(Carboxybenzyl)imidazole-5-acrylic Acids", J. Med. Chem., 36, 1880-1892 (1993).
    
    
45. D. Rogers and A. J. Hopfinger, "Applications of Genetic Function Approximation (GFA) to Quantitative Structure-Activity Relationships (QSAR) and Quantitative Structure Property Relationships (QSPR)", J. Chem. Inf. Comp. Sci., 34, 854-866 (1994).
    
    
46. Peter J. Goodford, J. Med. Chem., 28, 849 (1985).
    
    
47. A.T. Balaban, Applications of Graph Theory in Chemistry, J. Chem. Inf. Comput. Sci., 25, 334 (1985).
    
    
48. Lowell H. Hall and Lemont B. Kier, The Molecular Connectivity Chi Indexes and Kappa Shape Indexes in Structure-Property Modeling, in Reviews in Computational Chemistry, Volume 2, K.B. Lipkowitz and D.B. Boyd (Eds.), VCH Publishers, Inc., New York, 1991.
    
    
49. I.B. Bersuker and A.S. Dimoglo, The Electron-Topological Approach to the QSAR Problem, inReviews in Computational Chemistry, Volume 2, K.B. Lipkowitz and D.B. Boyd (Eds.), VCH Publishers, Inc., New York, 1991.