Molecular Modeling for Screening Environmental Chemicals for

Aug 31, 2009 - E-mail: [email protected]., † ... Anders Olofsson , A. Elisabeth Sauer-Eriksson , and Patrik L. Andersson ... John D. Sivey , ...
9 downloads 0 Views 282KB Size
Download PDF
1594

Chem. Res. Toxicol. 2009, 22, 1594–1602

Molecular Modeling for Screening Environmental Chemicals for Estrogenicity: Use of the Toxicant-Target Approach James R. Rabinowitz,* Stephen B. Little, Susan C. Laws,† and Michael-Rock Goldsmith‡ National Center for Computational Toxicology, U.S. EnVironmental Protection Agency, Research Triangle Park, North Carolina ReceiVed April 8, 2009

There is a paucity of relevant experimental information available for the evaluation of the potential health and environmental effects of many man made chemicals. Knowledge of the potential pathways for activity provides a rational basis for the extrapolations inherent in the preliminary evaluation of risk and the establishment of priorities for obtaining missing data for environmental chemicals. The differential step in many mechanisms of toxicity may be generalized as the interaction between a small molecule (a potential toxicant) and one or more macromolecular targets. An approach based on computation of the interaction between a potential molecular toxicant and a library of macromolecular targets of toxicity has been proposed for preliminary chemical screening. In the current study, the interaction between a series of environmentally relevant chemicals and models of the rat estrogen receptors (ER) was computed and the results compared to an experimental data set of their relative binding affinities. The experimental data set consists of 281 chemicals, selected from the U.S. EPA’s Toxic Substances Control Act (TSCA) inventory, that were initially screened using the rat uterine cytosolic ER-competitive binding assay. Secondary analysis, using Lineweaver-Burk plots and slope replots, was applied to confirm that only 15 of these test chemicals were true competitive inhibitors of ER binding with experimental inhibition constants (Ki) less than 100 µM. Two different rapid computational docking methods have been applied. Each provides a score that is a surrogate for the strength of the interaction between each ligand-receptor pair. Using the score that indicates the strongest interaction for each pair, without consideration of the geometry of binding between the toxicant and the target, all of the active molecules were discovered in the first 16% of the chemicals. When a filter is applied on the basis of the geometry of a simplified pharmacophore for binding to the ER, the results are improved, and all of the active molecules were discovered in the first 8% of the chemicals. In order to obtain no false negatives in the model that includes the pharmacophore filter, only 8 molecules are false positives. These results indicate that molecular docking algorithms that were designed to find the chemicals that act most strongly at a receptor (and therefore are potential pharmaceuticals) can efficiently separate weakly active chemicals from a library of primarily inactive chemicals. The advantage of using a pharmacophore filter suggests that the development of filters of this type for other receptors will prove valuable. Introduction The potential risks due to exposure to man-made chemicals in the environment must be evaluated in order to protect human health and the environment. Little or no experimental information is available about the potential biological effects of as many as 75% of these chemicals (1). In addition, the data that are available are often insufficient to adequately evaluate the potential hazards of each of these chemicals. Therefore, there is a compelling need to develop information that would enable the screening of the potential health and environmental effects of large numbers of man-made chemicals (2). While animal studies have historically been preferred for risk assessment, such data sets are often difficult and expensive to obtain, and yet still require significant extrapolation to be applied to this task. The National Academy of Science report, Toxicity Testing in the Twenty-First Century presents a vision of a new paradigm * Corresponding author. E-mail: [email protected]. † S.C.L. is a member of the Reproductive Toxicology Division, NHEERL, U.S. EPA. ‡ During this study, M.-R.G. was a cross Office of Research and Development postdoctoral fellow working in the NCCT. He is presently at the National Exposure Research Laboratory, U.S. EPA, Research Triangle Park, NC.

10.1021/tx900135x CCC: $40.75

for toxicity testing (3). This vision is based on the identification and characterization of toxicity pathways. The capacity of specific chemicals to participate in these pathways would then be interrogated using rapid (in vitro or perhaps computational) test methods. Integration of the data using informatics approaches that consider the causal nature of these pathways could then be applied for the evaluation of the risks posed by a specific chemical. In support of this vision, efforts are underway to devise approaches that combine more rapidly obtained experimental data and data computed from chemical structure with informatics approaches to perform toxicity evaluations (4). In the initial stage, these approaches may be seen as extrapolations from more readily obtainable data through the more traditionally applied animal data to the evaluation of potential health and environmental effects. Until more experience is obtained and success is demonstrated, these approaches should be considered as methods to screen chemicals and develop testing priorities. After sufficient experience, the support of this approach through the more traditionally accepted animal data may become unnecessary and this type of evidence used directly to evaluate chemical toxicity.

This article not subject to U.S. Copyright. Published 2009 by American Chemical Society. Published on Web 08/31/2009

Molecular Modeling: Screening for Estrogenicity

