5702
J. Med. Chem. 2006, 49, 5702-5709
Virtual Screening for Aryl Hydrocarbon Receptor Binding Prediction Elena Lo Piparo,*,† Konrad Koehler,‡ Antonio Chana,† and Emilio Benfenati† Istituto di Ricerche Farmacologiche “Mario Negri”, Via Eritrea 62, 20157 Milano, Italy, and Karo Bio Computational & Medicinal Chemistry NoVum, SE-141 57 Huddinge, Sweden ReceiVed May 5, 2006
The overall goal of this study has been to validate computational models for predicting aryl hydrocarbon receptor (AhR) binding. Due to the unavailability of the AhR X-ray crystal structure we have decided to use QSARs models for the binding prediction virtual screening. We have built up CoMFA, Volsurf, and HQSAR models using as a training set 84 AhR ligands. Additionally, we have built a hybrid model combining two of the final selected models in order to give a single operational system. The results show that CoMFA, VolSurf, HQSAR, and the hybrid models gives good results (R2 equal to 0.91, 0.79, 0.85, and 0.82 and q2 0.62, 0.58, 0.62, and 0.70, respectively). Since the techniques analyzed show a good correlation and good prediction also for an external test set, particularly the HQSAR and the hybrid model, we can conclude that these models can be used for predicting AhR binding in virtual screening. Introduction Nuclear receptors are involved in the regulation of critical cellular processes such as regulation of cell growth, differentiation, and metabolic processes.1 They constitute an important super family of transcription regulators that include the dioxin/ aryl hydrocarbon receptor (AhR). Free AhR is located in the cytoplasm, associated with heat shock proteins. Ligand binding to AhR is presumed to produce conformational changes in the AhR protein, causing the translocation of the whole complex into the nucleus.2,3 Within the nucleus, the AhR-ligand complex dissociates from associated proteins and dimerizes with ARNT (its nuclear partner) to reconstitute an active transcription factor that binds specific DNA sequences.4 Other AhR ligands such as dibenzo-p-dioxin (TCDD) and coplanar polychlorinated biphenyls (PCBs) are potent toxicants widespread in the environment. Their resistance to metabolic breakdown along with their lipophility causes them to accumulate in the food chain, bringing about their relevant effects on human health.5,6 The enormous number of compounds within the human food supply makes it impracticable to screen all of them for nuclear receptor binding experimentally. However computational procedures are available that can rapidly assess the likelihood of a given compound to bind a given receptor.7-9 These rapid in silico methods can be used to prioritize compounds for followup experimental verification of nuclear receptor binding. Computational methods for affinity prediction may be broadly classified into two categories.10 When a detailed 3D structure of the protein receptor is known, then receptor fitting approaches can be done by docking a candidate ligand into the receptor cavity and using either molecular mechanics or an empirical scoring function to estimate the interaction energy, hence the affinity between the ligand and the receptor. The receptor structure can be obtained experimentally (e.g., X-ray crystallographic or NMR). Alternatively, if the structure has not been determined, but experimental structure(s) of closely related * To whom correspondence should be addressed. Present address: Bioinformatics Group, Department of Bioanalytical Science, Nestle Research Center, P.O. Box 44, CH-100000 Lausanne 26, Switzerland. Telephone: (0) 21 785 9530. Fax: (0) 21 785 9486. E-mail:
[email protected]. † Istituto di Ricerche Farmacologiche “Mario Negri”. ‡ Karo Bio Computational & Medicinal Chemistry Novum.
proteins are available, a homology model can be created by threading the sequence of the target receptor through an experimental template and mutating the corresponding amino acid residues in the template to match those of the target receptor. If no experimental structure or homology model is available, then a second method, called receptor mapping, can be used, which attempts to build a model of the receptor based on what binds to it. A variant of receptor mapping is 3D-QSAR (quantitative structure-activity relationship) in which a series of ligands with known affinity is aligned and then the strengths of the electrostatic and steric potentials of each ligand at regular grid points surrounding the molecule are correlated with the affinity of the compound. Once a QSAR regression equation has been generated for a training set of molecules, it can be used to predict the affinity of molecules not included in the training set. Another alternative to receptor mapping is to use other QSAR techniques that do not need structures optimization and superimposition. They are fast, user-friendly, and do not need human supervision, so they fit the purpose of virtual screening very well. To identify unknown endocrine disruptors in the food supply, the AhR is a critical receptor, since no crystallographic structure of this receptor is available upon which to base a virtual screen. For this reason we created a homology model of the ligand binding domain (LBD) of AhR based on the NMR structure of the C-terminal PAS domain of human HIF-2a (PDB code 1P97).11 However because of the relatively low sequence homology (∼25%) between the target AhR and the experimental template HIF-2a, virtual screening using this homology model is not likely to be accurate enough. As an alternative, we used QSAR models based on the experimentally determined binding affinities of dioxin and other families of AhR ligands to predict the activity of new ligands. Using published data and 3D-QSAR models by Waller and McKinney,12 we added further descriptors and we explored more sophisticated and automatic methods for ligands alignment (maximizing overlap of steric and electrostatic fields), which is the critical phase that determines the quality and the utilization of the resultant 3D-QSAR model for the prediction of a big library of compounds. Additionally, we explored alternative tools for virtual screening. One of the proposed within this work, VolSurf,13 does not require the alignment of the molecules; the other alternative, Hologram
10.1021/jm060526f CCC: $33.50 © 2006 American Chemical Society Published on Web 08/22/2006
Virtual Screening for Receptor Binding Prediction
Journal of Medicinal Chemistry, 2006, Vol. 49, No. 19 5703
Table 1. Dibenzo-p-dioxins
Table 2. Dibenzofurans
compda
R1
R2
R3
R4
R6
R7
R8
R9
pIC50
log P
SE
compd
R1
R2
R3
R4
R6
R7
R8
R9
pIC50
log P
SE
1 2 3 4 5 6* 7* 8 9 10 11 12 13 14 15* 16* 17* 18* 19 20 21 22* 23 24 25
H Cl H H Cl Cl Cl Cl H H Cl Cl Cl Cl H H H H Br Br Br Br H H H
Cl Cl Cl Cl Cl H Cl Cl Cl Cl Cl Cl Cl H Br Br Br Br H H Br Br Br Br Br
Cl Cl Cl Cl Cl Cl H Cl Cl H Cl H Cl H Br Br Cl Cl Br Br H Br Br H H
H H H H Cl H Cl Cl H H Cl Cl Cl H H H H H H H Br H H H H
H H Cl Cl H H H H H H H H Cl H H H H H H H H H H H H
Cl Cl Cl H Cl Cl Cl H Cl H Cl H Cl H Br Cl Cl Cl Br Br Br Br Br Br H
Cl Cl H H Cl Cl Cl H H Cl H H Cl H Br Cl Br Cl H Br Br Br H H H
H H H H H H H H H H H H Cl H H H H H Br H H H H H H
9.144 8.118 7.768 7.61 7.49 6.975 6.811 6.728 8.171 6.281 5.937 5.585 5.715 4.572 10.086 10.093 10.687 9.074 8.038 9.943 8.881 9.35 10.209 8.927 7.464
6.35 6.84 6.22 5.74 7.32 6.43 6.92 5.84 5.74 5.12 6.71 5.43 8.30 4.10 6.99 6.67 6.67 6.51 7.15 7.07 7.72 7.64 6.22 5.45 4.41
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63
H H H H H H Cl Cl H H H H H Cl H Cl Cl Cl Cl Cl Cl H Cl Cl Cl H H Cl Cl Cl H Cl H H Cl Cl H Cl
Cl H H Cl Cl Cl H H Cl Cl Cl Cl Cl H Cl Cl Cl Cl Cl Cl Cl Cl Cl Cl Cl Cl Cl Cl Cl H Cl Cl H Cl Cl H Cl Cl
H Cl H Cl H H Cl Cl Cl Cl H Cl Cl Cl Cl H H H Cl Cl H Cl Cl Cl H Cl Cl Cl Cl Cl Cl Cl H Cl Cl Cl Cl H
H H Cl H H H H H Cl H H Cl Cl H H Cl Cl Cl Cl H Cl Cl Cl H Cl Cl H H H Cl Cl H H Cl H Cl Cl Cl
H H H H Cl H Cl H H H Cl Cl H Cl H H Cl H H H H H H Cl Cl Cl Cl Cl H H H H H H H H H Cl
H H H H H H H H H H Cl H H H Cl H Cl Cl H Cl Cl Cl Cl Cl Cl Cl H H Cl Cl Cl Cl H Cl Cl Cl Cl H
H H H H H Cl H Cl H Cl H H Cl Cl Cl Cl H H Cl Cl Cl Cl Cl Cl Cl Cl Cl H H Cl H H H H H Cl H Cl
H H H H H H H H H H H H H H H H H Cl H H H H H H H H H H H H Cl Cl H H H H Cl H
4.061 5.003 3.429 6.088 4.125 4.103 6.123 4.653 5.396 6.858 7.255 7.379 7.657 7.61 8.444 5.715 8.194 5.371 7.911 8.147 6.728 8.943 7.587 7.508 5.808 8.376 7.61 7.379 7.954 7.657 7.657 7.313 3.429 8.689 7.954 7.623 7.623 6.297
4.48 4.48 4.48 5.04 5.12 5.12 5.82 5.82 5.59 5.67 5.67 6.23 6.23 6.53 6.23 6.42 6.98 7.13 6.83 6.79 6.98 6.79 7.39 7.34 7.53 7.34 6.38 6.23 6.23 6.98 6.93 6.93 3.84 6.23 6.23 6.98 6.93 7.13
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
a
The compounds of the test set are marked with an asterisk.
QSAR, is based on molecular fragments and also does not require structure optimization. Experimental Details Data Set. We used a set of 93 AhR ligands, splitted into a training and test sets, with experimental binding affinities as in the original article by Waller and McKinney,12 considering the negative logarithm of the chemical molecular concentration necessary to displace 50% of radiolabeled 2,3,7,8-tetrachlorodibenzo-p-dioxin (TCDD) from the Ah receptor and reported as pIC50. These data came from three different laboratories using 2,3,7,8-tetrachlorodibenzofuran (TCDF) as internal standard, to provide normalization for interlaboratory variability. All pIC50 values used to build up the models were normalized to a value of 8.444 for TCDF. From the original data set of 99 compounds it was decided to remove five of them, for which precise binding data were not available, and their removal greatly improved also the Waller et al. model.12 Two compounds in the data set used by Waller et al. have identical structures and affinities (dibenzofurans called 57 and 63 in the original article), so we used only one, eliminating the duplication. Our models were through and applied for the screening of a list of artificial chemicals with high exposure risk identified given by the European Commission within CASCADE (EU contract no. FOODCT-2004-506319). Within the list of compounds to screen only nine have known binding activity, and we recognized that they were already included in the data set selected. Therefore, we excluded them from the training set and they have been used as the external set to validate our virtual screening method. Thus, in this work we used 84 compounds as a training set, plus the nine compounds in the external test set. They include the dibenzo-p-dioxins, dibenzofurans, biphenyls, naphthalenes, indolocarbazoles, and indolocarbazoles derivatives listed in Tables 1-7. Molecular Modeling and Alignment Rules. The 3D atomic coordinates of the compounds were extracted from the SMILES using CORINA software.14 A rough geometry optimization was performed using Schrodinger premin and Schrodinger Bmin MMFF/ MCMM for a stochastic conformational search including full geometry optimization to find the global energy minimum.15
We also considered some constrains of the central torsion angle to 0° for nonplanar ligands (e.g., biphenyls) using MacroModel/ MMFF and Maestro 7.0 graphical user interfaces and setting the torsion constraint to 4.182 kJ/mol (1 kcal/mol). Charges used within this work were based on the density functional B3LYP ab initio calculations at the 6-311G+ level. The SEAL program that maximizes the overlap of steric and electrostatic fields16 was used for aligning the ligands. Aligning the ligands using SEAL with the electrostatic parameter reduced the default value of 1.0-0.33, increases the alignment quality of the smaller ligands such as biphenyls that achieve better results, since steric components are now more heavily weighted than electrostatic ones. Calculation of Descriptors (CoMFA, VolSurf, Hologram, log P, SE). The CoMFA analysis was done on a Linux-based PC workstation using the software package SYBYL version 7.1.17 The aligned molecules were placed in a three-dimensional grid space with the dimensions automatically set by the program and 1.5 Å (x, y, z) grid stepping. An absolute maximum of 30 kcal/mol for the steric and electrostatic energy calculated at each grid point was established experimentally. The CoMFA descriptors in terms of van der Waals (steric) and Coulombic (electrostatic) interactions were calculated using an sp3 carbon atom with a +1 charge as a probe. Equal weights were assigned to steric and electrostatic fields using the CoMFA standard scaling procedure implemented in SYBYL. The important issue of the reduction of the number of descriptors was considered.18 In the case of the CoMFA model,
5704
Journal of Medicinal Chemistry, 2006, Vol. 49, No. 19
Lo Piparo et al.
Table 3. Biphenyls
compd
R2
R3
R4
R5
R6
R2′
R3′
R4′
R5′
pIC50
log P
SE
64 65 66 67 68 69 70 71 72 73 74 75 76 77
H H H Cl Cl Cl Cl Cl Cl Cl Cl Cl Cl Cl
Cl Cl Cl H Cl H Cl Cl H Cl H H Cl H
Cl Cl Cl H Cl Cl Cl Cl Cl Cl Cl Cl Cl Cl
H Cl Cl H H Cl Cl Cl Cl Cl H Cl Cl H
H H H H H H H H H H H H H Cl
H H H H H H H H H H Cl Cl H H
Cl H Cl Cl Cl Cl H Cl Cl Cl H H H Cl
Cl Cl Cl Cl Cl Cl Cl Cl Cl Cl Cl Cl H Cl
H H H Cl H H H H Cl Cl H Cl H Cl
7.028 5.204 7.871 5.584 6.134 5.762 6.157 6.057 5.482 5.885 4.442 4.689 4.405 4.577
6.41 6.37 6.99 6.24 6.82 6.99 6.78 7.40 7.57 7.98 6.41 7.57 6.10 7.57
1 1 1 1 1 1 1 1 1 1 1 1 1 1
Table 4. Naphthalenes
Table 6. Indocarbazole Derivates
compd
R1
R2
R3
R4
R5
R6
R7
R8
pIC50
log P
SE
compd
X
pIC50
log P
SE
78 79 80 81 82
H H Br Br Br
Br Br Br Br Br
Br Br H Br Br
H H Br Br H
H H H H Br
H Br Br Br Br
H Br Br Br Br
H H H H H
5.616 7.668 7.465 7.608 7.996
4.67 6.19 7.10 7.66 7.62
0 0 0 0 0
90 91
C N
8.602 6.863
6.97 4.84
0 0
Table 7. Indocarbazole Derivates
Table 5. Indolocarbazoles
compd
X
R1
R2
R3
R4
pIC50
log P
SE
83 84 85 86 87 88 89
N S N N N N O
CH3 H H CH2CH3 COCH3 H H
H H H H H CH3 H
H H H H H H H
H H H H H H H
8.921 8.482 8.444 8.051 7.951 7.721 7.538
5.55 6.57 5.14 6.48 3.64 6.26 5.74
0 0 0 0 0 0 0
where thousands of descriptors are calculated during the analyses, the performance of different models considering the balance between grid stepping and molecular filtering was assayed. Finally, we selected a small grid stepping of 1.5 Å and an high column filtering of 2.5 to take into account most of the information and to avoid the problem of the computational time being intolerably long. This “column filtering” technique reduced the number of columns in the QSAR Molecular Spread Sheet to 435. We have calculated as well VolSurf descriptors, which do not need structures to be aligned. The program calculates energetically favorable interaction sites around the molecules to produce a 3-D molecular interaction fields (MIF)19 grid map that is compressed into a few quantitative 2-D numerical descriptors that are physi-
compda
R1
pIC50
log P
SE
92* 93*
H CH3
7.319 6.857
5.67 6.67
0 0
cochemically meaningful.13,20 We used five probes (water, hydrophobic, carbonyl oxygen, carboxy oxygen, and amphipathic) to characterize the interaction sites around target molecules. Threedimensional molecular field maps were transformed into 118 descriptors by VolSurf 4.0. Such descriptors were molecular volume (V), surface (S), molecular weight (MW), critical packing (CP), size of the hydrophilic (W) and hydrophobic (D) region, hydrogenbonding properties (HB), integy moments and hydrophobic integy moment, and local interaction energy minima, which represent the energy of the best three local minima of interaction energies between the water probe and the compound. Integy moments are vectors pointing from the center of the mass to the center of hydrophilic and hydrophobic regions, respectively. The HQSAR approach uses as molecular descriptors holograms that encode a fixed length array containing counts of a priori defined substructural fragments. This method uses only 2D structure information, thus avoiding the usual conformational flexibility and structure alignments problems. Holograms were generated using the standard parameters implemented in Sybyl 7.1. Molecular fragments were generated using the fragment size default (minimum 4, maximum 7) and the following fragment distinctions: atoms,
Virtual Screening for Receptor Binding Prediction
Journal of Medicinal Chemistry, 2006, Vol. 49, No. 19 5705
bonds, and connections. The HQSAR analysis was done by screening the 12 default series of hologram length values from 53 to 401 bins. The fragment patterns counts from the training set compounds were then related to the measured biological activity, and the best HQSAR model hologram length found was 257 bins. The logarithm of octanol-water partition coefficient, log P, was calculated using the Pallas 3.021 package. An indicator variable that takes into account strain energy (SE) was also added as descriptor for all ligands that were not torsionally constrained and in fact are, in the global energy minimum, set SE ) 0.0 and, for torsionally constrained ligands, set SE ) 1.0. Statistical Analysis. Statistical analysis was done using the partial least squares (PLS) method as employed in the QSAR module of SYBYL 7.1 and VolSurf 4 running on a Linux-based PC workstation. PLS is based on linear transformation of the descriptors’ space, producing a new variable space based on a small number of orthogonal factors (latent variables), so there is no correlation. This method is particularly useful when the number of variables equals or exceeds the number of compounds (data points), because it leads to stable, correct, and highly predictive models, even for correlated descriptors.22,23 To establish the maximal number of components contributing to the model with the lowest standard error of estimate (SEE), PLS analysis was combined with leaveone-out (LOO) cross-validation (CV). LOO-CV implies exclusion of each compound of the training set and the prediction of its activity by the model developed using the remaining compounds. To assess the goodness of the model we used the cross-validated coefficient q2, which expresses the model’s ability to reproduce the training set. It was calculated as follows: q2 ) 1 -
∑(Y ∑(Y
pred
- Yobsd)2
obsd
- Ymean)2
(1)
where Ypred, Yobsd, and Ymean are predicted, actual, and mean pIC50, and ∑(Ypred - Yobsd)2 is the predictive sum of squares known as PRESS. For each model, the LOO-CV predictions were examined. The models were subsequently validated also using cross-validation with 10 groups. In this way, the data set is randomly divided into 10 sets with approximately equal size and class distributions. The model is trained afterward using all but one of the 10 groups and then tested on the unseen group. This procedure is repeated for each of the 10 groups. The cross-validation score is the average performance across each of the 10 training runs. After establishing the optimal number of components, the PLS procedure was repeated without cross-validation while being given the exact number of components contributing to the final model as input. When additional columns (SE and logP descriptors) were added to the molecular spread sheet, the QSAR standard scaling method was set during the PLS analysis to weight the SE and log P columns as heavily as the CoMFA and Hologram descriptors. Since the high LOO q2 is the necessary condition but not a sufficient one for a model to have a high predictive power,24 we also used cross-validation with 10 groups and an external test set of compounds never seen by the models. For the ideal model, the slope and the correlation coefficient is equal to 1 and intercept is equal to 0. Since a good QSAR model may have a high predictive ability if it is close to the ideal one, we have set the intercept of the test set plot to 0. A hybrid model was also developed in order to improve the final results. To build combined models we used in-house software built as a PC-Windows Excel macro. We have selected the rule-based approach that consists of dividing the results interval into three main areas, where some of the values maximum, minimum, or mean from the selected models will rule the trend of the final model. In this way we have obtained a noncontinuous function that can be expressed as combinations of simple linear equations such as: pIC50 calc) kn[Min,Mean,Max(models to combine)] + an (2) The result can be regarded as a new set of rules whose final
expression is a hybrid system able to combine different models. The final expression found can be expressed in the following way. If mean(models to combine) > 6.642 pIC50 calc ) 1.049[Min(models to combine)] + 0.002 If mean(models to combine) > 3.604 pIC50 calc ) 1.005[Min(models to combine)] - 0.009 Otherwise pIC50 calc ) 1.021[Mean(models to combine)] - 0.106 A systematic variation of the combinations of maximum, minimum, and mean values, a decision on which models to combine, and a later optimization of the values of V1, V2, k1, a1, k2, a2, k3, a3 to give a better value for r2 were carried out. The optimization has been performed by means the downhill simplex method modified to search uphill for the higher r2 values. A cross-validation has been also performed to obtain the q2 value as described in eq 1 for the new hybrid model.
Results Molecular Modeling and Alignment Rules. The constraint of the central torsion angle to 0° for nonplanar ligands (biphenyls) had very little effect on the torsion angle (
