Determinants of molecular reactivity as criteria for predicting toxicity: problems and approaches.

We discuss the physicochemical basis for mechanisms of action of toxic chemicals and theoretical methods that can be used to understand the relation to the structure of these chemicals. Molecular properties that determine the chemical reactivity of the compounds are proposed as parameters in the analysis of such structure-activity relationships and as criteria for predicting potential toxicity. The theoretical approaches include quantitative methods for structural superposition of molecules and for superposition of their reactivity characteristics. Applications to polychlorinated hydrocarbons are used to illustrate both rigid superposition methods, and methods that take advantage of structural flexibility. These approaches and their results are discussed and compared with methods that afford quantitative structural comparisons without direct superposition, with special emphasis on the need for efficient automated methods suitable for rapid scans of large structural data bases. Quantum mechanical methods for the calculation of molecular properties that can serve as reactivity criteria are presented and illustrated. Special attention is given to the electrostatic properties of the molecules such as the molecular electrostatic potential, the electric fields, and the polarizability terms calculated from perturbation expansions. The practical considerations related to the rapid calculation of these properties on relevant molecular surfaces (e.g., solvent- or reagent-accessible surfaces) are discussed and exemplified, stressing the special problems posed by the structural variety of toxic substances and the paucity of information on their mechanisms of action. The discussion leads to a rationale for the use of the combination of theoretical methods to reveal discriminant criteria for toxicity and to analyze the initial steps in the metabolic processes that could yield toxic products.