Knowledge of the potential mechanisms of toxicity provides a rational basis for extrapolation from chemicals where there is a great deal of data on toxicity to chemicals for which little data exists. The differential step in many mechanisms of chemical toxicity may be generalized as the interaction between a small molecule (a potential toxicant) and one or more macromolecular targets. The small molecule may be the chemical itself or one of its descendents. Modeling the potential of a molecule for specific interactions of this type is a source of insight into its potential toxicity. In a previous paper, the application of an approach based on computation of the interaction between a potential molecular toxicant and a library of macromolecular targets of toxicity was proposed for chemical screening (5). In order to use a library of this type to assess the potential for untested chemicals to be toxic and to determine the most likely pathways for toxicity, a rapid method to evaluate interactions between the molecule and the target is needed. Molecular docking (6-8) has been developed to screen large libraries of chemicals for molecules that interact with specific sites on proteins and therefore are potential pharmaceutical agents (9-11). It also has been used to identify xenoestrogens (12) but has been infrequently applied to investigate the potential toxicity of weaker agents (5, 13). It is a rapid method that depends on a surrogate for the interaction energy (a scoring function) that is determined heuristically from data determined from many interactions between small molecules and proteins. In the current study, this computational approach is applied to a set of environmentally relevant chemicals interacting with the estrogen receptor (ER1) and the computational results compared to results from a recent experimental study that determined the relative binding affinities (RBA) of these molecules to the rat (ER) in a complex preparation (14-16). In the experimental study, the RBAs were determined by the ability of the natural ligand, to bind to the estrogen receptors in a rat tissue preparation in the presence of increasing amounts of the test chemical. This assay is currently being evaluated for inclusion in a test battery to be used for determining the potential of environmental chemicals to disrupt endocrine function (17, 18). It simulates the initiating step of a process that results in estrogen receptor mediated physiological responses. Chemicals that are active in this assay may have tissue dependent agonist and/or antagonist properties. Measuring the response as a function of the concentration of the test chemical allows true binders to be separated from chemicals that interfere with the binding of the natural ligand through other mechanisms. These results are analogous to what is modeled in molecular docking computer experiments and as such form an excellent system for exploring the use of computational docking as a tool in a toxicity screen. The data set has the additional advantage that most of chemicals that bind are found to bind weakly, and docking methods can be tested for their capacity to find chemicals that bind weakly to a receptor and are not drug like. 1 Abbreviations: RBA, relative binding affinity; ER, estrogen receptor; KIERBL, the label in DSSTox (15) for the primary experimental data set used in this study (14); IC50, concentration of a test chemical that inhibits the maximal specific binding of 0.33 nM radiolabeled (3H) 17-β-estradiol (E2) to rat uterine cytosolic ER by 50%; Ki, inhibition constant for the test chemical, i.e., micromolar concentration that will bind to half the ER at equilibrium in the absence of radioligand or other competitors; TPR, true positive ratio, the ratio of the positives discovered divided by the total number of positives; FPR, false positive ratio, the ratio of the number of false positives divided by the number of negatives; HPTE, 2,2-bis-(phydroxyphenyl)-1,1,1-trichloroethane.

Chem. Res. Toxicol., Vol. 22, No. 9, 2009 1595

Materials and Methods Molecular modeling software tools, which have been developed for the discovery of novel pharmaceutical agents, were applied to identify chemicals that bind weakly to the rat estrogen receptor. These tools attempt to order a set of chemicals by their capacity to bind to a specific site in a protein (7). This is accomplished by computationally docking each chemical into the ligand binding pocket of the receptor. In molecular docking, the most favorable interaction between the chemical-binding pocket pair (toxicant-target in this case) is identified and quantified. The protein targets used in this study were derived from crystallographic structures of the estrogen receptors, each with a ligand bound. The atoms of the ligands were computationally removed from crystal structures to create computational targets for docking. In order to compare the results from this computational study to the experimental results from the cytosolic preparation (14, 16), four targets were created. Multiple targets were required because the rat uterine cytosol used in the experimental study contained both ERR and ERβ, and the experimental method does not distinguish between agonist-like and antagonists-like binding. Crystal structures were not available from the Protein Data Bank (PDB) (19) for the rat ERR. Therefore, the rat ERR targets were derived from human ERR (20, 21) by homology models. There is 96.6% homology between the ligand binding domain of the rat ERR and the ligand binding domain of the human ERR. The homology models were developed in Molecular Operating Environment substituting the ligand binding domain of the rat receptor sequence for the similar human receptor sequence (22-24). The model targets derived in this way were subsequently optimized by molecular mechanics using the AMBER94 force-field (25). Crystal structures of the rat ERβ are available directly from the PDB for both agonist (26) and antagonist (27) cocrystallized conformations. The AMBER94 force-field was also used to optimize the structures of these targets. In this manner, four targets were created, ERR (agonist), ERR (antagonist), ERβ (agonist), and ERβ (antagonist). A structural data file containing each of the 281 chemicals included in the experimental study (14) was obtained from DSSTox (16). The label for this data set in DSSTox is KIERBL. The structures were imported into the Molecular Operating Environment software (22) and the molecular geometries were optimized using the MMFFx force field (28). For the purposes of this study, a chemical was considered active (e.g, a true competitive binder for the ER) if it had an experimentally determined Ki and an IC50 less than 100 µM in the experimental study (14). Thus, of the 281 chemicals, 15 were classified as active ER binders; 266 were classified as inactive for this study. Of those classified as inactive for this study, the experimental study left some uncertainty about 11 of those chemicals. Of those 11 chemicals, 9 had limited competitive binding curves because of poor solubility and had IC50 values that were uncertain in the experimental study but were greater than 100 µM, while the other 2 chemicals had partial competitive binding curves that were deemed sufficient to obtain IC50 values but those values were greater than 100 µM, more than 5 orders of magnitude greater than the natural ligand. All chemicals used in this study are reported in the U.S. EPA’s Distributed StructureSearchable Toxicity (DSSTox) Public Database Network where all experimental results and a Structure Data File with chemical identifiers are available for public access (15, 16). Two different molecular docking approaches, using two different software packages were employed. The first package, eHiTS (29, 30) relies on a decomposition of the optimized structure of potential ligands into fragments. Then individually docking each fragment into the binding site, obtaining multiple fragment poses and using a graph matching scheme to reconstruct the structure of the potential ligand based on connectivity tables. The reconstructed molecular poses are subsequently locally optimized, and a score (a surrogate for the binding energy) for each pose is obtained from a heuristic scoring function (29). The second algorithm, FRED (31, 32), applies a systematic exploration of the rotational and translational space of each putative ligand in the target space to

1596

Chem. Res. Toxicol., Vol. 22, No. 9, 2009

Rabinowitz et al.

Figure 1. Idealized example of the number of chemicals as a function of docking score. In this example the scoring function imperfectly separates the positive chemicals from the negative chemicals. A line Q-R is drawn and used to make a prediction of the activity of each chemical. The predictions made on the basis of the docking scores are a function of the position of Q-R as it moved horizontally.

