Introduction
A key step in predicting the toxic effects of polychlorinated biphenyls (PCBs), polychlorinated dibenzofurans (PCDFs), and polychlorinated dibenzo- p-dioxins (PCDDs) is the estimation of their binding to a common intracellular cytosolic protein called the aryl hydrocarbon receptor (AhR) (1, 2). This receptor controls the induction of hepatic cytochrome P4501A1 and associated aryl hydrocarbon hydroxylase (AHH) and 7-ethoxyresorufin O-deethylase (EROD) activities (3-5). Moreover, the relative affinity of individual PCBs, PCDFs, and PCDDs for the receptor has been correlated with many toxic responses such as thymic atrophy, body weight loss, immunotoxicity, and acute lethality (1-8).
Because of the importance of the AhR in determining toxicity, there has been a number of attempts to model the relationship between receptor binding and structure of xenobiotic chemicals. The simplest models are based on planar chemicals, such as PCDDs, and propose that planar ligands which fit into a 3 by 10 Å rectangle bind more tightly than nonplanar chemicals (1). Recent comparative molecular field analyses by Waller and McKinney (9, 10) suggest that AhR ligands could be approximately 12-14 Å in length and 5 Å in depth (z-axis). Based on these models (1, 9, 10), PCDFs are expected to behave similarly to PCDDs; however, except for planar PCB congeners, the affinity of most PCBs should be much lower than that of PCDDs.
Building on this modeling construct, there have been a number of attempts to develop reasonable quantitative structure activity relationship (QSAR) models for the binding affinity of chemicals to the AhR. Safe and co-workers (2, 8) developed linear free-energy relationships involving substituent constants and indicator variables for PCB affinity data (4, 8) and concluded that steric factors probably do not play a significant role in the interaction at the AhR binding site, except in the case of large substituents such as C6H5, n-C4H9, and t-C4H9. Moreover, these authors were the first to suggest that the absolute planarity of PCBs might not be required for effective binding. This hypothesis was supported by ab initio quantum mechanical calculations and crystallography studies that suggested polarizability and electron-acceptor properties of the ligands, rather than size and planarity (or coplanarity), can control the affinity of binding to the AhR (11-18).
More extensive efforts to develop QSARs for AhR binding have incorporated dispersion interactions between the receptor and aromatic xenobiotics (10, 19-21). In these models, the equilibrium binding constants for AhR interactions are related to the molecular polarizability and chemical-to-receptor separation distances. The importance of a stacking process in binding to the AhR has also been proposed to explain the binding of PCBs when the phenyl ring with the greatest degree of chlorination is assumed to be parallel with the receptor (approximated by a porphyrin moiety) and the plane of the other phenyl ring is rotated to a minimum energy conformation consistent with quantum-mechanical calculations. This model, developed by varying the overall distance from the receptor to optimize the correlation with experimental data, gives a good qualitative picture of the binding process.
Finally, it also has been hypothesized (9, 20, 21) that chemicals such as PCBs act as electron acceptors in charge-transfer interactions with the AhR, not only through the aromatic rings but also through chlorine or other substituents. In considering that the most polarizable chlorine atoms in PCBs are those in the lateral positions, a second charge transfer interaction, or cleft binding, between these lateral chlorine atoms and putative electron donor regions that likely constitute the interior pockets in the receptor has been proposed (7, 22).
The purpose of this study was to attempt to use the mechanistic concepts of AhR binding to develop a robust QSAR model applicable to multiple classes of chemicals. In conducting this analysis, we hypothized that QSARs developed with only planar structures (e.g., PCDDs) or only highly hydrophobic chemicals, cannot delineate all the important parameters controlling binding because the range of variation in chemical structure is too small. Thus, we used literature data of AhR binding for a variety of PCDDs, PCDFs, and PCBs with different substituent patterns. However, incorporation of flexible nonplanar structures, such as certain of the PCBs, requires rigorous assessment of stereoelectronic indices, which can vary more among the various conformers of the same molecule than among the conformers of different chemicals in a series under investigation (18, 25). Because there is little reason to expect that flexible molecules would assume the lowest energy (gas phase) conformation when binding to a macromolecule, we used an approach to exhaustively explore the importance of all energetically feasible conformations for each chemical (23). This technique permitted the examination of AhR binding affinities for planar chemicals as well as those that could assume nearly planar conformations.
Data. Previously reported data for the relative affinity for binding of PCBs, PCDDs, and PCDFs to the AhR are reproduced in Tables 1 and 2 (2, 4, 8, 24). The PCBs fall into relatively distinct groups; group A (congeners 1-14) consists of PCBs with varying chlorination patterns on both phenyl rings, while group B (congeners 7,13,15-30) is composed of a wide range of 2,3,4,5-tetrachlorobiphenyl derivatives in which the para-substitutent on the nonchlorinated ring was varied. The concentration of the test chemical necessary to reduce specific binding of [3H]TCDD (2,3,7,8-tetrachlorodibenzo- p-dioxin) to 50% of the maximal value in the absence of the competitor defines the EC50; log(1/EC50) values are indicated in Tables 1 and 2. (Table 2 also includes AhR binding data for several brominated dioxin derivatives; however, for the sake of simplicity, we refer to this data set collectively as PCDDs).
Generation of conformations. Using the OASIS system and the 3DGEN algorithm with conformer screening techniques, we derived sets of conformations that were then used to explore relationships between specific electronic properties and AhR binding (23, 25-27). The conformation generation technique was described in detail by Ivanov et al. (23) and consists of a combinatorial procedure to generate conformations in the context of certain steric constraints (e.g., distances between nonbonded atoms, ring-enclosure limits, torsional resolution) and expert rules (e.g., likelihood of intramolecular hydrogen bonds, cis/ trans or +/- isomers). The approach incorporates the conformational flexibility of saturated cyclic molecular fragments, as opposed to other conformational analysis techniques, which only explore conformational space formed by rotations around acyclic single bonds. As expected, for the relatively rigid PCDDs and PCDFs, all generated conformations approximated the minimum energy (planar) structures. However, for the PCBs of Group A, 102 conformations were generated using a torsional resolution of 30° and a nonbonded cutoff of 1.8 Å. For those of group B, 859 conformations were found.
Selection of conformations. Molecular descriptors of interest were computed either after geometry optimization or by direct single point (ISCF) calculations of the unoptimized conformations. For the purpose of screening both optimized and nonoptimized conformations, we introduce a new steric parameter for planarity (P1) because previous work had shown the importance of planarity in some binding interactions. P1 was defined as the normalized sum of the absolute values of torsional angles between all nonhydrogen atoms in the molecule:
where ntors is the total number of torsional angles, with |
tors|=|tors| for |tors|<90° and |tors|=180-|tors| for |tors|>90°. Although P1 is not intended to be used as an independent variable in developing QSARs, the planarity parameter permitted the screening of conformations with respect to the possibility of their assuming a planar conformation within the designated energy constraints discussed below.
To successfully derive robust QSARs, care must be taken to ensure that molecular conformations included in the models are energetically reasonable, i.e., there must be a systematic way to screen and exclude high energy structures. For the PCDFs and PCDDs, this was not an issue because conformers generated for these chemicals approximated the minimum energy planar structures. In the case of the geometry-optimized PCB conformers from group A, the range in formation enthalpy values (
hf) for all conformations included in the QSAR analysis was within 20 Kcal/mol of the hf value for the lowest energy conformer for a particular molecule. This cut-off was selected under the assumption that energy provided during binding to the receptor is sufficient to facilitate interconversion among lower energy conformations, as has been suggested for the binding of ligands to the estrogen receptor ( 28, 29). In the case of the optimized PCBs from group B and nonoptimized PCBs for both groups, consideration of all conformations generated resulted in a range of hf values significantly greater than 20 Kcal/mol. To address this, we systematically excluded the most unstable (highly energetic) 20-30% (corresponding, respectively, to inclusion of conformers populating the 80 and 70% E-levels) from the QSAR analysis. This typically resulted in a range of hf values of less than 8 Kcal/mol for the conformations ultimately included in the QSARs and a corresponding significant increase in model statistics.
It should be noted that in generating our QSAR models there was no weighting/scaling of the various conformers, other than exclusion of structures based upon (lack) of planarity and energetics.
Selection and calculation of molecular descriptors. MOPAC 6 ( 30), running on a VAX 4500 (Digital Equipment Co., Maynard, MA), was used to generate physicochemical and reactivity parameters for the various conformers. The PM3 Hamiltonian was used for electronic structure assessment. A broad set of geometric and physicochemical parameters, as well as electronic descriptors, could potentially have been calculated for the compounds under study; however, the actual parameters used were restricted to those hypothesized to be associated with AhR binding affinity, based upon mechanistic considerations from earlier studies. First, as discussed above, we used the planarity index to obtain conformers that were most planar. Following this, we systematically eliminated the most energetic structures from the analysis. The physicochemical descriptor log P (octanol-water partition coefficient) was calculated using the method of Kellog et al. ( 31). Further selection of regressors was based upon the work of McKinney and co-workers ( 10, 32), who showed that chemicals which have greater ability to accept electron density through charge-transfer interaction should bind to the AhR with greater affinity than those with lower electron acceptor properties. We hypothesized, therefore, that these stronger electron acceptors should have a lower energy unoccupied frontier orbital (ELUMO), lower energy for the occupied frontier orbital (EHOMO), and a lower energy difference in these frontier orbitals (EHOMO -LUMO), which can be related to molecular reactivity ( 33, 34). The charge-transfer interaction should increase in more planar conformers because aligning the planes of the aromatic rings increases electron acceptor properties and because the planar configuration permits closer proximity to the putative electron donor region of the receptor. For binding of more nonplanar congeners, charge-transfer and hydrophobic interactions of the lateral substitutents would lead to the importance of steric parameters related to entropy. Global steric indices, such as the geometric analogue of Wiener topological index (GW), its informational counterpart (GIW), and the maximum distance (Lmax), should reflect the steric influence of these substituents. Greater values of GW and Lmax and smaller values of GIW should also reflect an increase in molecular polarizability if more polarizable substituents were added to the structures. Additional descriptors examined for the QSAR analysis included acceptor superdelocalizability indices (SiE and SiN), frontier donor and acceptor superdelocalizability indices, and polarizabilities 
i (i denotes a specific atom in a molecule).
Results and Discussion
PCBs in Group A. As discussed above, attempts to model AhR affinity of PCBs using the lowest energy conformations have been unsuccessful. For example, the coefficients of determination ( r2) for the best models derived for energetically optimized PCBs from group A using ELUMO and EHOMO -LUMO as descriptors were 0.384 and 0.410, respectively. These results were not surprising because we also observed a large variation of many of the parameters among the different conformers of a single PCB congener. For example, the range of ELUMO for 2,3,3',4,4',5'-hexachlorobiphenyl was almost 0.35 eV over all conformers generated, whereas the range over all planar PCB conformers in group A was less than that amount (Table 1). In fact, the energy variation on frontier orbitals of PCB conformers for some congeners was greater than that among all planar congeners within the series.
|
|
Figure 1. Variation of observed aryl hydrocarbon receptor (AhR) binding affinity versus energy of the lowest unoccupied molecular orbital (ELUMO) for PCBs of group A: (A) lowest energy conformations; (B) optimized most planar conformations (see Eq. 2); and (C) nonoptimized most planar conformations (see Eq. 8). EC50, median effective concentration.
|
|
To illustrate this finding, Figure 1A presents the AhR binding affinity of PCBs as a function of ELUMO for the lowest energy conformers of the congeners in group A. The data reveal a poor correlation coefficient for two reasons. First, as we have demonstrated previously ( 35), biological activity is orthogonal to both hydrophobic and electronic descriptors except in special cases where the mechanism involved or the chemical series studied precludes the significance of one of the descriptors. Second, we suggest that the use of the lowest energy conformers, commonly chosen in computational packages, unnecessarily biases QSARs toward the inclusion of conformers that might be found in gas phase. There are many other plausible conformers for flexible structures at room temperature and in the presence of solvation forces and binding energy. For example, non- ortho PCBs have minimum energy structures that include torsional angles of approximately 42° derived from ab initio calculation. Yet, the rotational barriers for non- ortho PCBs are only 1-2 Kcal/mol, indicating almost free rotation at room temperature. Moreover, the planar conformers have ELUMO values that are significantly different than the lowest energy structure.
Figure 1B illustrates the improvement in the correlation when P1 was used to select the most planar conformers. For the ortho-substituted PCBs, the conformers were clustered in two regions of minimum energy. For example, the 10 conformers generated for 2,3,3',4,4'5'-hexachlorobiphenyl (#10) have two primary regions of minimal energy corresponding to torsional angles of 55.5°-57.5° and 91.5°-93.5° between the two phenyl rings.
The best monoparametric QSARs for optimized, most planar conformers of PCB group A were obtained after the selection of conformers with the lowest values of P1. As hypothesized, ELUMO and EHOMO -LUMO were found to be best descriptors for the most planar cluster of conformers (Fig. 1B) as follows:
where n is the number of conformers (the number of actual compounds from which the conformers were generated is given in parentheses), r2 is the coefficient of determination, s2 is the variance, F is Fisher's criterion, and rcv2 and scv2 are the corresponding statistics for the "leave one out" cross-validation analysis ( 36).
The negative coefficient of the above equations supports the concept that the most active conformers are those with the highest electron acceptor properties and reactivity (lowest values of ELUMO and EHOMO -LUMO, respectively). We also found a high correlation between P1 and the above descriptors (e.g., r2(EHOMO -LUMO/P1) = 0.83), which provides additional support for the concept that planarity increases charge transfer and stacking.
Biparametric models also were derived for PCBs from group A. In addition to the parameters describing charge-transfer interactions, steric indices such as GW, GIW, and Lmax; physicochemical parameters such as log P; and local electronic parameters such as acceptor superdelocalizability indices (SiN), charges (qi), and polarizability (i) of meta- and para-positions of the second phenyl ring were evaluated for statistical significance. Local indices that describe the contribution of the lateral substituents were incorporated into the model at a lower level of statistical significance (90%) and have smaller statistical weight, wi. These results are illustrated by the following equations:
The negative coefficient between GIW and the positive coefficients with log P and Lmax are consistent with the assumption that the cleft type of binding is influenced by interaction with the lateral ( meta- and para-) chlorine atoms ( 7, 9). Moreover, the positive coefficient with SiN indicates that PCBs with greater electron acceptor properties at the lateral positions in the second phenyl ring have relatively stronger AhR affinities.
The receptor binding of the PCBs from group A also was modeled using nonoptimized structures. The straightforward selection of the planar conformations from the initial set did not yield statistically significant monoparametric models: r2 = 0.240 and 0.130 for ELUMO and EHOMO -LUMO, respectively. The high energy of some of the conformers was the likely explanation for these results. Further modeling confirmed this hypothesis. When the most planar 3D-isomers were selected from among the energetically more stable conformations (i.e., populating the 70% E-level-the 30% most energetic were excluded), statistically significant correlations were derived with parameters assessing the capacity of PCBs for the stacking type of interaction (Fig. 1C).
These results are noteworthy in the sense that correlations using nonoptimized but energetically reasonable planar conformations were comparable to those obtained using computationally intensive minimum energy techniques. The elimination of additional higher energy conformations (e.g., populations at 60% and 50% of E-level) without geometry optimization did not significantly change the model results (data not shown).
PCBs in group B. Group B contains tetrachlorobiphenyl with an array of para substituents on the second ring having widely varying polarity and electron acceptor properties. The substituents include H, Cl, Br, I, OH, OCH3, NO2, COCH3, NHCOCH3, C6H5, CH(CH3)2, F, CF3, CH3, C2H5, I-C3H7, n-C4H9, and t-C4H9. The variance in hydrophobicity in group B is large enough that log P should be an important independent variable. Indeed, in the subset containing the most planar, optimized conformers, log P and local meta- and para-electron acceptor properties of the second biphenyl ring were significant descriptors. The correlations obtained for these models were higher than the r2 = 0.55 associated with the models reported by Safe and co-workers ( 2, 8). Models with log P and acceptor superdelocalizability at the 3' and 5' meta positions (S3'N and S5'N) resulted in an r2 = 0.63 and s2 = 0.128. The QSARs were significally improved upon elimination of the 30% most energetic conformers (i.e., inclusion of those populating the 70% E-level) from the subset of the most planar structures. This is summarized in equations 10 and 11, and illustrated in Figure 2.
Figure 2. Variation of observed versus predicted aryl hydrocarbon receptor (AhR) binding affinity derived from the optimized most planar PCBs in group B (see Equation 10). EC50, median effective concentration.
While equations 10 and 11 seem different from those derived for PCBs in group A in that ELUMO is not as important, this was not unexpected because the constant highly chlorinated ring for PCBs in group B suggests that the stacking interaction also would be constant. Hence, the shift in descriptors in these equations to the local electron acceptor properties for lateral positions supports the concept that, for molecules that cannot achieve total planarity, charge transfer at those positions influences binding through putative cleft type interactions with AhR.
The electronic and hydrophobic effects were found to be essentially orthogonal, with a covariance of 0.28 (log P/S3'N). Thus, although substituents such as I and C2H5 have low electron acceptor properties (average S3'N values of all conformers are 0.300 and 0.290, respectively), their high log P values (7.3 and 6.85, respectively) compensate, with the net result being that these PCBs have relatively large AhR binding values [log(1/EC50) = 5.82 and 5.46, respectively]. Conversely, strong electron acceptors, such as NO2 (S3'N = 0.325), have relatively low binding capacity due to their high polarizability and low hydrophobicity (log P = 2.49). Finally, the highest binding affinity observed was for the CF3 derivative [log(1/EC50) = 6.43] because this substituent increases the electron-acceptor properties of the meta positions (S3'N = 0.314) and, at the same time has a high hydrophobicity (log P = 6.86).
Results similar to Equation 11 were obtained with nonoptimized PCB conformations from group B. The best model, based on log P and S3'N or S5'N, was obtained after the selection of the most planar conformations (e.g., r2 = 0.62 and s2 = 0.264). The models improved further upon elimination of the highest energy conformations. For example, elimination of 30% of the highest energy conformations resulted in the following model statistics: r2 = 0.75 and s2 = 0.112 for 79 isomers representing 18 PCBs. Again, these results may seem surprising because they were derived without geometry optimization, clearly one of the most costly aspects of modeling. We do not advocate eliminating optimization from QSAR; rather, we show only that optimization should not give the pretense that the structure derived is anything more than one optimized structure.
Because the planarity of PCB conformers appeared to be important in AhR binding, we combined the subsets of the most planar conformers from series A (37) and B (167), selected at 10 quintals. The best QSARs derived either for optimized or nonoptimized structures were unsatisfactory (e.g., r2 = 0.55 for the best models obtained from optimized conformers). This was probably due to an insufficient screening based only upon planarity. In further analyses, the most planar isomers occupying 20% of the highest E-levels were eliminated; this significantly improved the modeling results. The variation in binding data was explained satisfactorily with biparametric equations that included descriptors for the stacking type of interaction (ELUMO, E HOMO-LUMO) in combination with log P. An example from the set of geometry-optimized conformers is shown in equation 12.
PCDFs and PCDDs. The AhR binding affinity of PCDFs and PCDDs was initially investigated separately for each class. Unlike the multitude of conformations found for the flexible PCBs, only two types of conformations--bent and planar--were generated by the 3DGEN system for PCDDs, and only planar conformations were generated for PCDFs.
Statistically robust monoparametric models were not derived for predicting the AhR binding of PCDFs or PCDDs for either optimized or nonoptimized structures. For example, the coefficient of determination was only 0.55 for the correlation with ELUMO and the binding of PCDDs to the AhR. Because the results obtained for optimized PCDDs were relatively poor when compared to nonoptimized structures, we decided to use only the nonoptimized structures for calculating stereoelectronic parameters from the 1SCF calculations to derive biparametric models. These modeling results for AhR binding of PCDFs are summarized by Equations 13-16 and illustrated in Figure 3:
Figure 3. Variation of observed versus predicted aryl hydrocarbon receptor (AhR) binding affinity for PCDFs (see Eq. 13) EC50, median effective concentration.
The negatively correlated GIW parameter combines best with -ELUMO and +log P, whereas +Lmax combines with +GW and ELUMO. The maximum covariance between regressors was r2
The best biparametric QSAR for all nonoptimized conformations of PCDDs resulted in an r2 = 0.71. The correlation was significantly improved after the selection of only the planar conformations. The results can be summarized by the following equations and in Figure 4.
Figure 4. Variation of observed versus predicted aryl hydrocarbon receptor (AhR) binding affinity for PCDDs (see Eq. 17). EC50, median effective concentration.
The results emphasize the importance of steric, hydrophobic, and electron acceptor properties in modeling the binding affinities.
Despite the importance of log P and electron acceptor parameters, the two parameters could not be included together in an acceptable model for this group of PCDDs. Equally surprising, we were not able to generate significant models for AhR binding of PCDDs using geometry optimized structures. This may be attributable to problems encountered with the quantum-mechanical geometry optimization, especially for bromine-substituted dibenzo- p-dioxins. It was difficult to reach a convergence for almost all of the compounds, even using diverse optimization procedures for the different compounds from the test set. In fact, the degree of PCDD bending, corresponding to the initially derived structures (generated via 3DGEN), was strongly reduced after the geometry optimization. For almost all structures, the final optimized structure was completely planar.
PCBs, PCDDs, and PCDFs. The final modeling effort involved the combined group of the PCBs, PCDDs, and PCDFs. From the group of PCBs, the most planar conformers were selected, after the elimination of 20% of the highest energy (nonoptimized) conformations. Octachlorodibenzo- p-dioxin (#13, Table 2) was removed from the correlation sample because through-field interactions between multiple neighboring Cl atoms are not accounted for by the quantum chemical method that was used.
When all of the chemicals are included in the correlation sample, the basic model indicated that electron acceptor and steric descriptors were most important. However, the regression coefficient for mono- and biparametric expressions fell below 0.70. Triparametric models resulted in slightly better prediction of the AhR binding, with the best model given in Equation 20 and presented graphically in Figure 5.
Figure 5. Variation of observed with the predicted aryl hydrocarbon receptor (AhR) binding affinity for the most planar conformations of PCBs, PCDFs, and PCDDs (see Eq. 20). EC50, median effective concentration.
From this model, it can be seen that parameters directly related to the stacking type of interaction, such as EHOMO -LUMO, along with parameters that ostensibly predict both stacking and cleft types of interactions, such as GIW and Lmax, are important descriptors of AhR binding for the combined classes of planar PCBs, PCDDs, and PCDFs. The negative correlations between EHOMO-LUMO as a measure of molecular electronegativity and log(1/EC50) are consistent with those observed for ELUMO and EHOMO in previous equations because all parameters show that increases in electron acceptor properties will increase binding affinity.
Summary and Conclusions
This work demonstrates that the relative affinity of PCBs, PCDFs, and PCDDs for the AhR can be predicted from molecular descriptors reflecting the electron acceptor capability of the frontier orbitals and lateral substituents, as well as the hydrophobicity and polarizability of the chemicals. The QSARs were consistent with mechanisms proposed in the literature in which the ligands bind to a planar surface that has strong electron donor properties. One of the possible reasons QSARs using these mechanistically based stereoelectronic descriptors generally have not been successful in the past is that most studies have assumed that the most appropriate conformers on which to focus are represented by spectroscopic data, or minimum energy structures. Our analysis suggests that minimum energy structures may not be those of the chemicals when bound to the planar receptor, and the stereoelectronic descriptors derived from structures would be inaccurate of the real potential for charge-transfer interactions.
The analysis has also shown that developing quantitative mechanistic models for the binding affinity of even modestly heterogeneous sets of chemicals is limited by the ability of current methods to quantify structure. The fact that optimized structures often give less statistically robust results in QSARs than nonoptimized structures suggests that considerable effort needs to be made in understanding the geometry optimization algorithms within large families of chemicals. The importance of solution and receptor effects on selecting conformations for deriving stereoelectronic descriptors is critical and also requires further study. Finally, the QSARs presented herein have shown subtle preferences among the steric and the electron acceptor parameters, most of which are sensitive to the chemicals selected for study. Thus, we also plan to extend this QSAR analysis to include AhR binding data generated for other chemicals.