Introduction
The actions of xenobiotic agents in biological systems can be produced by a variety of specific molecular interactions and combinations of nonspecific processes. Interacting specifically with receptors and enzymes, the xenobiotics can mimic endogenous substances or interfere with the actions of hormones and neurotransmitters. Other effects can be caused by their penetration into cell membranes or organelles where the physicochemical properties of the agents can cause changes in the local environment (e.g., pH, viscosity) to trigger or modulate normal processes. All these actions reflect a close relation between the physicochemical properties encoded in the molecular structure of the agents and the responses they evoke in the biological systems. It becomes possible, therefore, to seek an understanding of the structure-activity relationship for the toxic actions of xenobiotics by trying to reveal the relation between molecular structure of chemicals and their effects on biological systems based on the molecular properties encoded in their structure. These properties determine the molecular reactivity of the agents, and are responsible for their recognition at biological acceptors and receptors, as well as for the triggering of molecular mechanisms that lead to tissue response.
The special difficulties presented by the analysis of molecular determinants for toxic effects, as opposed to specific physiological and pharmacological actions of drugs, stem from the enormous chemical variety of the agents and the fact that they were not designed to have selective biological effects. In addition, there are often very few experimental data available to evaluate the spectrum of their biological effects. The large number of agents that must be evaluated within short periods of time further reduces the feasibility of a complete evaluation of their effects on health and the environ-ment; in many cases only the chemical structure and the simplest physicochemical properties are available for the characterization of the agents. These difficulties emphasize the importance of theoretical criteria for the evaluation of potential toxic effects, and the special role that considerations based on chemical structure and molecular reactivity properties must play in these evaluations.
The basic assumption in the analysis of the relation between molecular structure and biological effect is that the interaction between exogenous compounds and biological targets is dependent on the same molecular properties that determine chemical interactions and reactions. Based on this assumption, a description of actions, interactions, and mechanisms should follow the formal patterns of chemistry, with the special parameters of biomolecular interactions, such as size and complexity of environment, determining a special class for these phenomena. As indicated by documented successes in the theoretical study of such biological interactions, a first step in the analysis of these complex phenomena is their separation into physically meaningful steps that are amenable to analysis by methods of theoretical and quantum chemistry (1). For example, in the complex phenomena of protein-protein interaction, entropic stabilization was shown to be the major factor in protein association due to the very large number of protein-water interactions replaced by the formation of a protein-protein interface; but the specificity of the interaction is entirely dependent of the complementarity of the interfacing protein surfaces (2). The onset of the stabilizing bulk effect of the solvent is thus determined by the discrete interactions between properly positioned residues, a mechanism dependent on the same factors that govern interactions of small molecules. The selectivity of the protein-protein interaction, like that of protein-ligand interactions, should therefore be describable by the formalisms of chemical reactivity; this is illustrated in many studies in which the interactions are treated as a superposition of discrete physicochemical components (3).
Another documented example is the study of enzymesubstrate interactions which are being investigated not only for the direct information they can provide on the particular systems, but also as heuristic models for mechanisms in other, structurally or functionally related biological systems. Thus, the theoretical studies of carboxypeptidase, (4)(5)(6)(7)(8)(9)(10)(11) and carbonic anhydrase (12,13), and their interactions with ligands such as substrates and inhibitors, are aimed at an understanding of the structural basis for specific mechanisms of action. These studies are complemented by structural and theoretical investigations of other enzymes such as chymotrypsin (13)(14)(15), and other serine proteases, as well as papain (16,17) and thermolysin (18,19), in which generalizations of the findings are sought to explain the actions of yet other classes of biological macromolecules (13,(20)(21)(22), including receptors (23). A variety of results from such studies yield inferences on the molecular mechanisms of recognition and selectivity of biological macromolecules. These mechanisms illustrate the molecular basis for biological structure-activity relationships expressed in terms of the reactivity characteristics of the external agent, represented by the ligands (e.g., substrates and inhibitors of enzymes), and of its biological target, represented by the recognition site (e.g., the active site of enzymes).
The ability to generalize conclusions obtained from the study of enzyme-substrate interactions to other biological mechanisms rests on the causal relationship between the properties of the molecules and their effects on biological systems (4). When the causal relationship between the physical properties of a specific molecule and the response of a biological system to that molecule has been sufficiently elucidated, that relationship may be generalized, and the capacity of other (perhaps structurally unrelated) molecules to elicit a similar response may be assessed, in a straightforward manner, from the physical properties of those molecules. However, for many toxic substances, detailed information about the mechanism for biological response is not available to perform this straightforward assessment. In these cases, physical properties of molecules may usefully be applied as reactivity characteristics, in a correlative approach to the prediction of activities of untested chemicals. The advantage of reliable physical properties as reactivity characteristics is that: (1) existing information about mechanisms of action (either of a general or of a specific nature) may be incorporated directly; (2) physical properties that relate to the potential of a molecule to interact with other biomolecules may be used with or without inferences for a specific mechanism; (3) when correlations are established, they can be used to provide insight into the mechanisms of action. The insight that stems from the solid physicochemical basis of the reactivity characteristics used as criteria provides some assurance that the correlations are not just fortuitous and can be used for the assessment of activity outside the initial set of molecules.
Since the molecular properties that are directly responsible for the molecular interactions leading to the biological effect are encoded in the entire molecular structure, the description of the chemical reactivity characteristics must also rest on a representation of the properties of the entire molecule. This consideration reveals one of the major drawbacks of the approaches in which structure-activity relationships (SAR) are studied on the basis of parameters representing contributions from certain structural components (functional groups, substituents). These parameters used in "quantitative SAR" (QSAR) often fail to provide a reliable relationship between molecular structure and biological activity for they seldom constitute an adequate representation of molecular properties (24). Useful predictive criteria for biological activity should therefore be constructed from molecular properties derived from the entire molecular structure because their role in determining various aspects of molecular activity can be understood on the basis of clear physicochemical principles. The possible use of such determinants of molecular reactivity as criteria for explaining and predicting toxicity must further be evaluated on the basis of mechanistic hypotheses linking molecular structure to molecular properties and to the types of molecular interactions determined by these properties.

Evaluation of Molecular Properties Determining Reactivity
Structural Considerations: Superpositions and Comparisons Methodological Framework. Some relationships between structure and function in molecules exhibiting biological activity can be revealed from correlations of activity with similarity in characteristics that describe structure -.g., size, shape, topology or conformation-or with similarity in geometrically defined elements of reactivity, e.g., dipole moments, electrostatic orientation vectors. Because physicochemical methods for the analysis of molecular structures, and especially X-ray crystallographic studies, have provided detailed structural information, including the atomic coordinates of many compounds [more than 35,000 compounds are now available through the Cambridge Crystallographic Data Base (25)], it becomes possible to obtain meaningful comparisons of both the structural and the geometrical parameters of the elements of reactivity of molecules that share biological properties. Results from such comparisons should yield answers to some fundamental questions underlying predictions from the structure-activity considerations: (1) given molecules of similar structure, how similar are their physical properties that might be responsible for recognition, specificity and reactivity? and (2) given that a particular physical property, or combination of more than one such property, can be established to be the determinant for recognition, specificity or reactivity, can two molecules, that may be structurally dissimilar, be conformationally altered to become functionally similar by adopting conformations in which their determinants for recognition are made to coincide?
Answering these questions requires the application of techniques based on structural superposition, either of rigid molecule or with conformational freedom, in addition to the computation of the molecular properties of the various structural forms. Some of the approaches for structural superposition are outlined and exemplified in this section, while some of the relevant molecular properties are presented below.
The methods used in structural superposition may be classified as operating on rigid or on flexible molecules, but in either case certain common aspects of the methods are retained. Thus, in both approaches it is necessary to first select a subset of the list of atoms of molecule 1 (it may be the entire list of atoms) that are to be directly equivalenced with a list of atoms in molecule 2. These two subsets form the "equivalenced atom list" for each molecule, and the superposition technique attempts to derive a rotation-translation procedure that will effect the optimization of their superposition as monitored by a function, most commonly the "rootmeans-square deviation in positions" of the two equivalenced atom lists. According to the criteria used in these comparisons, further classification is possible, into methods of strict structural comparison, including the Distance Matrix analysis (26,27) and the Partitioned Distance Matrix analysis (28), and methods comparing parameters related to properties depending directly on structure, including the Excluded Volume analysis (29) and the analysis of Surface Accessibility (30).
Illustration of the Approach: DDT Analogs. The polychlorinated hydrocarbons 1,1,1 ,2-tetrachloro-2,2bis(p-chlorophenyl)ethane (Cl-DDT) and 1,1,1-tribromo-2,2-bis(p-chlorophenyl)ethane (Br-DDT), are analogs of DDT that exhibit differences in their respective activity: Br-DDT is moderately toxic to mosquito larvae, while Cl-DDT is not (31). In addition, Br-DDT is dehydrohalogenated in vivo more rapidly than DDT, because it is an excellent substrate for DDT-ase (32). The structures of these molecules have been examined by X-ray crystallography to high resolution (33), and the comparisons presented below will be based on these structure (using only the coordinates of molecule A in the crystal of Br-DDT). RigidMolecule Superposition. The methods for determining the appropriate rotation-translation matrix have been described in detail by Cox (34), Rao and Rossman (35), Hendrikson and Ward (36), and Argos and Rossman (37). Attempts have been made to improve the sensitivity of this approach to enable it to discriminate structural subgroups that might superimpose more exactly even when the superposition of the entire molecules is not very good, e.g., when the compared mol-ecules are composed of very similar fragments that differ only in their spatial position. This limitation of the rigid molecule rotation is exemplified in the comparison of Br-DDT and Cl-DDT.
To effect the superposition, a matrix was generated by using he common carbon and hydrogen atoms of rings A and B, and the ethyl that links them, a total of 22 common atoms. The root-mean-square (RMS) fit of this equivalence set is 2.27 A (mean = 1.89 A) with the statistical analysis of the distribution of the fit among the set of equivalenced atoms indicating a single outlier, C14. Removal of this outlier and redetermination of the rotation-translation matrix yielded a final fit of 1.95 A (mean = 1.70 A). In this superposition the shortest distance between equivalenced atoms occurs between ring A, atom C3, with a deviation of 0.54 A, and ring B atom C8, also with a deviation of 0.54 A. The average deviation for all atoms in ring A is 1.28 A, and in ring B 1.38 A, with the deviation for the two carbon atoms in the ethyl group being 3.2 A for C13 and 5.8 A for C14. As stated above, the only conformational difference among those atoms that are common to Br-DDT and Cl-DDT result from the difference in the torsional angles, C8-C7-C13-C14 and C2-C1-C13-C14, with no significant alterations occurring in the conformations of the rings, themselves. This observation is readily verified by the superposition of rings A and B independently, yielding RMS and mean deviations of 0.04 A,0.04 A for ring A and 0.04 A, 0.03 A for ring B, respectively.
The conformational difference between the torsional angles thus yields an average deviation of 4.47 A in the position of the second ring, not included in the single ring superposition. This comparison clearly indicates the reasons for the inability of the rigid molecule superposition to reveal the high degree of structural analogy that exists between the components of the two molecules. This is a limitation that might obscure the reason for the difference in the biological properties of these DDT analogs by making it appear to be due to dissimilar structures. The alternative use of techniques of flexible superposition, discussed below, reveals the elements of structural similarity and thus focuses attention on the possible differences in the physicochemical properties of the molecules that are more likely responsible for the difference in their actions.
Molecular Superposition with Internal Group Rotation. In the flexible superposition of molecules (i.e., of rigid substructures joined by rotatable bonds), free rotation around the bonds linking rigid substructures is permitted during the superpositioning procedure. Comparison of two molecules in this manner does not determine if the molecules are superimposable in their original conformations, but rather if they could be made to be structurally superimposable by use of available conformational flexibility. This method of structural comparison can also identify the similarity present in nonbonded (through-space) orientations between noncontiguously linked substructures, as might be required for recognition of functional groups at biological recep-tors or enzyme active sites. Consequently these techniques permit the search for similar topological features that may be attainable in structures whose internal chemical and structural linkages may differ significantly.
The traditional approach to structural superposition of conformationally flexible molecules performs rotations about bonds connecting sequentially linked rigid substructures. In this manner, conformational space is sampled by the progressive rotation of each rotatable bond; following each rotational step, superposition is performed by generating a rotation-translation matrix, and the quality of the superposition is evaluated and compared with the global fit (or of a segment of the structure). Such a technique, the Burlesque approach (38), has been used to compare distances from atoms within the conformationally flexible molecule to a fixed external point. Frequently, this reference point bears some physical or mechanistic significance, or is derived from experimental observations (e.g., paramagnetic resonance shifts, etc.) for the nonaltered structure. The practical limitation of this approach lies in the need to sample an extremely fine grid to adequately approximate a continuous conformational search. The number of conformations to be generated in this manner is a function of both the size of the search interval and the number of bonds to be searched, and has been shown to grow rapidly: a three-angle search using a 10-degree sampling interval on each requires 46656 conformations to be individually generated and tested. This approach assumes that the conformational space will be adequately sampled by this search grid although it is possible that the region that satisfies the experimental criteria will be narrower than the search grid and therefore no solution will be found. Improvements to this approach use segmented, overlapping searches and have been described (39).
In another approach, proposed by Barino (40), rotation matrices are derived for each of the substructure fragments and combined with the necessary translations to superimpose two molecules. This method attempts to generate only those conformations of the flexible molecule that are superimposable onto the molecule of fixed conformation. This application of flexibility-linked rigid segments has been adapted to permit superposition with noncompletely equivalenced lists of atoms (i.e., groups differing in total number of atoms). Refinement of the respective rotation-translation matrices is based on the same root-mean-square deviation statistic described above for the rigid molecule superposition. This approach can be extended to include energy considerations in the weighting of conformations, along with the root-mean-square statistic for emphasis on realistic, energetically feasible structures (40).
An even more general approach is based on the use of distance geometry where the distance matrix construction readily permits incorporation of a variety of structural characteristics into the superposition algorithm (41). This approach does not require the direct superposition of one molecule onto another but rather attempts to generate a conformation that is in agreement with the distance constraints imposed by the interatomic distances within the first molecule to which the second molecule is being fit. The primary advantage of this technique is the ability to generate only those structures that fit the exact structural criteria within preset boundary conditions. It is thus possible to weight these characteristics rather than generate the large set of conformations that poorly span conformational space, and to select for analysis only those that fit specific criteria as produced by the Burlesque approach described above. The adaptation of this approach to incorporate functional descriptors (e.g., molecular properties such as the electrostatic potential and its derivatives) reaches beyond the analysis of structural or through-space orientations alone. An illustration of the conformational search based on such reactivity criteria is presented in a subsequent section.
In our general example, which compares Br-DDT with Cl-DDT, it is apparent that any of the three techniques described above should permit the determination of the rotations about C13-C1 and C13-C7 that are necessary to superimpose the equivalent atoms in the two analogs.
It is to be noted that the gridded rotational search as carried out in the Burlesque technique would require the use of a relatively fine search grid to achieve optimal structural superposition. All the methods described in this section would require the definition of rigid substructures, but unlike the distance geometry approach which does not require direct superposition, these methods are not well suited for integration into an automated algorithm for structural comparison. In the development of a computer-based system for correlation of structure/function/toxicological properties, it is especially interesting to be able to utilize the methods of pattern recognition and automated structural comparison in data base searching. Work in this area is ongoing in several other laboratories (42), in addition to ours.
Other techniques designed to provide a solution to the limitations inherent in superposition procedures are presented in the following section, where we discuss alternative methods for comparison of structure based on the distance matrix (i.e., pattern recognition procedures) that circumvent the need for structural superposition.
Structural Comparison without Direct Superposition. The example of Br-DDT and Cl-DDT indicates that the single statistic, root-mean-square deviation of superposition is inadequate to fully describe the quality of the fit of two compared molecules. Even with the added considerations and improvements of the superposition procedures outlined in the section above, a higher resolution form of comparison and analysis is needed, particularly one that may prove useful for automated structural search algorithms required for the search of large data bases to reveal compounds related in their biological or toxicological activities. Distance matrix analysis (27,43) provides such an alternative and has proven useful in both small molecule and protein structure analysis and comparison. In this representation, a given molecular structure is transformed into a square matrix of order N, where N is the total number of atoms. The atoms are first ordered into a list, and the resultant matrix element, ij, represents the distance between atoms i and j of the list. The distance matrix for a given molecule is symmetric about the diagonal as the distance between atoms i and j is equivalent to that between atoms j and i; the diagonal elements, i =j, are all 0. The most important aspect of the distance matrix for use in structural pattern matching is the property that it is invariant to translations and rotations of the whole molecule. This-property enables the direct comparison of two molecules, or substructures, without the need to translate and rotate one onto the other. This technique has proven powerful in its ability to reveal conformational perturbation, and ascribe it either to changes in substructure conformations or to changes in the through-space orientations between these substructures (29,44,45). This approach simultaneously and independently compares the configuration and conformation within the two molecules, to provide quantitative structural comparisons that are suitable for structure-activity correlation analysis (44,45).
For the comparison of Br-DDT and Cl-DDT, the individual distance matrices were generated as well as the difference distance matrix, as shown in Figure 1. This difference distance matrix was computed as an element-wise difference between the matrices for the two compounds, and is annotated in Figure 1 to indicate the regions concerned with the conformation of the individual rings, A and B, the ethyl group linking the rings, as well as the regions concerned with the respective orientations of these substructures. The distance matrix is symmetric, so we can use the upper half matrix to represent two other significant metrics, derived from the lower half matrix elements: the absolute and signed sums of deviations within the respective partitions. This metric was generated by summing the signed values within the partition, and also by summing the absolute values of the elements in each partition. From their definition, it is clear that when the magnitudes of the sums are comparable, then the conformational change is concerted (i.e., a majority of elements within the partition contract or expand); when the magnitudes are such that the signed sum is approximately equal to 0.0 while the absolute sum is significantly greater than 0.0, then the conformational difference is of a symmetrical origin (e.g., rotation about a bond of symmetrical group, etc.).
Results in Figure 1 show that the elements within either ring A or ring B deviate little from 0.0 (i.e., absolute/signed sum are 0.4/0.4 A for both ring A and ring B), whereas the partition concerned with the orientation of ring A vs. ring B yields an absolute sum of 15.4 and a signed sum of 1.4, indicating the symmetrical nature of the conformational difference that arises from the difference in rotation of the phenyl rings with respect to the ethyl. This rotational difference, which we 6    know involves the C13-C1 and C13-C7 bonds, produces approximately alternating lengthening and shortening of distance pairs within the matrices and gives rise to the observations made above. Thus, the partitioned distance matrix technique made it immediately apparent that the respective conformations of the substructures in Br-DDT and Cl-DDT remain constant although the through-space orientation of their rings A and B are altered. A more detailed description of this form of analysis and its application to other DDT analogs is presented elsewhere (45).
Comparison of Surface Accessibility. The surface area of a molecule that is susceptible to interactions with a reactant agent, and especially the steric accessibility of certain key atoms or groups to such an attack are additional elements in the set of molecular reactivity characteristics that could determine biological activity. Their use in the study of structure-activity relationships is well documented (30). The molecular reactivity criteria define reactive regions in the molecules and the likely orientation for certain types of molecular interactions (e.g., electrostatic, stacking, etc.). In order to add the consideration of spatial reactivity criteria we have implemented a computational procedure, SUR-VOL, to evaluate the total surface area surrounding a particular atom and to quantitate the relative accessibility of an atom or a molecular region. This rapid and efficient procedure is based on a Monte Carlo simulation of space filling within a box of enclosure. The box surrounds the molecule and is constructed de novo for each structure examined with this procedure. To construct the box, the list of atomic coordinates of the molecule is searched to determine the minimal and maximal values of the x, y, and z coordinates. The atoms belonging to these extrema are identified and the boundaries of the box are generated so that it will enclose the molecule within walls that have a clearance of two times (2 x ) the radius from the boundary atom. The volume of the box is then computed by filling it with a set of randomly generated points at a fixed density of 50 points per cubic angstrom (50 p/A3). This density was found to yield statistically significant reproducibility of molecular volumes in our examination of the effects of varying the density of points used to fill the box on the resulting calculated volume.
The list of coordinates of the atoms is then searched again to determine if the generated points lie within the van der Waals radius of any atom. Counting those points that fall within van der Waals radii and taking the ratio of this number to the total number of all points generated yields a ratio of the volume of the enclosed molecule to that of the entire box. In this manner, a molecular van der Waals volume can be calculated; it may be further extended to yield the solvent (or reagent)-accessible volume by the simple addition of the radius of the solvent (or reagent) probe to that of each atom.
The exposed surface area of each atom and the resulting molecular surface area are computed in an analogous manner to the volume computation described above. The van der Waals surface of each atom is simulated by a set of points randomly generated at uniform density at the fixed van der Waals radius. These points are then compared to the atom list to determine if they occur within the van der Waals radius of any other atom (i.e., if they are buried in an atom-atom interface), or if they are exposed. The ratio of points exposed to the total number of points yields a ratio of the exposed to theoretically available surface areas for each atom. This parameter can be measured rapidly for each atom or functional group in structurally dissimilar molecules to provide a profile that reflects the effects of structural differences on the molecular reactivity (i.e., the accessibility of the reactive regions in the molecule). The results obtainable with this approach are illustrated with a calculation of the reactive surface areas of the atoms in Br-DDT and Cl-DDT ( erties, molecular interactions, and biological activities has been documented and reviewed critically in a recent symposium dedicated to this topic (46). The symposium successfully presented the fundamental and applied research in electrostatic potentials and illuminated the strengths and limitations of this approach. In our contribution to the proceedings of the symposium we have illustrated the use of MEP in the prediction of biological activity (47), here we review some new approaches, and some of the problems in the use of MEP for the characterization of molecular reactivity properties related to the biological activities of molecules. The molecular properties that determine reactivity depend on the spatial distribution of molecular electronic and nuclear charges, which is an averaged, time independent distribution within the working approximations of quantum mechanics. A direct description of molecular properties relevant to biological recognition should therefore be obtainable as a description of the molecular reactive sites and interaction ability based on long-range intermolecular electrostatic forces. The electrostatic potential field generated by the distribution of charge in the molecule generates a "picture" of the forces that the molecule can apply on another charge or on the continuous electron density distribution of another molecule. Regions of space in which the molecule generates negative potentials (i.e., where the contributions from all the electrons in the molecules are being felt stronger than those from all its nuclei) should attract positively charged species or fragments. A positive region (where the nuclear contribution is stronger) would repel positive, and attract negatively charged fragments. A molecular species approaching the investigated molecule should therefore adopt an orientation that will bring its charged fragments into matching attractive interaction with the positive and negative potential pattern. On this basis, the calculation of molecular electrostatic potentials has been successful in first order predictions of the relative reactivity of functional groups in ionic-type reactions by the characterization of the shape and localization of regions representing attractive and repulsive zones of varying potency (46). As discussed in detail below, the spatial disposition and relative strength of these regions forms a pattern characteristic of a given molecule, in a certain geometric conformation, and can be analyzed to predict with some accuracy the reactivity of the whole molecule towards reagents of varying complexity (4 7,48).
Other reactivity criteria often used in chemistry and applicable to the study of molecular determinants for biological activities are properties such as net atomic charges, bond orders and bond polarities, bond polarizabilities and "frontier orbital" parameters. However, these reactivity criteria are often artificially constructed parameters with no obvious formal relationship to recognition and activity, whereas the electrostatic potential generated by a molecule is an experimentally observable quantity from high energy electron scattering (46)]. Thus, it is documented that quantum chemical calculations generate electrostatic potentials that agree with the "static potential" component measured in scattering experiments (49)(50)(51). We have therefore defined an interaction pharmacophore based on the MEP (52,53) and have demonstrated its analytical and predictive power in studies of drug actions in a variety of systems related to cholinergic muscarinic receptors (52)(53)(54), the actions of phencycidine (54,55), of hallucinogens (24,56), of histamine (57,58) and models for metalloenzymes (4-7).
For example, in a series of studies on the molecular determinants for recognition and activity at the receptors of the neurotransmitter serotonin (5-hydroxytryptamine; 5-HT) a major recognition element is the electrostatic potential (24,48). We showed that the MEP over the indole portion of 5-HT congeners, and of congeners of lysergic acid diethylamide (LSD), has a characteristic directionality that can be represented by an electrostatic orientation vector. By using 5-HT and LSD as templates it was possible to define geometrical requirements for recognition at the high affinity 5-HT receptor based on the spatial relationship between the electrostatic orientation vector in these molecules, and the site of attachment at the anionic part of the receptor (47,48). These recognition criteria transcend structural comparisons of the kind described above and constitute a new kind of predictive criteria for biological activity (47).
The large number of similar studies using MEP to elucidate determinants for drug actions (59-65) make a comprehensive review impossible in the present context. The common approach is, however, noteworthy: In all these examples some characteristics of the MEP are compared among active and inactive species to infer on the common properties that might be responsible for the activity. We have recently developed an approach in which the techniques of structural analysis and comparison, described in the previous section, are used to optimize and quantitate the comparison of molecular reactivity characteristics such as the MEP. This approach is illustrated here by a comparison of the spatial orientation of the major determinants for recognition of 5-HT analogs at serotonin receptors. As described above, these determinants were shown to include the cationic alkylamine side chain that is common to all congeners, and the MEP over their indole portion. Using the new technique we evaluated the constraints for positioning the analog 6-hydroxytryptamine (6-HT) in accordance with our hypothesis for 5-HT-like recognition at the receptor (48), by aligning the electrostatic vectors and allowing the side chains to interact with the same anionic site. To use the conformational search procedures described above, we selected a point along an N-H bond in the cationic side chain, at a distance of 2.85 A from the nitrogen, to approximate the position of the anionic site, and two other points chosen from the molecular reactivity properties: the minimum in the electrostatic potential near the hydroxyl oxygen, and a third point chosen along the electrostatic orientation vector at 5 A from the second one. The distances between the three points were: 1-2 = 8.3 A; 1-3 = 6.4 A; 2-3 = 5 A. The constraints used in the analysis permitted full rotation about the C(3)-C(beta) bond and the C(alpha)-C(beta) bond, at 10 degrees intervals, to search for a suitable conformational region which was then further searched at intervals of 1 degree. The maximal errors in the three distances mentioned were set to be 0.3 A to the points along the orientation vector and 0.1 A to the nitrogen. With these constraints, the structure of 6-HT was rotated about the equivalent bonds, i.e., the N-C(alpha) and the C(alpha)-C(beta) bonds. (The conformational search of the rotatable bonds is performed with the goal to superimpose the fixed parameters in the investigated molecule on those of the template.) The search was performed by sequentially stepping about each of the rotatable bonds, and subsequent evaluation of potential van der Waals interactions which act as constraints to prohibit configurations in which a preset tolerance is violated. Initial tolerances were set to permit 0.95 penetration (i.e., no pair of atoms may be closer than 0.95 times the sum of their van der Waals radii). A threedimensional translational minimization was then performed to satisfy the placement of the parameters within the distance constraints produced by the comparison with the template. The calculations were applied exhaustively until all conformations examined were either retained or eliminated. For the retained conformations, the comparison with the values found in 5-HT in the extended conformation was carried out to assess the success of the conformational search. The search yielded a conformation of 6-HT and associated anionic binding site representation which approximated the topological feature in 5-HT. A subsequent superposition of the 5-HT and 6-HT models using the two points along the electrostatic orientation vector and the anionic binding site yielded a root mean square deviation of 1.1 A overall. The superposition is shown in Figure 2.
This analysis shows that the criteria used to superimpose 5-HT and 6-HT based on the interaction of their cationic side chain with a common anionic binding site, and the requirement for parallel orientations of the electrostatic vectors, does not yield a full superposition of the indole portion of the two tryptamines. This indicates the importance of superposition criteria based on molecular properties that are likely to determine the interaction with biological targets, rather than criteria based on mere structural considerations as discussed above. Nevertheless, such superpositions of flexible molecules require support from a large number of comparisons and, whenever possible, reference to rigid analogs. Thus, the results from the search illustrated here were obtained for a side chain kept in the (tau 1 = 90, tau 2 = 180) conformation which is close to a minimum energy conformation, but is not necessarily the functional conformer for recognition (47). This conformation does not place the amine nitrogen of 5-HT at a position equivalent to that in LSD, the rigid molecule that is recognized with high affinity at 5-HT receptors (see above). A simultaneous fit to the rigid analog is therefore required to refine this functional superposition. Problems related to the environmental effects on molecular reactivity merit separate discussion (1), but it seems appropriate to mention here a specific remark related to the use of the interaction pharmacophore which is based on the calculated MEP. Thus, Hopfinger has suggested that an "obvious reservation to the interaction pharmacophore concept is the change in the electrostatic potential about a drug because of the "bound water on the drug" (3). Since it is difficult to assume that molecules could interact specifically "through" their specific hydration shell, Hopfinger's consideration would seem to negate specific recognition in solution, whether described by the electrostatic potential or any other descriptor. One therefore concludes that the interaction pharmocophore, like any other descriptor of "complementarity" or "matching site interaction" or "matching hydrogen bonding site" (such as between DNA base pairs), refers to the situation in which the outcome of the interaction iE decided by (possibly instantaneously) "naked" sites. This does not assume total desolvation of the molecules but rather the formation of a common solvent-cage for the molecules interacting through specific sites (2). Even at this stage, the distance between the centers of interacting sites would still be large enough (e.g., 4 A) for the electrostatic potential to serve as a reactivity criterion. This would be the final, and probably the decisive step in the interesting "path of unique interaction pharmocophores" suggested by Hopfinger (3). It is obvious from simple theoretical considerations that only the initial steps of the interaction can be represented by the interaction pharmacophore approximation. Other theoretical methods are continually being developed and applied selectively to appropriate problems of drug-receptor interactions. For example, we have constructed a formalism to investigate the consequences of molecular interactions occurring simultaneously at several sites, by studying the pertubation of the reactivity and of the electronic structure following such interactions (66)(67)(68). Applied to the study of multimolecular complexes-such as the simultaneous interaction of adenine with several water molecules (68) or the interaction between 5-HT and imidazolium cation (69)-applications of this formalism revealed the nature of a possible mechanism for the "epreferential activation of (certain) secondary sites when a large molecule is already interacting at a primary active site" (68). We have discussed the importance of such reactivity induction mechanisms in drug-receptor interactions for muscarinic cholinergic agonists (52) and antagonists (54) and also for the mechanism of activation of the histamine H2-receptor (57). In a separate section below we describe the essential procedures related to the calculation of these multiple pertubation corrections and their use in the characterization of molecular activity.
Calculation of Electrostatic Potentials and Force Vectors: General Considerations. Classical electrostatics define the potential generated at any point in space by a given charge distribution. According to the well-known Hellman-Feynman theorem (70), the charge distribution calculated from the molecular wave function gives the description of the electronic charge distribution of a molecule as a classical distribution of charges. Consequently, it is possible to obtain the electrostatic potential at any point in space from the continuous electronic charge distribution and the nuclear charges of a molecule. This description of the electrostatic potential is independent of the quality of the wave functions and is therefore not related to the assumptions inherent in the generalized Hellman-Feynman theorem regarding the cancellation of forces on the nuclei. The electronic potential is also the first-order energy (El) in the pertubation expansion of the interaction between the molecule and a point charge. It therefore represents the energy of the interaction between the point charge and the molecular charge distribution before the latter is perturbed by the interaction. Truhlar and collaborators have shown (49)(50)(51) that this physical interpetation is borne out by scattering experiments with high energy electrons. The electrostatic potential is measured in these experiments as the "static" part of the potential surface from which the electrons are scattered. Measurements agree with the results of calculations using the formulation given above (48).
The major obstacle in the calculation of accurate MEP maps is an economical algorithm for computation of this property at many points on surfaces relevant to molecular reactivity. An analytic formulation based on the Poisson equation of electrostatics yielded an efficient procedure for the computation of MEP maps (71,72). Taking advantage of the fact that most current LCAO methods for the calculation of molecular wave functions use gaussian expansions of the basis sets, this procedure makes possible the calculation of electrostatic potentials on extensive surfaces for large molecules (47). However, the choice of surfaces relevant to molecular reactivity needs attention, as discussed below.

New Approaches to the Calculation of MEP and the Analysis of Its Characteristics
The Choice ofMolecular Surfaces. To represent a property of a compound relevant to its potential toxicity, the MEP has to be calculated and analyzed in a manner that will elicit the best representation of the reactivity characteristics of the molecule. For example, the molecular reactivity that is responsible for the interactions causing the toxic effects, either directly or via transformation into a potentially toxic compound, could be related to a strictly localized area of the molecule (e.g., attack by a charged species in electrophilic or nucleophilic reactions); it would then suffice to examine the molecular reactivity expressed by the MEP in that localized area alone. If, however, the toxicity of the molecule is related to its ability to act as a whole, e.g., by intercalation,then a more relevant representation of the MEP would be on a surface around the molecule. These two simple examples illustrate the manner in which MEP maps can be used, and have been used in the past [e.g., see (57) for localized values, and (24,47,48) for values on surfaces] to obtain information about the reactivity characteristics of drugs in relation to specific mechanistic propositions. As a screening procedure, however, this approach is strongly dependent on information about the mechanism and is therefore less general than required for a selection of MEP characteristics relevant to broadly defined toxicity. A methodology required for computing, displaying and analyzing the MEP on molecular surfaces relevant to a wide range of molecular interactions is clearly needed.
One such surface was defined by Richards (73,74) as the surface accessible to a water molecule. In practice, the probing water molecule is usually approximated by a sphere with an appropriate radius (14). An example is shown in Figure 3, in which the molecular surface surrounding the molecule of the neurotransmitter 5hydroxytryptamine is viewed from six different directions that represent the six faces of a box of enclosure. A potential problem with this approach is the possibility that due to the finite size of the probe, molecular enclaves will be formed which will not be sampled. This caveat is especially important in conformationally flexible molecules in which such molecular enclaves can often be exposed by a simple conformational change.
There are, however, ways of changing the molecular surface on which the MEP will be calculated. One is to change the size of the sampling probe, and the other is to change the distance at which the probe is allowed to approach the atoms that constitute the molecule. The two modulations do not have the same consequence for the investigation of molecular properties: changing the size of the probe results essentially in a differently shaped surface, but at the same distance from the molecule. Increasing the size of the probe will tend to smooth the surface so that the transition between atoms will be less pronounced; however, more molecular enclaves become possible as the size of the probe increases. Decreasing the size of the probe eliminates the problem of enclaves but may results in a rugged surface. This mode of modulation is equivalent to changing the steric properties of the reagent with which the molecule interacts. On the other hand, changing the distance at which the probe is allowed to approach the molecule introduces relatively small changes in the shape of the surface but the resulting surface is positioned at a different distance from the molecule and therefore the molecular properties, such as MEP, calculated on this surface will have different values. For example, increasing the allowed distance will reduce the values of the MEP and will reduce the contrasts between the different regions on the molecular surface. Both modes of modulation of the computed molecular surface should be investigated in order to determine the sensitivity of the MEP to such changes, and in order to identify the most relevant ap-  proach for the prediction of molecular reactivity related to toxic effects.
Display of MEP: Considerations and Methods. Closely related to the calculation of the MEP on a molecular surface is the questions of its display. Here, the major methodological problem is the display of a nonuniform property on a three-dimensional surface. Several approaches can be delineated to display the maximal amount of information which will be conducive to quick and meaningful analysis. The first approach to the display of the MEP consists of generating isoenergetic contours on the three-dimensional surface. In the past we have used the method of isoenergetic contouring of MEP planes that were calculated at some fixed orientation with respect to the molecule (47). A complication in isoenergetic contouring on three-dimensional surfaces is the isoenergetic interpolation between calculated points. The interpolation on planes assumes a polynomial (usually quadratic or cubic) dependence of the potential on the cartesian coordinates. This assumption can also be applied to the MEP on the three-dimensional surface, but the interpolation has to be constrained by the shape of the surface in the vicinity of the interpo-lated point. This constraint does not usually present difficulties on regular surfaces (e.g., spherical sections that surround an atom), but surfaces with regions of inflection are more complicated and therefore more difficult to deal with. One can adopt an operational scheme in which the contour can be calculated on regular surfaces using the usual constraints and then joined over the problematic regions. Since these regions are often small, the approximation is likely to be very good. The major problem with displays of three-dimensional surfaces and the isoenergetic contours on them is the fixed point of view and the resulting hidden portions of the surface. Thus, no single viewpoint can provide a display that affords a complete inspection of the MEP on the molecular surface. Since several viewing orientations are usually required, one needs to rotate the molecule with its associated surface and redisplay the information in the new orientation.
The second approach consists of a projection technique in which the MEP calculated on the molecular surface will be projected on a "reporter" plane (75). This approach has several advantages. The projection of the MEP on a plane eliminates the problem of constrained interpolations as described above and allows the use of simple modes of representation: either the usual contouring technique (60) or the recently implemented threedimensional display in which the third dimension, perpendicular to the plane, is used to represent the values of the potential. The use of six "reporter" planes is achieved by enclosing the molecule in a box, and using the six faces of the box as the "reporter" planes, in the same way as described for the molecular surfaces in Figure 3, the third dimension (height) being the value of the electrostatic potential.
Analysis of MEP. The two major aspects of the analysis of MEP are its characterization and the computation and characterization of properties derived from the MEP. Implicit in the analysis is the problem of display (see above) because the MEP depends on the choice of surfaces on which it is displayed, as well as on the distance from the various parts of the molecule.
CHARACTERIZATION OF MEP. The characterization of MEP can be divided into two stages: the localization of extrema on the surface and the formulation of a recognition pattern which is discriminant on the molecular level and therefore might be related to the toxicity of the molecule. The identification of local extrema is performed by a search through the complete set of MEP values on the molecular surface. The identification of a local extremum is usually simple because it has to be surrounded by values that are higher than the value at that point (for a local minimum), or values that are lower for a local maximum. This search, performed on the entire surface, also yields the absolute extrema of the MEP.
In the past, we have used this procedure to identify local minima in MEP generated on planes positioned in parallel to the molecules; e.g., the results from this procedure were used to identify directionality vectors in relation to recognition elements in drugs that bind to the 5-HT receptor, as discussed above. Drugs such as 5-HT (24), LSD and Br-LSD (47,48), and various neuroleptics such as spiroperidol, haloperidol, benperidol, and azaperone (76) were analyzed by this procedure. We have also applied this procedure to the identification of local minima on electrostatic interaction surfaces (not MEP) between two molecules (47), e.g., the electrostatic interaction between 5-HT, or 6-HT, and imidazolium cation (47).
The identification of local extrema is the first step in the characterization of the MEP. In our studies we have shown that the two local minima in the MEP of 5-HT and related compounds determine the directionality of the electrostatic orientation vector that is related to the recognition of these two molecules at the 5-HT receptor (47). As discussed above, this vector and its spatial relation to the positively charged nitrogen in these molecules serve as discriminant properties by which the drugs are specifically recognized at receptors and also identify their orientation in molecular complexes in solution (77,78). Such spatial relationship between extrema in the MEP are related to the molecular structure and to the charge distribution in the molecule. They also are the representation of the way in which the molecule presents itself towards a site with which it interacts, e.g., the biological target of a toxic xenobiotic. DIRECTIONALITY OF MEP. Beyond the simple characterization of maxima interacting with negative sites and minima with (relatively) positive sites the disposition of these extrema in space is imposing a geometrical constraint that defines the orientation of the molecule in a nearly homogeneous electrostatic field such as that produced at relatively large distances from charged centers. In molecular interactions, this geometrical constraint rests on the requirements of maximization of the interaction between the molecule and a complementary site in the reagent (e.g., receptor, membrane recognition site, etc.). It is therefore useful to characterize the dispositions of local extrema in terms of the commonly accepted physical nomenclature including dipolar and homopolar directionalities. For molecules with distinct separation between a maximum and a minimum in the MEP, a dipolar directionality with respect to an internal reference in the molecule is an appropriate characterization. This can also be quantified according to the values of the potential and the distance between the extrema. The dipolar directionality can be compared to the dipole moment of the molecule. Molecules with no distinct separation between maxima and minima (e.g., one type of the extrema does not exist) but separation between several local extrema of the same type can be characterized by a homopolar directionality. This can be derived by connecting the local extrema by a straight line and specifying the directionality with respect to an internal coordinate system in the molecule. (A more rigorous algorithm is described below.) It is assumed that the local extrema that will be used for such homopolar directionality will be approximately positioned on the same side of the molecule. In terms of the "reporter" planes used to display the MEP calculated on three-dimensional surfaces (see above), the homopolar directionality is expected to be either on the same "reporter plane" or on adjacent planes. The importance of such homopolar directionalities is that they avoid the need for a specific reference to a molecular skeleton when different molecules are compared. In fact, alignment of the homopolar directionalities regardless of molecular structure can provide an indication of how molecules that are structurally different can exhibit the same reactivity and show similar patterns of toxicity, an exceedingly important capability for the evaluation of potential toxicity of chemicals before they are synthesized and tested [for a related example see the comparison of 5-HT and LSD (79)].
THE GEODESIC AS A RIGOROUS DESCRIPTOR OF DI-RECTIONALITY. Given the importance of the directionality of the MEP for the characterization of molecular reactivity, there exists a need for a rigorous definition of directionality vectors. In our experience with the calculation and characterization of MEP we found that a simple connection of the local minima by a straight line may sometimes require passage through a region that is sterically obstructed or has a very different nature from the local minima (e.g., connecting two minima through a maximum). The general directionality is then better defined by a geodesic, which is the shortest path between two points, e.g., two extrema, on a given surface. This definition is subject to the requirement that the absolute value V of the average MEP on this path is maximal. The average MEP can be defined in the following way: Let s be the path between two extrema and V, the value of the MEP along this path. The average MEP is then: This constraint has the advantage that it selects a meaningful geodesic as well as provides an algorithm for calculating the path s.
The physical significance of this path for molecular interactions can be understood from the following system. When a line of positive charges is distributed on a given surface above a molecule, the absolute value of the energy of interaction of these charges with the molecule will be maximal if they are distributed along the geodesic as defined above, because in this orientation the charges will be positioned along the greatest values of MEP when compared to any other orientation. Thus, this geodesic defines a general concept of the directionality of the MEP that is related to optimal molecular interactions. This concept can be easily generalized to MEP which have only one minimum by defining a certain length of the path s between two arbitrary points chosen according to an alternative reactivity criterion. THE GRADIENT OF MEP (THE ELECTRIC FIELD). Another important property in the analysis of MEP is the gradient of the MEP. The gradient of the MEP at each point on the molecular surface represents the electric field at this point. There exist two alternative ways of representing the gradient of the MEP. One is to draw the electric field lines that converge to the local minima (76). An alternative way of displaying the gradient is to compute it at the same points as the MEP and then display it in a vectorial mode. The orientation and magnitude of the force vector at each point can be represented by an arrow whose direction will indicate the direction of the force acting on the point charge and will have a length proportional to the magnitude of the force (80).
Higher Order Contributions to Molecular Interactions Obtained from Pertubation Expansions. Many of the chemicals that are toxic and exhibit a variety of biological activities have large molecular structures that often incorporate several reactive sites. These sites could compete for specific interaction at biological targets (e.g., similar chemical bonds competing for metabolizing enzymes) or could reinforce each other's action by interacting at complementary sites in membranes, structural proteins, or genetic material (e.g., intercalating portions preparing alkylation site for attack of a biological target). In our work on the molecular mechanisms of drug-receptor interactions we have found that some sites in the drug molecule are activated by a primary interaction of another group; this was the case of the side chain in acetylcholine (52,54), and in 5-hydroxytryptamine (77)(78)(79) of the cationic amine in morphine (59), and the triggering of the activation mechanism of histamine H2 receptors (57). Theoretical methods are therefore needed to deal with criteria for multiple site reactivity. In addition, since the molecules could be dissolved in a solvent which is itself active, a useful theoretical approach must make it possible to deal effectively with the consequences of solvation at one site on the interaction with a second active site.
In dealing with structure-activity relationships of chemicals we would like to associate the activity with the properties of the studied molecule and thereby define criteria that are exclusively dependent on stereoelectronic characteristics as defined above for the reactivity characteristics derived from MEP. An excellent approach is perturbation theory, where the various stages of the interaction can be calculated separately. The total interaction energy can be broken down into physically identifiable components, and the changes caused by the interaction in the molecular wave function, the electronic density, and the MEP, can be calculated. We have developed such a multiple perturbation formalism within the frame of uncoupled Hartree-Fock formalism (66)(67)(68). In this formalism, the multiple perturbations (e.g., the effect of the biological target on the xenobiotic, or the effect of the environment) are represented by an electrostatic operator of the type Q/(rro), which assumes that the main interactions are electrostatic and can be represented by discrete charges (Q). The total energy of the system is separated into the energy of the unperturbed molecule and the energy of interaction. The energy of interaction in turn consists of the electrostatic interaction of the perturbing point charge with the nuclei and with the electrons. The electronic interaction energy of the xenobiotic with its target can then be expanded in a double perturbation expansion representing the multiple effects of the molecules on each other; each term in this expansion has a clear physical meaning (66). The firstorder terms in the expansion (i.e., E1 (, and E,1) represent the interaction of a point charge with the static charge distribution of the molecule and are therefore equal to the electronic component of the electrostatic potential described above. The second-order terms (i.e., E2,0, EO,2 and E1l1) represent the interaction of a point charge with a perturbed charge distribution, i.e., of the target interacting with the xenobiotic it already modified. In the case of E2,0 (or Eo,2), the interaction is between the point charge (i.e., the target) and the perturbed charge distribution that it induced. Therefore, this term is equivalent to the polarization energy that the point charge induced in this molecule. The meaning of E,,1 is more complex because it represents the interaction of a point charge with a perturbed charge distribution induced by another (e.g., another component of the target or the environment); this term represents therefore the change in reactivity of a molecule produced by the presence of a perturbing interaction. The physical meaning of these terms makes them useful in assessing the degree of molecular polarizability and induced reactivity, properties that have often been connected specifically to carcinogenicity and toxicity (81,82) and can be calculated explicitly within this formalism.
Using these terms we showed that if the attacked molecule (e.g., the biological target or the solvent) can be represented by a set of point charges that reproduce some of its physical properties (such as the dipole and quadrupole moments, or the total electrostatic potential), then we can use the same multiple perturbation formalism to describe interactions with neutral molecules, such as water (53). Our more recent improvements make possible the specific treatment of molecular interactions without need for the point charge expansion (67,68). We have also generalized this formalism with diagramatic perturbation theory (68), and have shown that the pairwise treatment of the interaction between any number of molecules is well represented within the framework of that approximation.
In other studies we examined the relation between quantitative SAR (QSAR), molecular structure, and biological activity (24,48,77,78), and showed that reactivity characteristics obtained from MEP and from the higher-order terms (El l and E2,0) provide a link between the various forms of SAR studies. For example, in the case of drugs acting at the biological receptor of the neurotransmitter 5-hydroxytryptamine, maps ofthe polarizability (ES ,,) indicated the molecular regions that are most readily polarizable and are therefore most likely to contribute to the stabilization of electron donor-acceptor complexes (77,83). The relation between the E2,0 values and the localization of the highest occupied molecular orbital (HOMO) of all the tryptamines derivatives that we calculated ab initio showed that HOMO has always a disproportionately large contribution to the toal ES,(, value and that the polarizability calculated from HOMO coincides with its sites of localization (77,78). These results are compatible with the hypothesis that intermolecular forces enhanced by polarization (such as in electron donor-acceptor complexes), are involved in the interaction of tryptamines with their receptor, and support the use of HOMO localization as a reactivity criterion in cases where such a mechanism of interaction with a biological target is suspected. We therefore made use of E2,,( maps in the analysis of changes in the degree and localization of site polarizability in the various molecules we studied (84).
The applications of these methods to the study of criteria for toxicity lie in their ability to represent the effects of sequential and secondary interactions on the activation of new reactive sites. By simulating the effects of charge redistribution in a compound after its interaction with hypothetical reactants or biological targets such as membrane fragments or active sites of enzymes, this approach makes possible the consideration of very complex mechanisms for the action of xenobiotics. It could bring the theoretical study of discriminant criteria for toxicity closer to the description of the possible steps of the molecular mechanisms that determine the interaction of the chemicals with the biological targets where the toxic effects are produced; notably, it could even provide a representation of the initial steps in the enzymatic mechanisms that convert harmless xenobiotics into toxic metabolic products.