discover the configurations for the putative ligand-target interaction with the best score. It has been shown to be sufficiently accurate and rapid for virtual screening for leads for potential pharmaceuticals (33, 34). Both of the docking methods used in this study have been shown to be successful at finding drug-like (molecules that bind strongly to the receptor) estrogenic molecules (29, 33). Their capacity to separate weak binders from nonbinders is a primary goal of this computational study.

Results The experimental results of Laws and co-workers (KIERBL) (14) provide an excellent data set for investigating the application of the toxicant-target approach for the prioritization of chemicals for potential toxicity. The advantages of comparing computational molecular docking results to these experimental results are as follows: first, there are a number of excellent crystal structures of both R and β estrogen receptors available in the Protein Data Bank (20, 21, 26, 27) and others) that can be used to synthesize macromolecular targets for computer docking; second, a library of 281 chemicals was tested in the same laboratory with the same protocol for their capacity to compete with radiolabeled 17-β-estradiol for their binding to the rat ER; and third, the experiments yield a relatively direct measurement of what is modeled in computational docking, the energy of interaction between the test chemicals and the receptor compared to the energy of interaction of the receptor with 17β-estradiol. In the current study, each of the 281 chemicals was computationally docked into 4 rat ER targets. The large structural differences between targets obtained from cocrystallization with an agonist or an antagonist chemical were verified by observation (results not shown). The targets created from the agonist crystal structures are smaller and enclosed, while the targets created from the antagonist crystal structures are larger, more open, and will accommodate larger molecules. The best scores for each molecule interacting with each target are reported for each computational approach and used for further analysis in this study. The chemical names and the experimental determination of binding from the criteria described in the Materials and Methods section and the measured values along with the docking scores for the interaction of that

chemical with each of the four targets using both methods are in Table 1 of the Supporting Information. For this data set, each chemical that binds (with the exception of 17-β-estradiol) has an experimental Ki that is 3 to more than 5 orders of magnitude greater than that of 17-β-estradiol. The rat uterine cytosol used as the source of receptors in the binding assays contained both the R and β subtypes, with the predominant type being ERR (35). Models that attempted to realistically combined ERR and ERβ scores were used, but they did not significantly improve the results when compared to those from a consideration of ERR alone (data not shown). Thus, in the remainder of this analysis only the scores for ERR with both molecular docking methods are used. Since the experimental results do not distinguish between molecules that bind like an agonist and molecules that bind like an antagonist, scores for both targets have been combined to give a composite score for each chemical. The composite score for each chemical is chosen as the score that indicates the most stable interaction for that chemical. Comparison of the results for agonist and antagonist targets might indicate whether the chemical was likely to be an agonist or an antagonist. The composite score for each chemical was considered along with the separate scores for the agonist and antagonist targets. The results from computational docking are scores (a surrogate for the interaction energy) for each potential ligand-protein pair. The scores are determined by the internal scoring functions of the computational method used. The parameters in the scoring functions have been determined from multiple experimental measurements of protein-small molecule structures (30, 32). They are not adjusted to optimize the results for comparison to these specific receptors or experimental results. For the purpose of predicting potential activity, a demarcation between likely active and likely inactive chemicals must be chosen. An example of this demarcation is shown in Figure 1. The only adjustable parameter in this approach is the choice of demarcation between the predicted likely active chemicals and the predicted likely inactive chemicals (line Q-R in Figure 1). In Figure 2, the predicted true positive ratio (TPR) is plotted as a function of the predicted false positive ratio (FPR) as the position of the demarcation line between predicted positive and negative (Q-R) is varied. Table 1 shows a synopsis of these results for ERR

Molecular Modeling: Screening for Estrogenicity

Chem. Res. Toxicol., Vol. 22, No. 9, 2009 1597 Table 2. Summary of the Results for Each Approach When 2 Constraints from the Pharmacophore Are Applieda number of true positives identified

5

10

14

15

FRED agonist FRED-antagonist FRED-composite eHiTS agonist eHiTS antagonist eHiTS composite

6 6 7 7 10 10

14 18 15 14 19 20

22 (8) 30 29 20 41 29

39 (14) 37 (13) 48 (17) 23 (8) 42 (15) 32 (11)

a Each entry is the number of chemicals that would be called positive when Q-R is adjusted to yield the number of true positives heading that column. For example, when using the eHiTS method and the rat ERR agonist target to obtain the first 5 true positives, 7 chemicals are called positive. In that case, there are 2 false positives, 10 false negatives, and 264 true negatives. The numbers in parentheses are the percentage of the database that number represents. For the last column, there are no false negatives.

Figure 2. (a) Docking results using the method FRED with no constraints. The true positive rate (TPR) is shown as a function of the false positive rate (FPR) for agonist, antagonist targets, and a composite score constructed from selecting the best score of the two for each chemical. They are compared to random selection. (b) Docking results using the method eHiTS with no constraints. The TPR is shown as a function of the FPR for agonist, antagonist targets, and a composite score constructed from selecting the best score of the two for each chemical. They are compared to random selection.

Table 1. Summary of the Results for Each Approach without Constraints Applieda number of true positives identified

5

10

14

15

Fred agonist Fred-antagonist Fred-composite eHiTS agonist eHiTS antagonist eHiTS composite

16 20 23 8 18 18

51 57 59 18 50 51

76 (27) 105 104 36 95 68

153 (54) 109 (39) 119 (42) 46 (16) 97 (35) 74 (26)

a Each entry is the number of chemicals that would be called positive when Q-R is adjusted to yield the number of true positives heading that column. For example, when using the eHiTS method and the rat ERR agonist target to obtain the first 5 true positives, 8 chemicals are called positive. In that case, there are 3 false positives, 10 false negatives, and 263 true negatives. The numbers in parentheses are the percentages of the data base. For the last column, there are no false negatives.

targets for both docking methods. Using the eHiTS method, we observed that the results for the agonist target are superior to the results for the antagonist target and the composite score. Using the FRED method, we observed that the results for the

agonist target are superior to the results for the antagonist target and composite score for identifying 14 of the 15 active chemicals, but finding one of the active chemicals (4,4′sulfonyldiphenol) is more difficult; that chemical is best found by the antagonist target. All of the active chemicals appear in the first 16% of the data set when the eHiTS, fragment based approach is used with the agonist target and 39% (27% when one positive chemical is excluded) when FRED is used. For the preceding results, the best score for each chemical was chosen without consideration of the geometry of binding between the toxicant and the target. Many experimental studies of estrogen receptor binding have shown that specific interactions between the atoms in the ligand and atoms in the receptor play a key role in molecular recognition. As a result, a pharmacophore for binding to the estrogen receptor has been developed (36). In theory, docking algorithms should ensure the proper orientation of the ligand and the target so that the requirements of this pharmacophore are satisfied because the scoring functions should include the interactions that favor the pharmacophore recognition process. However, the development of scoring functions is a continuing process (37, 38), and it has been found that including pharmacophores, in the form of constraints on the allowed docking poses, improves results in some cases (12, 39, 40). In the docking program FRED, constraints of this type may be applied to the geometry of the interacting molecules. The set of chemicals was docked again under the constrained condition that only geometries that produced two appropriate hydrogen bonds between the toxicant and the target (equivalent of the hydrogen bonds made by the hydroxyl group of the A-ring of 17-β-estradiol to an arginine and glutamate in the binding pocket of the receptor) were permitted. The resulting scores using this approach in FRED are shown in Table 1 of the Supporting Information. The plots of TPR as a function of FPR for these results are shown in Figure 3. Table 2 shows these results for various choices of demarcation. All positive chemicals are now found in the first 14% (8% if the same difficult chemical is omitted) of the molecules, and while the results for finding all 15 chemicals is better for the antagonist target than the agonist target, they are quite similar, and the agonist target finds the first 14 more rapidly. The eHiTS program does not yet contain the option of applying externally determined pharmacophore driven constraints directly, but it is possible to use the insight from the results generated by FRED to enhance the results obtained from eHiTS. In the library of possible ligand-receptor poses generated by eHiTS, all poses that do not satisfy the constraints have been eliminated. Those results are also shown in Table 1 of the Supporting Information, and the plots of TPR as a function of

1598

Chem. Res. Toxicol., Vol. 22, No. 9, 2009

Figure 3. Docking results using the method FRED with 2 constraints. The TPR is shown as a function of the FPR for agonist, antagonist targets, and a composite score constructed from selecting the best score of the two for each chemical. They are compared to random selection.

Rabinowitz et al.

is capable of separating chemicals that bind weakly to the estrogen receptor from those chemicals that do not bind to the receptor at all, at least for the example of this data set. In order to explore the activity domain of the approach used in this study, a set of chemicals that are primarily potent estrogen receptor binders was docked into targets constructed from the human estrogen receptor (20, 21). The activities of these chemicals at the estrogen receptors have been measured in a manner similar to that used to obtain the measurements in the KIERBL data set in a single study (41). Chemical information and the experimental determination of RBA (from ref 41) along with the docking scores for the ERR receptor are shown in Table 2 of the Supporting Information. The demarcation lines (Q-R), established previously for the KIERBL data set with each docking method, were then used to separate the active chemicals from the inactive chemicals in this additional data set of primarily strong binders. With the eHiTS approach, all chemicals that compete with E2 for the estrogen receptor are identified by the agonist target with the exception of coumestrol, methoxychlor (a very weak binder that is not a pharmaceutical-like chemical), and the known antagonists. The FRED method identified all of the active chemicals in this data set except for the known antagonists. The known antagonists are too large for the agonist target but are discovered by the antagonist target. There is only a single false positive with eHiTS and at least two false positives with FRED. Figure 5 shows a comparison of the docking scores for these pharmaceutical-like chemicals and the chemicals in the KIERBL.

Discussion

Figure 4. Docking results using the method eHiTS with 2 constraints. The TPR is shown as a function of the FPR for agonist, antagonist targets, and a composite score constructed from selecting the best score of the two for each chemical. They are compared to random selection.

FPR for these results are shown in Figure 4. Table 2 shows the summary of these results for various choices of demarcation (Q-R). All positive chemicals appear in the first 8% of the molecules. For the eHiTs method with the constraints, the agonist target consistently yields the best results. The addition of a simplified pharmacophore constraint significantly decreases the number of false positives that would result from a demarcation scheme that minimizes the false negative rate for all approaches and eliminates no positive chemicals in this data set. For all of the docking results in this study, molecules that could not be docked were inactive in the experimental study. With the exception of 17-β-estradiol, all of the chemicals in this data set that displace the natural ligand, bind only weakly to the estrogen receptor. The computational methods used in this study were designed to increase the discovery rate of molecules that bind efficaciously to specific protein targets and have been shown to enrich databases for potential pharmaceutical agents by identifying chemicals that bind strongly to receptor targets, including the estrogen receptor (33). Applying the same set of computational tools, an approach has been described that

In a previous paper (5), the long-term goal of developing a library of molecular targets for chemical toxicity was introduced. In order to use a library of this type to assess the potential for untested chemicals to be toxic and to determine the most likely pathways for toxicity, a rapid method to evaluate interactions between the molecule and the (macromolecular) targets is needed. Molecular docking has been used to screen large libraries of chemicals for molecules that interact strongly with specific sites on proteins and therefore are potential pharmaceutical agents. In this study, we evaluate the capacity of two docking methods to discover chemicals that interact weakly (3-5 orders of magnitude less than the natural ligand) with the estrogen receptor. In the KIERBL data of weak binders used in this study, only 5% of the chemicals have any measured capacity to compete with 17-β-estradiol for the ER in the tissue preparation. It is difficult to estimate how many of the untested industrial chemicals are likely to interact with particular biological receptors. However, it is likely that the amount for each receptor is small and that their interactions with the receptors weak when compared to the natural ligand. (Chemicals that interact strongly should be more easily discovered.) For this reason, the experimental data set used in this study is more relevant to environmental circumstances than a data set composed of pharmaceuticals or other chemicals that have been designed to exhibit pharmaceutical-like biological activity. The choice of the ERs as a model target for this exploration of the application of docking as a tool for exploring potential toxicity has added advantages. First, the target is of environmental interest and has specific relevance for endocrine disruptor screening. Second, there are a relatively large number of crystal structures of the ERs cocrystallized with various agonist and antagonist ligands, in the literature. The availability of many similar crystal structures supplies leverage to counteract a short

Molecular Modeling: Screening for Estrogenicity

Chem. Res. Toxicol., Vol. 22, No. 9, 2009 1599

Figure 5. Comparison of docking scores for the KIERBL data set and a set of known strong estrogen agonists (see ref 41). See Tables 1 and 2 of the Supporting Information for scores for all chemicals considered. Only a single active chemical has a score that is above the demarcation determined for the weak binders.

coming of the docking methods used in this study, that is, the protein targets are not considered flexible. The computational targets do not respond to the potential ligand and do not relax to provide a better fit for each ligand. While there are methods that include protein flexibility (42, 43), they are currently too computationally intensive to be considered for this type of application. In examples like the estrogen receptor where there are a number of crystal structures available, flexibility may be included by constructing multiple targets for the same receptor and accepting the best score for each toxicant-target pair. In this study, we have utilized targets for both the agonist and antagonist configuration of the receptor but only a single target for each of these major receptor configurations. Two different docking packages with two different approaches and different scoring functions have been used. The most useful results were obtained when a simplified pharmacophore filter was applied to constrain the allowed docking results. By eliminating toxicant-target poses that do not satisfy a predetermined pharmacophore, the filter significantly reduces the number of chemicals that would be classified as false positives for each approach used, yet its application had no effect on the discovery of true positives in this data set, that are up to 5 orders of magnitude weaker than 17-β-estradiol. The constraints of this pharmacophore based filter could be satisfied by poses generated with 102 of the molecules considered in this study. The pharmacophore filter was derived from decades of experience with the estrogen receptors (36) and is used in this study in a simplified form. It is composed of only two of the three hydrogen bonds found in the interaction between strong binders and the ER. More complex pharmacophores have been proposed (44, 45). There is a balance between minimizing the false discovery rate and discovering all potential positive chemicals that must be carefully considered when increasing the complexity of the filter. The same pharmacophore filter used to discover chemicals that bind strongly to a receptor may not be appropriate in a chemical screen designed to also find chemicals that bind weakly. The pharmacophore filter used in this study contains the geometry for a single hydrogen bond donor interacting with a glutamate and a single hydrogen bond acceptor interacting with an arginine. For other potential receptor targets, the available literature may not be as extensive, and a similar approach for the development of a pharmacophore filter may not be feasible. An alternative approach for the development of pharmacophore filters is to use computational methods

to identify a pharmacophore from crystal structures of the target protein cocrystallized with an array of ligands (46, 47). The pesticide methoxychlor could not fully satisfy even this simple pharmacophore. It contains the necessary hydrogen bond acceptor but does not contain a potential hydrogen bond donor. Experimental evidence suggests that it interacts weakly and in a complex manner with the estrogen receptors (48-50). HPTE (2,2-bis-(p-hydroxyphenyl)-1,1,1-trichloroethane), a metabolite of methoxychlor, does satisfy the filter and is a more active estrogen receptor binder than its parent (48, 49). (See Table 2 of the Supporting Information for the docking results for HTPE.) An even simpler pharmacophore may provide the best approach for including all potential weakly estrogen molecules. The degree of improvement in the results of this study through application of the pharmacophore filter was surprising because the specific interactions characterized by the filter should already be in the scoring functions used in docking. This result suggests that docking methods could be improved (at least relative to their capacity to discover weak binders) by improving scoring functions. Molecular docking with a pharmacophore filter might be incorporated into a tiered or other structured approach for chemical screening (see for instance (51)). The application of a pharmacophore filter could be incorporated into one of the initial tiers, and all molecules that could not satisfy the filter in any docking pose would be eliminated before docking. Each individual pose would still be required to satisfy the filter for scoring during the docking phase, which would be part of a later tier. The docking results obtained from the ERR agonist target are consistently better than the docking results for all other individual targets and the composite scores for the KIERBL data set of weak binders. Part of the explanation is that the tissue preparations in the experimental study contain much more ERR than ERβ (35). All of the active chemicals are discovered most rapidly by using only the agonist target, and the addition of results for the antagonist target by creating composite scores increases the false discovery rate. This may result from something specific about the KIERBL data set (that perhaps it does not contain antagonists). The decrease in the discovery rate when composite scores (constructed from agonist and antagonist results) are used may result from the crystal structures used to create the targets. As a result of a more favorable crystal structure interactions in the antagonist target, there might be a constant favorable offset. This would introduce false positives

1600

Chem. Res. Toxicol., Vol. 22, No. 9, 2009

through the antagonist target when compared to that in the agonist target. The current results suggest that this is the case, and an adjustment of the scores when comparing results from different targets is needed. Clearly, the discovery rates of this study (with the KIERBL data set) would not be changed by including docking at an antagonist target and combining it with the results from the agonist target with an offset, but for data sets that contained (bulky) antagonists, the results would be improved. Most of the molecular identification features are the same for the agonist and antagonist targets, but the antagonist target can accommodate larger molecules. The crystal structure of ERR cocrystallized with hydroxy-tamoxifen was used to create the antagonist target in this study. Hydroxy-tamoxifen is too bulky to dock well into the agonist target. Similar bulky known antagonist molecules were not discovered with only an agonist target when docking the series of pharmaceutical-like molecules. It is likely that a general screen would need both targets. As an element in a tiered or another type of ordered approach for determining chemicals that bind to the estrogen receptor, agonist and antagonist targets might be used separately after prefiltering steps based on various molecular descriptors for each target. The result of the application of a docking tool to the toxicant-target paradigm is a list of chemicals ordered by their score. The score for a molecule is a surrogate for the interaction energy between the potential ligand and the target. These scores are derived from scoring functions that are obtained from consideration of the general form of the interaction between a protein and a small molecule and experimental data describing many protein small molecule interactions. The scoring functions do not depend on the specific chemicals or the specific protein being considered. The score for a specific protein-ligand interaction depends on the geometry of the interaction and the description of the general scoring function. The development of the ordered list is independent of and does not require any data on other chemicals interacting with this specific target. This is different from many other computational methods for evaluating potential toxicity where data for a training set of chemicals interacting with the specific receptor is required. The use of docking tools requires only a crystal structure of the target or a similar protein if homology models are used. The variable line Q-R (see Figure 1) used to dichotomize the set of molecules is the only adjustable parameter in this approach. Optimizing the placement of Q-R allows the determination of the best separation between active and inactive chemicals in the data set and facilitates identification of all the positive chemicals in the smallest subset of molecules. This is a reasonable approach where the necessary relevant data exists, but in cases where that data is not available or in cases where there are other rationales (perhaps resource driven) for choosing the number of chemicals to be tested, the ordered list provides a justification for the order of chemicals to be tested. When chemicals in the initial data set are compared to a large set of known estrogen receptor binders (52), and three physicochemical properties that are relevant for ER binding (but probably not a complete set of relevant properties relevant to ER binding) are computed (53, 54) and employed for comparison, the KIERBL data set is found to be more diverse than the set of known estrogen binders (see Figure 6). To a major extent, this is the case because most of the chemicals are inactive. When only the active chemicals are considered, 4 of the 15 active chemicals are found to be within the physicochemical space of the external strong binding chemicals, 8 are adjacent to that space, and 3 of the active chemicals are outside that space.

Rabinowitz et al.

Figure 6. Comparison of the KIERBL data set of weak binding chemicals to a larger set of known ER binding chemicals. It is shown in the three dimensions of relevant (but probably not complete) computed (22) physical properties, log P (54), total polar surface area (TPSA) (53), and molecular weight.

About 10% of the nonbinding chemicals are contained in the space of external strong binding chemicals.

Conclusions In the primary experimental library considered in this study, only 14 of the 280 (excluding 17β-estradiol) chemicals have measured Ki values and an IC50 values less than 100 µM. Computational molecular docking methods are able to identify all of these chemicals using targets constructed from ER crystal structures. For the best approach applied in this study, the 14 positive chemicals are identified with only 8 false positive assignments. All of the computational molecular docking approaches and comparisons used in this study give reasonable results. All chemicals from the experimental library that could not be docked in the biomolecular targets were inactive. The molecular docking algorithms applied in this study were developed to aid in pharmaceutical discovery. As such, they were designed to find the molecules with the greatest activity. However, the active chemicals in this study (with the exception of 17β-estradiol itself) all have activities 3-5 orders of magnitude smaller than that of 17-β-estradiol and are discovered by these docking algorithms. Ninety-five percent of the chemicals in the library considered in this study are not active. This is not typical of computational studies of environmental chemicals where the chemical library often contains more positive chemicals, and attempts are made to balance the database. This data set is most likely more similar to the problem seen in environmental circumstances where for any biomolecular target the majority of the chemicals will be negative. The relevant problem is to identify the few chemicals that are positive. This approach is based on modeling the forces that determine the interaction between a small molecule and a biomolecular target. A score (a surrogate for that interaction energy) is obtained for each molecule. The actual determination of likely positive and negative chemicals depends on the choice of a line of demarcation (called Q-R in Figure 1) imposed on the continuous spectrum of scores for the model interaction. This

Molecular Modeling: Screening for Estrogenicity

value may be adjusted to provide a more or less conservative evaluation. In order to use this approach, there must be high quality structural descriptions of the biomolecular targets or similar macromolecules. Previous data on the activity of other chemicals acting at this target are not required. In this approach, a single step in a complex path for toxicity is modeled. It is an indicator of suitability of the potential toxicant to take part in that single facet of more general mechanisms for toxicity. There are undoubtedly other pathways and perhaps other steps in a single pathway for toxicity that may be influenced by a potential chemical toxicant. As a library of targets is developed, interactions with many targets may be combined to provide a more complete picture of potential chemical toxicity. The results generated by the toxicant-target approach may be integrated with other types of data, chemical information, pharmacophores, and bioassay data. These different types of data may be best combined in a tiered approach that recognizes the strengths of each type of data and the requirements for obtaining it. Acknowledgment. M.-R.G. was supported during a portion of this work by NHEERL-DESE Training Agreement, EPA CT829471. We thank Drs. Chang, Dix and Kavlock for reading and commenting on this manuscript and Dr. Richard for helpful advice. Supporting Information Available: Table 1 contains the name, CAS number, SMILES code, measured IC50, and the activity calls based on the criteria described in this study for each chemical. It also contains all the docking data (scores) used in this study, four targets ERR, ERβ, from both agonist and antagonist crystal structures and four methods, FRED without constraints, FRED with 2 constraints, eHiTS without constraints, and eHiTS results with the constraints from FRED superimposed on the scores. Table 2 contains the name, CAS number, SMILES code, and RBAs for ERR from ref 41. It also contains the docking data (scores) for both ERR, the agonist, and antagonist targets using both FRED and eHiTS without constraints.This material is available free of charge via the Internet at http://pubs.acs.org. This work was reviewed by the EPA and approved for publication but does not necessarily reflect official Agency policy.

References (1) Judson, R., Richard, A., Dix, D., Houck, K., Martin, M., Kavlock, R., Dellarco, V., Henry, T., Holderman, T., Sayre, P., Tan, S., Carpenter, T., and Smith, E. (2009) The toxicity data landscape for environmental chemicals. EnViron. Health Perspect. 117 (5), 685– 695. (2) Richard, A. M. (2006) Future of toxicology-- predictive toxicology: An expanded view of “chemical toxicity”. Chem. Res. Toxicol. 19 (10), 1257–1262. (3) National Research Council (NRC). (2007) Toxicity Testing in the 21st Century: A Vision and a Strategy. National Academy Press, Washington, DC. http://www.nas.edu. (4) Dix, D. J., Houck, K. A., Martin, M. T., Richard, A. M., Setzer, R. W., and Kavlock, R. J. (2007) The ToxCast program for prioritizing toxicity testing of environmental chemicals. Toxicol. Sci. 95, 9–12. (5) Rabinowitz, J. R., Goldsmith, M.-R., Little, S. B., and Pasquinelli, M. A. (2008) Computational molecular modeling for evaluating the toxicity of environmental chemicals: Prioritizing bioassay requirements. EnViron. Health Perspect. 116, 573–577. (6) Kuntz, I. D. (1992) Structure-based strategies for drug design and discovery. Science 257 (5073), 1078–1082. (7) Halperin, I., Ma, B., Wolfson, H., and Nussinov, R. (2002) Principles of docking: an overview of search algorithms and a guide to scoring functions. Proteins 47 (4), 409–443. (8) Sousa, S. F., Fernandes, P. A., and Ramos, M. J. (2006) Protein-ligand docking: current status and future challenges. Proteins: Struct., Funct., Bioinf. 65, 15–26.

Chem. Res. Toxicol., Vol. 22, No. 9, 2009 1601 (9) Abagyan, R., and Totrov, M. (2001) High-throughput docking for lead generation. Curr. Opin. Chem. Biol. 5 (4), 375–382. (10) Bissantz, C., Folkers, G., and Rognan, D. (2000) Protein-based virtual screening of chemical databases. 1. Evaluation of different docking/ scoring combinations. J. Med. Chem. 43, 4759–4767. (11) Cozzini, P., and Dottorini, T. (2004) Is it possible docking and scoring new ligands with few experimental data? Preliminary results on estrogen receptor as a case study. Eur. J. Med. Chem. 31, 601–609. (12) Amadasi, A., Mozzarelli, A., Meda, C., Maggi, A., and Cozzini, P. (2009) Identification of xenoestrogens in food additives by an integrated in silico and in vitro approach. Chem. Res. Toxicol. 22 (1), 52–63. (13) Ekins, S. (2004) Predicting undesirable drug interactions with promiscuous proteins in silico. Drug DiscoVery Today 9 (6), 276– 285. (14) Laws, S. C., Yavanhxay, S., Cooper, R. L., and Eldridge, J. C. (2006) Nature of the binding interaction for 50 structurally diverse chemicals with rat estrogen receptors. Toxicol. Sci. 94, 46–56. (15) Richard, A. M., and Williams, C. R. (2002) Distributed structuresearchable toxicity (DSSTox) public database network: A proposal. Mutat. Res. 499, 27–52. (16) Laws, S., Kariya, J., Wolf, M., and Richard A. M. (2009) DSSTox EPA Estrogen Receptor Ki Binding Study (Laws et al.) Database (KIERBL): SDF file and documentation, Launch version: KIERBL_v1a_ 278_17Feb2009, www.epa.gov/ncct/dsstox/sdf_kierbl.html. (17) http://epa.gov/endo/pubs/edspoverview/finalrpt.htm (accessed Mar 24, 2009). (18) NIH NTP/ICCVAM Report (2002) Current Status of Test Methods for Detecting Endocrine Disruptors: In Vitro Estrogen Receptor Binding Assays, NIEHS, RTP, NC. http://iccvam.niehs.nih.gov/docs/ endo_docs/final1002/erbndbrd/ERBd034504.pdf (19) Berman, H. M., Westbrook, J., Feng, Z., Gilliland, G., Bhat, T. N., Weissig, G. H., Shindyalov, I. N., and Bourne, P. E. (2000) The Protein Data Bank. Nucleic Acids Res. 28, 235–242. (20) Wa¨rnmark, A., Treuter, E., Gustafsson, J. A., Hubbard, R. E., Brzozowski, A. M., and Pike, A. C. (2002) PDB ID: 1GWR: Interaction of transcriptional intermediary factor 2 nuclear receptor box peptides with the coactivator binding site of estrogen receptor alpha. J. Biol. Chem. 277 (24), 21862–21868. (21) Shiau, A. K., Barstad, D., Loria, P. M., Cheng, L., Kushner, P. J., Agard, D. A., and Greene, G.L.,. (1998) PDB ID: 3ERT: The structural basis of estrogen receptor/coactivator recognition and the antagonism of this interaction by tamoxifen. Cell 95 (7), 927–937. (22) MOE, Chemical Computing Group, Inc., Montreal, Quebec, Canada, 2008. (23) Levitt, M. (1992) Accurate modelling of protein conformation by automatic segment matching. J. Mol. Biol. 226, 507–533. (24) Fechteler, T., Schomburg, U., and Dengler, D. (1995) Prediction of protein three-dimensional structures in insertion and deletion regions: A procedure for searching data bases of representative protein fragments using geometric criteria. J. Mol. Biol. 253, 114–131. (25) Cornell, W.D., Cieplak, P., Bayly, C. I., Gould, I. R., Merz, Jr, K. M., Ferguson, D. M., Spellmeyer, D. C., Fox, T., Caldwell, J. W., and Kollman, P. A. (1995) A second generation force field for the simulation of proteins, nucleic acids, and organic molecules. J. Am. Chem. Soc. 117, 5179–5197. (26) Pike, A. C., Brzozowski, A. M., Walton, J., Hubbard, R. E., Thorsell, A. G., Li, Y. L., Gustafsson, J. A., and Carlquist, M. (2001) PDB ID: 1HJ1: Structural insights into the mode of action of a pure antiestrogen. Structure 9 (2), 145–153. (27) Pike, A. C. W., Brzozowski, A. M., Hubbard, R. E., Walton, J., Bonn, T., Thorsell, A.-G., Engstrom, O., Ljunggren, J., Gustaffson, J.-A., and Carlquist, M. PDB ID: 2J7X: Structure of agonist-bound estrogen receptor beta LBD in complex with Lxxll motif from Ncoa5, unpublished work. (28) Halgren, T. (1996) Merck molecular force field. I. Basis, form, scope, parameterization and performance of MMFF94;. J. Comput. Chem. 17, 490–519. (29) http://www.simbiosys.ca/ehits/index.html (accessed Feb 22, 2009). (30) Zsoldos, Z., Reid, D., Simon, A., Sadjad, B. S., and Johnson, A. P. (2006) eHiTS: An innovative approach to the docking and scoring function problems. Curr. Protein Pept. Sci. 7, 421–435. (31) http://www.eyesopen.com/products/applications/FRED.html (accessed Feb 22, 2009). (32) McGann, M. R., Almond, H. R., Nicholls, A., Grant, J. A., and Brown, F. K. (2003) Gaussian docking functions. Biopolymers 68, 76–90. (33) Schultz-Gasch, T., and Stahl, M. (2003) Binding characteristics in structure-based virtual screening:evaluation of current docking tools. J Mol Model 9, 47–57. (34) Kellenberger, E., Rodrigo, J., Muller, P., and Rognan, D. (2004) Comparative evaluation of eight docking tools for docking and virtual screening accuracy. Proteins: Struct., Funct., Bioinf. 57, 225–242.

1602

Chem. Res. Toxicol., Vol. 22, No. 9, 2009

(35) Shughrue, P. J., Lane, M. V., Scrimo, P. J., and Merchenthaler, I. (1998) Comparative distribution of estrogen receptor-R (ER-R) and β (ER-β) mRNA in the rat pituitary, gonad, and reproductive tract. Steroids. 63 (10), 498–504. (36) Anstead, G. M., Carlson, K. E., and Katzenellenbogen, J. A. (1997) The estradiol pharmacophore: Ligand structure-estrogen receptor binding affinity relationships and a model for the receptor binding site. Steroids 62 (3), 268–303. (37) Stahl, M., and Rarey, M. (2001) Detailed analysis of scoring functions for virtual screening. J. Med. Chem. 44 (7), 1035–1042. (38) Jain, A. N. (2006) Scoring functions for protein-ligand docking. Curr. Protein Pept. Sci. 7 (5), 407–420. (39) Muthasa, D., Sabnis, Y. A., Lundborga, M., and Karle´n, A. (2008) Is it possible to increase hit rates in structure-based virtual screening by pharmacophore filtering? An investigation of the advantages and pitfalls of post-filtering. J. Mol. Graphics Modell. 26 (8), 1237–1251. (40) Peach, M. L., and Nicklaus, M. C. (2009) Combining docking with pharmacophore filtering for improved virtual screening. J. Cheminf. 1 (6), 1-13. (41) Kuiper, G. G. J. M., Carlsson, B., Grandien, K., Enmark, E., Hagblad, J., Nilsson, S., and Gustafsson, J.-A. (1997) Comparison of ligand binding specificity and transcript tissue distribution of estrogen receptors R and β. Endocrinology 138 (3), 863–879. (42) Carlson, H. A. (2002) Protein flexibility and drug design: how to hit a moving target. Curr. Opin. Chem. Biol. 6 (4), 447–452. (43) Sherman, W., Day, T., Jacobson, M. P., Friesner, R. A., and Farid, R. (2006) Novel procedure for modeling ligand/receptor induced fit effects. J. Med. Chem. 49, 534–553. (44) Mukherjee, S., Nagar, S., Mullick, S., Mukherjee, A., and Saha, A. (2008) Pharmacophore mapping of arylbenzothiophene derivatives for MCF cell inhibition using classical and 3D space modeling approaches. J. Mol. Graphics Modell. 26 (6), 884–992. (45) Yang, J.-M., and Shen, T.-W. (2005) A pharmacophore-based evolutionary approach for screening selective estrogen receptor modulators. Proteins: Struct., Funct., Bioinf. 59, 205–220.

Rabinowitz et al. (46) Ekins, S., and Erickson, J. A. (2002) A pharmacophore for human pregnane X receptor ligands. Drug Metab. Dispos. 30 (1), 96–99. (47) Podolyan, Y., and Karypis, G. (2009) Common pharmacophore identification using frequent clique detection algorithm. J. Chem. Inf. Model. 49 (1), 13–21. (48) Gaido, K. W., Leonard, L. S., Maness, S. C., Hall, J. M., McDonnell, D. P., Saville, B., and Safe, S. (1999) Differential Interaction of the methoxychlor metabolite 2,2-bis-(p-Hydroxyphenyl)-1,1,1-trichloroethane with estrogen receptors R and β. Endocrinology 140 (12), 5746– 5753. (49) Blair, R. M., Fang, H., Branham, W. S., Hass, B. S., Dial, S. L., Moland, C. L., Tong, W., Shi, L., Perkins, R., and Sheehan, D. M. (2000) The estrogen receptor relative binding affinities of 188 natural and xenochemicals: Structural diversity of ligands. Toxicol. Sci. 54, 138–153. (50) Laws, S. C., Carey, S. A., Ferrell, J. M., Bodman, G. J., and Cooper, R. L. (2000) Estrogenic activity of octylphenol, nonylphenol, bisphenol A and methoxychlor in rats. Toxicol. Sci. 54, 154–167. (51) Fang, H., Tong, W., Welsh, W. J., and Sheehan, D. M. (2003) QSAR models in receptor-mediated effects: The nuclear receptor superfamily. J. Mol. Struct: THEOCHEM 622, 113–125. (52) Liu, T., Lin, Y., Wen, X., Jorrisen, R. N., and Gilson, M. K. (2007) BindingDB: A Web-accessible database of experimentally determined protein-ligand binding affinities. Nucleic Acids Res. 35, D198–D201. (53) Ertl, P., Rohde, B., and Selzer, P. (2000) Fast calculation of molecular polar surface area as a sum of fragment-based contributions and its application to the prediction of drug transport properties J. Med. Chem. 43, 3714-3717 (as implemented in 22). (54) Wildman, S. A., and Crippen, G. M. (1999) Prediction of physiochemical parameters by atomic contributions J. Chem. Inf. Comput. Sci. 39 (5), 868-873 (as implemented in 22).

TX900135X