# Combining Molecular Docking and Pharmacophore Models Predicts Ligand Binding of Endocrine-Disrupting Chemicals to Nuclear Receptors

Asma Sellami, Matthieu Montes, Nathalie Lagarde

► To cite this version:

Asma Sellami, Matthieu Montes, Nathalie Lagarde. Combining Molecular Docking and Pharmacophore Models Predicts Ligand Binding of Endocrine-Disrupting Chemicals to Nuclear Receptors. Chemical Research in Toxicology, 2025, 38 (4), 10.1021/acs.chemrestox.4c00455. [hal-05034696](https://hal.science/hal-05034696v1)

HAL Id: hal-05034696

https://hal.science/hal-05034696v1

Submitted on 15 Apr 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
1 Combining molecular docking and pharmacophore

2 models predicts ligand binding of endocrine

3 disrupting chemicals to nuclear receptors

4 Asma Sellami¹, Matthieu Montes¹,² and Nathalie Lagarde¹*

5 1 Laboratoire GBCM, EA 7528, Conservatoire National des Arts et Métiers, 2 rue Conté,

6 75003 Paris, France

7 2 Institut Universitaire de France, Paris, France

8 Corresponding Author: * nathalie.lagarde@lecnam.net

9

10

11 KEYWORDS. Nuclear receptors, endocrine disrupting chemicals, docking, pharmacophores

12
13 'For TOC only'

<figure><img src="image_2.png" /></figure>

14

15
ABSTRACT. Nuclear Receptors form a family of proteins capable of accommodating a wide variety of small molecules in their ligand binding domain, ranging from therapeutic compounds to endocrine-disrupting chemicals. The rapid identification of these compounds, especially within the latter category, is of paramount importance. Using data extracted from the CompTox Dashboard, an Environmental Protection Agency initiative, we assessed the effectiveness of a combination of molecular docking and pharmacophore models in identifying ligands binding to six nuclear receptors: androgen receptor, estrogen receptor alpha, estrogen receptor beta, glucocorticoid receptor, peroxisome proliferator-activated receptor gamma and thyroid hormone receptor alpha. For each nuclear receptor, we selected a specifically designed and optimized *in silico* protocol that, in conjunction with experimental assays, can prioritize compounds for further evaluation to detect any potential toxicological concerns.

## Introduction

Nuclear Receptors (NR) are a family of transcription factors involved in the regulation of numerous biological functions¹. These proteins share a similar organization in domains including a DNA (deoxyribonucleic acid) Binding Domain (DBD) and a Ligand Binding Domain (LBD)². The binding of small molecules into a buried pocket of the LBD is associated with conformational changes that activate or inhibit the DBD binding to specific DNA sequences and modulate genes expression. Endogenous lipophilic NR ligands include steroids and retinoids, but exogeneous compounds can also target NR LBD. This is notably the case of numerous therapeutic compounds and marketed drugs³⁻⁵, but also environmental compounds suspected to interfere with the hormonal system and cause harmful effects on human health, known as the Endocrine Disrupting Chemicals (EDCs)⁵,⁶.
EDCs are defined by the World Health Organization⁷ as an exogenous substance or a mixture that alters the functions of the endocrine system and consequently causes adverse effects in an intact organism, or its progeny, or (sub)populations. Being able to rapidly identify EDCs and thoroughly investigate their potential toxicity is of major importance, especially for the cosmetic, the pharmaceutical and the food industries. This allows for implementing risk management measures and therefore guarantee the protection of humans and of the environment from possible harmful exposure. Unfortunately, there is no standard experimental procedure that can predict with confidence whether a compound is a potential EDC⁸ because of their structural diversity and the variable, unclear mechanisms of action⁸. Various initiatives have been conducted in order to gather and share experimental data on EDCs and their mechanisms of action using high throughput assays ⁹,¹⁰. In particular, *in vitro* binding assays can be used to study the so-called “direct” mechanism of action of the EDCs. According to this mechanism, compounds act through direct binding to the active site of hormones in NR, and consequently disrupt endogenous ligand binding and modulate the biology of the endocrine system ¹¹. However, time and cost limitations make it impractical to use biological assays to test every chemical present in ecosystems. In this context, *in silico* modeling approaches are an asset to overcome such limitations ¹² and can be used in complement to experimental assays to identify potential EDC and gain insights about their possible mechanisms of actions⁵.

In a previous study¹³, we demonstrated that using a combination of molecular docking and pharmacophore models could be used to flag ligands binding to the estrogen receptor alpha (ERα) for further investigation of their potential EDC potency. In the present study, we used a similar approach to provide *in silico* protocols for 5 additional NR associated with direct mechanism of action of EDC: androgen receptor (AR), estrogen receptor beta (ERβ), glucocorticoid receptor (GR), peroxisome proliferator-activated receptor gamma (PPARγ)
and thyroid hormone receptor alpha (TRα), but also a modified protocol for ERα. For each of
these NR, molecular docking and pharmacophore modelling protocols were evaluated,
optimized and combined to assess their ability to identify compounds able to bind to one or
several of the 6 studied NR.

69 **1. Material and methods**

70 **1.1 Compounds collection and preparation.**

71 Compounds and their associated biological data have been downloaded from the CompTox
72 Dashboard ¹⁰, an EPA (Environmental Protection Agency) initiative. This dashboard gathers
73 results of biological and cellular assays against several biological targets including NR ¹⁴ for
74 compounds suspected to be of risk for humankind and the environment. Using the data
75 binding assays outcomes available in the CompTox Dashboard, 6 “EPA” data sets were
76 generated for the 6 NR studied, namely AR, ERα, ERβ, GR, PPARγ and TRα. Each data set
77 includes compounds annotated as active in the binding assays that will be considered as
78 binding compounds, and compounds annotated as inactive in the binding assays that will be
79 considered as non-binding compounds. The agonism and antagonist assays were used to
80 annotate the pharmacological profile of binding compounds when available (Figure S1).

81 In addition to the EPA sets, 6 validation sets for AR, ERα, ERβ, GR, PPARγ and TRα were
82 prepared from the NR-DBIND ¹⁵. The corresponding data sets were directly downloaded
83 from the website (http://nr-dbind.drugdesign.fr/) in SMILES format. Compounds for which
84 no binding data were available were discarded and binding and non-binding compounds were
85 selected according to the binding profile annotated in the NR-DBIND. These data sets will be
86 referred to as the NR-DBIND sets.

87 Both EPA and NR-DBIND data sets were then submitted to the preparation protocol
88 previously described ¹³. All compounds were initially stored in the data sets using SMILES
89 format. The Standardizer application from the ChemAxon Chemical Structure Representation
Toolkit ¹⁶ was used to remove salts, fragments, duplicates, and molecules containing less than 5 atoms. The cxcalc plugin calculator was used to compute the protonated state at pH 7.4 using the majormicrospecies option. The i-Con conformer generation tool of LigandScout version 4.5 ¹⁷ with BEST settings provided up to 25 conformers for each compound. All compounds were finally converted in pdbqt format (for molecular docking computations) and ldb format (for pharmacophore modelling).

217 physicochemical descriptors were calculated with the RDkit python package¹⁸ for each data set. These descriptors were used in a principal component analysis (PCA) performed with scikit-learn ¹⁹ and for the applicability domain definition.

The number of principal components was selected to capture at least 80% of the total variance, resulting in 50 components. Following PCA transformation, K-means clustering was performed on the first 50 principal components. The optimal number of clusters (k = 2) was chosen based on the Elbow method and Silhouette analysis, ensuring a balance between intra-cluster cohesion and inter-cluster separation. Each compound was assigned to one of the two clusters, and the distribution of binding compounds and non-binding compounds across clusters was examined to determine whether clustering effectively differentiates these two classes. Statistical significance of the clustering with respect to binding activity was assessed using Chi-square test for all NRs except TR alpha for which the Fisher's exact test was used. Two separate Z-tests were performed to compare the proportions of binding compounds and non-binding compounds between the two clusters.

## 1.2 Receptor protein structures preparation

For each of the 6 studied NR, all available (1) human structures, (2) with no mutations nor residue deletions in the ligand binding domain and (3) referenced in the literature were retrieved from the Protein Data Bank (PDB). Protein structures were prepared using the DockPrep tool of the Chimera software ²⁰ to remove alternate side chains locations, correct
115

incomplete side chains and add hydrogens and charges. The protein structures were converted
to pdbqt format for molecular docking.

116

117

**1.3 Molecular molecular docking methods**

118

Molecular molecular docking aims at predicting the conformation and the orientation of a
ligand inside a protein. All molecular docking calculations were performed with the smina
software²¹. Smina is a fork of AutoDock Vina ²² that uses a genetic algorithm based
conformation sampling. Four built-in scoring functions are available within smina: vina ²²,
Vina RaDii Optimized (vinardo) ²³, dkoes_scoring ²¹, and ad4_scoring. All molecular
docking calculations were performed using default options for the maximum number of
binding modes and the exhaustiveness of the global search. The binding site and the
bounding box coordinates were determined based on the aligned crystal structures for each of
the 6 studied targets. The box parameters were chosen based on the co-crystallized ligand
position with a spacing of 1 Angstrom. For each of the 6 studied NR, the search space was
delimited by a cubic box defined by appropriate x, y, z grid center coordinates (see Table
3). The size of the search space was set to 20 Å x 20 Å x 20 Å.

119

120

121

122

123

124

125

126

127

128

129

130

Both single and ensemble structure molecular docking approaches were explored for each
NR and each scoring function enrichment analysis. In the single structure molecular docking
approach, the performance in enrichment was evaluated for each PDB structure using the area
under the ROC curve (AUC). In the ensemble molecular docking approach, for each NR and
each scoring function, the single structure molecular docking outcomes were post-processed
for each ligand to only keep the pose associated with the best molecular docking score among
all the PDB structures of the query. All ligands were then ranked according to the selected
scores and the corresponding AUC were computed. Only ensembles of 2 and 3 structures
were considered according to previous studies conducted on AR and ERα¹³,²⁴.

131

132

133

134

135

136

137

138
The predictiveness curve<sup>25,26</sup> was used to define a scoring threshold able to discriminate binding compounds from non-binding ones. The predictiveness curve enables associating a predictivity to each score, which is the probability for a compound with a given score to be a binding compound. For each predictivity, the associated sensitivity (Se) and specificity (Sp) can be computed and used to select the best threshold. Using this protocol, we highlighted potential scoring threshold for each NR. Additionally, the evolution of several metrics (described below in the metric performance section) was assessed in function of the scoring threshold selected for each database to make sure that the selected threshold combines the perfect equilibrium between all the computed metrics.

The final molecular docking model selected for each NR is defined by a (n ensemble of) protein structure(s), with optimized molecular docking parameters (including bounding box coordinates), a selected scoring function and an identified scoring threshold. Compounds screened using these molecular docking models will be predicted as binding compound if their molecular docking scores are higher than the threshold, and as non-binding compounds otherwise.

## 1.4 Pharmacophore models

A pharmacophore is defined by the IUPAC as "the ensemble of steric and electronic features that is necessary to ensure the optimal supramolecular interactions with a specific biological target and to trigger (or block) its biological response"<sup>27</sup>. Pharmacophore models thus do not describe the location of specifically functional groups or structural fragments of a molecule but of geometric entities (spheres, vectors and planes) representing the molecule chemical features involved in protein-ligand interactions<sup>28</sup>. Both ligand-based (LB) and structure-based (SB) pharmacophore models can be generated. LB pharmacophore models are derived from the structural features of reference active compounds, without prior knowledge of or consideration for the protein structure. Conversely, when either an apo or
164

holo protein structure is used to construct the pharmacophore, the resulting pharmacophore model is referred to as SB pharmacophore. In this study, only holo protein structures were used to generate SB pharmacophore models.

165

166

167

The EPA databases and the selected PDB structures were used to LB and SB pharmacophore models respectively, using the software LigandScout version 4.5¹⁷.

168

169

For LB pharmacophores, each data set was split into a train and a test set using a 75%/25% ratio. Binding compounds of the train set were clustered based on their structural similarities using the “icluster” tool of LigandScout with the default parameters. Pharmacophore models were generated for each cluster and used to screen the corresponding train set. To avoid bias caused by the data separation, the splitting step was performed 25 times, after which individual and general performance of the pharmacophores were evaluated and the iteration associated with the best performance was selected. 3D SB pharmacophores were automatically generated from the PDB holo structures selected for each NR.

170

171

172

173

174

175

176

177

Separate groups of LB and SB pharmacophores models were generated for each NR, these models were then optimized individually using a protocol defined according to several literature recommendations ²⁹⁻³¹ and presented in Figure 1. For each pharmacophore, with an initial number of pharmacophoric features T, a first screening was made with the number of omitted features (OMF) set to 0, i.e., a hit is a compound matching all the pharmacophore features. The pharmacophores associated with an empty hitlist after this screening were discarded at this step. Next, for each of the remaining pharmacophores, a second round of screening was performed with the OMF set to T-3 in order to generate all possible combinations of pharmacophoric features. Each pharmacophore corresponding to one of these feature combinations was generated and used to screen the training set with the max number of omitted features set to 0. If the associated hitlist had a binding compounds / non-binding compounds ratio equal or superior to 1, the pharmacophore model was kept,

178

179

180

181

182

183

184

185

186

187

188
otherwise it was discarded. Finally, redundant pharmacophores, i.e., pharmacophores retrieving compounds already identified as hits by other pharmacophores from the same collection, were discarded. The same optimization protocol was applied for both SB and LB pharmacophores except that for the LB, an additional preliminary step was made to set as “obligatory”, all the “optional” features of the initial pharmacophores. Optimized SB and LB pharmacophores were combined to form a unique SBLB pharmacophores group after removal of redundant pharmacophore models between both sets.

<figure><img src="image_3.png" /></figure>

196

<caption>197 **Figure 1.** Optimization protocol for pharmacophore model generation</caption>

198 **1.5 Combination protocols and models' selection**

The combination of both molecular docking and pharmacophore models was investigated.

199
200 The consensus and hierarchical protocols previously described ¹³ were used (Figure 2). In the
201 consensus protocol, a compound is predicted as a binding compound if it is predicted so by at
202 least one of the models. Conversely, using the hierarchical protocol, only the compounds
203 identified as binding compounds using both pharmacophore and molecular docking models
204 are tagged as binding compounds for the final prediction.
<figure><img src="image_4.png" /><figcaption>**Figure 2:** Illustration of consensus (left panel) and hierarchical (right panel) protocols. Docking and pharmacophore models prediction for a compound to be a binding compound (in green) or a non-binding compound (in red) are combined in a parallel (consensus) or sequential (hierarchical) way to obtain a final prediction (in gray).</figcaption></figure>

## 1.6 Applicability domain

The applicability domain (AD) is defined as "*the response and chemical structure space in which the model makes predictions with a given reliability*"³². In this study, the AD of our final prediction models was computed for each NR using a two-step approach adapted from the approach proposed by Trisciuzzi et al.³³ and using the 217 molecular descriptors computed with the RDkit python package¹⁸. In the first step, for each EPA data set, the bounding-box³⁴ for each descriptor was defined with its minimum and maximum values across the EPA data set. In the second step, a Mahalanobis distance-based approach³⁵ was applied in a reduced-dimensional space obtained via Principal Component Analysis (PCA). The dimensionality of the descriptor space was reduced to 50 principal components, capturing at least 80% of the total variance in the EPA data sets. For each NR-DBIND data set, the Mahalanobis distance of each test compound was then computed relative to the distribution of the corresponding EPA compounds in this PCA-transformed space. A chi-squared distribution threshold at the 95th percentile with 50 degrees of freedom was defined.
224

Finally, a test compound was considered inside the applicability domain if (i) all its descriptor values were within the bounding box and (ii) its computed Mahalanobis distance was below the chi-squared threshold at the 95th percentile.

225

226

227

## 1.7 Comparison with CoMPARA and CERAPP prediction models

228

229

230

231

232

233

The OPEn structure-activity/property Relationship App (OPERA) is a free and open source tool that gathers QSAR models dedicated to toxicological purposes and ADME properties prediction. OPERA was used to evaluate the CoMPARA (Collaborative Modeling Project for Androgen Receptor Activity) QSAR binding models performance on the AR NR-DBIND dataset and the CERAPP (Collaborative Estrogen Receptor Activity Prediction Project) QSAR binding models on the ERα and ERβ NR-DBIND dataset.

234

## 1.8 Metric performance

235

236

237

238

239

240

241

242

243

244

245

The performance of the models generated during this study was evaluated according to their ability to distinguish binding from non-binding compounds for the studied NR targets. During the evaluation process, compounds screened by the models were divided into 4 categories: true positive (TP, experimentally confirmed binding compound correctly predicted as a binding compound), true negative (TN, experimentally confirmed non-binding compound correctly predicted as a non-binding compound), false positive (FP, experimentally confirmed non-binding compound wrongly predicted as a binding compound) and false negative (FN, experimentally confirmed binding compound wrongly predicted as a non-binding compound). Different metrics were then computed (Table 1, Table S1).
<table><caption>**Table 1.** Metrics computed for performance evaluation</caption><thead><tr><th>Metric</th><th>Equation</th></tr></thead><tbody><tr><td>Sensitivity (or recall, Se)</td><td>Se = <sup>TP</sup>&frasl;<sub>TP + FN</sub></td></tr><tr><td>Specificity (Sp)</td><td>Sp = <sup>TN</sup>&frasl;<sub>TN + FP</sub></td></tr><tr><td>Matthews correlation coefficient (MCC)</td><td>MCC = <sup>TP * TN - FP * FN</sup>&frasl;<sub>&radic;(TP + FP) * (TP + FN) * (TN + FP) * (TN + FN)</sub></td></tr><tr><td>Positive Likelihood Ratio (+LR)</td><td>+LR = <sup>Se</sup>&frasl;<sub>1 - Sp</sub></td></tr><tr><td>Negative Likelihood Ratio (-LR)</td><td>-LR = <sup>(1 - Se)</sup>&frasl;<sub>Sp</sub></td></tr></tbody></table>

Positive and Negative Likelihood Ratios were computed as these metrics remains stable across different class distributions, and are especially well-suited for evaluating models in unbalanced datasets.

The Se and Sp were used to plot the ROC curve and computed the associated AUC with python 3.8 using the scikitlearn library¹⁹.

## 1.9 ScreeNR script

To facilitate the use of the final prediction models developed in this study, we provide a Python script, ScreeNR.py (Zenodo repository: https://doi.org/10.5281/zenodo.14974415).
This script enables the screening of small molecules, starting from a 3D SDF format, against one or multiple NR prediction models presented in this study. The script applies the pre-trained prediction models to classify compounds as either binding or non-binding compound for each selected NR. Additionally, the script evaluates whether the screened compounds fall within the applicability domain (AD) of the model using the two-step approach presented above (see paragraph 1.6). The output of the ScreeNR.py script is a CSV file containing three columns: the compounds ID, the predicted binding status (1 for predicted binding compounds and 0 for non-binding compounds) and the applicability domain status (True if the compound is within the AD, False otherwise), respectively.

## 2. Results

### 2.1 Compounds collection analysis

The CompTox Dashboard was filtered to extract only the compounds tested in binding assays for the 6 studied NR and 6 data sets were constituted for AR, ERα, ERβ, GR, PPARγ and TRα (see Table 2). The largest data set is the one gathering compounds experimentally tested for ERα binding with 2242 compounds. Conversely, the smallest data set is the one dedicated to TRα that includes only 222 compounds. The number of non-binding compounds exceeds the number of binding compounds by at least a factor of 3 for all data sets, except for GR, for which an inverted trend is observed.

<table><thead><tr><th rowspan="2">NR</th><th colspan="3">EPA data sets</th><th colspan="3">NR-DBIND data sets</th><th colspan="3">NR-DBIND data sets within the AD</th></tr><tr><th>B</th><th>NB</th><th>B/NB</th><th>B</th><th>NB</th><th>B/NB</th><th>B</th><th>NB</th><th>B/NB</th></tr></thead><tbody><tr><td>AR</td><td>371</td><td>1401</td><td>0.26</td><td>583</td><td>150</td><td>3.89</td><td>473</td><td>101</td><td>4.68</td></tr><tr><td>ERα</td><td>223</td><td>2219</td><td>0.10</td><td>534</td><td>164</td><td>3.25</td><td>404</td><td>119</td><td>3.39</td></tr><tr><td>ERβ</td><td>242</td><td>1495</td><td>0.16</td><td>676</td><td>81</td><td>8.35</td><td>469</td><td>51</td><td>9.20</td></tr><tr><td>GR</td><td>272</td><td>126</td><td>2.16</td><td>1200</td><td>63</td><td>19.05</td><td>859</td><td>56</td><td>15.34</td></tr></tbody></table>
<table><tr><td>PPARγ</td><td>108</td><td>1653</td><td>0.07</td><td>423</td><td>168</td><td>2.52</td><td>218</td><td>82</td><td>2.66</td></tr><tr><td>TRα</td><td>54</td><td>168</td><td>0.32</td><td>89</td><td>14</td><td>6.36</td><td>74</td><td>4</td><td>18.5</td></tr></table>

**Table 2:** Composition of the EPA data sets, NR-DBIND data sets and compounds of the NR-DBIND data sets falling in the computed AD in number of binding (B) and non-binding (NB) compounds and its ratio (B/NB).

For each EPA data set, a PCA was computed to compare the distribution of 217 physicochemical descriptors¹⁸ within binding and non-binding compounds and visualize their chemical spaces coverage (Figure 3).

To assess whether binding and non-binding compounds could be distinguished based on this analysis, we performed clustering on the first 50 principal components. The distribution of binding and non-binding compounds across the two resulting clusters was then evaluated (Table S2). This analysis revealed a statistically significant difference in the distribution of binding and non-binding compounds in the two clusters for AR, ERα, GR, and PPARγ, whereas no significant difference was observed for ERβ and TRα. Z-tests show that the proportion of binding compounds differ significantly between the two clusters only for GR and TRα, and that the proportion of non-binding compounds differ significantly between the two cluster for all NRs.
<figure><img src="image_5.png" /><figcaption>**Figure 3:** PC1 vs. PC2 plots obtained by PCA of physicochemical properties between binding compounds (in green) and non-binding compounds (in red) computed for each NR dataset separately.</figcaption></figure>

The distribution of binding compounds common between different NR was evaluated (Figure 4). The number of unique binding compounds for each NR is higher than the number of compounds able to bind to several NR for all NR except PPARγ and TRα. However, almost half of the binding compounds included in the EPA data sets are able to bind to at least 2 NR (307 compounds out of the 796 binding compounds). AR and ERα data sets present the largest number of unique binding compounds with 156 and 154 unique binding compounds, respectively. The highest number of binding compounds in common between different data sets is observed between the AR and GR data sets that share 60 binding
303

compounds. No common binding compound is present in the 6 data sets and only 7 binding compounds are shared by 5 data sets.

304

305

<figure><img src="image_6.png" /></figure>

306

<caption>**Figure 4:** UpSet plot of binding compounds in common between the 6 EPA data sets</caption>

307

308

The composition of the NR-DBIND validation sets (Table 2) shows that for all sets, the number of binding compounds largely exceeds the number of non-binding compounds.

309

310

## 2.2 Molecular docking

311

In order to select the best molecular docking protocol for each NR, we evaluated the results obtained for all structures of each NR with the 4 investigated scoring functions and the 2 molecular docking approaches, i.e., single structure and ensemble structures molecular docking (Table S3 and **Figure 5**). Using the single structure molecular docking approach, for each NR except TRα and each scoring function, the best AUC obtained is higher than 0.6. Conversely, TRα presents very low overall binding compounds prediction performance since the best AUC values obtained with the 4 scoring functions are lower than 0.5. Among the 6 NR, the highest overall performance is associated with the PPARγ structures for all 4 scoring

312

313

314

315

316

317

318
functions, with the highest AUC exceeding 0.7 and the mean AUC values for the 4 scoring functions being equal to 0.7. Focusing on the 4 scoring functions, the best AUC values were obtained using the dkoes_scoring scoring function for 4 out of 6 NR i.e. AR, ERα, GR and TRα. For ERβ and PPARγ, the highest AUC values were obtained with the vina and the vinardo scoring function respectively. Using the ensemble structure molecular docking approach, similarly to the performance obtained with the single structure molecular docking approach, for each NR except TRα and each scoring function, the best AUC obtained was higher than 0.6. Mean AUC values associated with the ensemble molecular docking approach of 3 structures are slightly higher than the mean AUC values associated with the ensemble molecular docking approach of 2 structures in all cases except ERα and ERβ using the dkoes_scoring scoring function and ERβ using the ad4_scoring scoring function. The highest AUC values obtained with ensembles of 2 and 3 structures approaches are equal in 6 cases. In most of the remaining cases, very similar highest AUC values are obtained for both, the ensemble of 2 structures and the ensemble of 3 structures presenting the best AUC values in respectively 7 and 11 cases. No major AUC value improvement with the ensemble molecular docking approach compared to the single structure molecular docking approach was noticed for any NR and any scoring function, not even for TRα for which the associated best AUC values remained lower or equal to 0.5.
AR

<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /></figure>

ERA

<figure><img src="image_11.png" /></figure>

<figure><img src="image_12.png" /></figure>

<figure><img src="image_13.png" /></figure>

<figure><img src="image_14.png" /></figure>

ERB

<figure><img src="image_15.png" /></figure>

<figure><img src="image_16.png" /></figure>

<figure><img src="image_17.png" /></figure>

<figure><img src="image_18.png" /></figure>
GR

<figure><img src="image_19.png" /></figure>

PPARG

<figure><img src="image_20.png" /></figure>

TRA

<figure><img src="image_21.png" /></figure>

340 **Figure 5**: Distribution of the AUC of the ROC curves associated with the molecular docking outputs obtained for each one of the 6

341 studied NR using smina and its associated embedded scoring functions vina, vinardo, dkoes_scoring and ad4_scoring with the single
342 structure molecular docking approach (colored in blue) and the ensemble structures molecular docking approach (colored in orange

343 for ensembles of 2 structures and in green for ensembles of 3 structures).
For each NR, an optimal molecular docking protocol (Table 3) was selected except for TRα for which all evaluated molecular docking protocols resulted in unsatisfactory prediction performance (AUC values did not exceed 0.5). To complete these protocols, we used the predictiveness curve approach to rationally identify a molecular docking score threshold to distinguish predicted binding and non-binding compounds for the 5 remaining NR. The identified molecular docking threshold is equal to -7 for AR and PPARγ and -6 for the 3 other NR (Table 3). The evolution of several performance metrics as a function of the molecular docking scores was assessed (Figure S2) to confirm the threshold selection. The identified thresholds were incorporated into the selected molecular docking protocols (Table 3). Each protocol is characterized by a PDB protein structure, the bounding box coordinates, a selected scoring function and a defined scoring threshold (Table 3).

<table><thead><tr><th>NR</th><th>PDB ID</th><th>Profile of the ligand bound</th><th>Docking box center coordinates (x,y,z)</th><th>Scoring function</th><th>Docking score threshold</th></tr></thead><tbody><tr><td>AR</td><td>2pnu</td><td>Ago</td><td>(1.725, 32.527, 4.659)</td><td>dkoes_scoring</td><td>-7</td></tr><tr><td>ERα</td><td>1qku</td><td>Ago</td><td>(107.175, 14.938, 96.009)</td><td>dkoes_scoring</td><td>-6</td></tr><tr><td>ERβ</td><td>3ols</td><td>Ago</td><td>(33.383, 80.731, -9.056)</td><td>vina</td><td>-6.5</td></tr><tr><td>GR</td><td>4mdd</td><td>Atgo</td><td>(32.811, 7.088, 3.008)</td><td>dkoes_scoring</td><td>-6</td></tr><tr><td>PPARγ</td><td>2prg</td><td>Ago</td><td>(17.836, -21.534, 10.007)</td><td>vinardo</td><td>-7</td></tr></tbody></table>
356 **Table 3 : Selected molecular docking protocols for each NR**

357 **2.3 Pharmacophores**

358 Using our pharmacophore generation protocol, from 9 to 37 SB pharmacophores and from 34
359 to 253 LB were generated for each NR (Table 4). These pharmacophores were then all optimized
360 using the same protocol, described in the Material and Methods section. The number of
361 pharmacophore models was reduced along this optimization step with a number of optimized
362 pharmacophores ranging from 1 to 18 for SB pharmacophores and from 3 to 15 for LB
363 pharmacophores.

364
<table><caption>**Table 4:** Evolution of the number of pharmacophoric models throughout the optimization protocol. For each NR, redundant pharmacophores were identified among the LB and SB optimized pharmacophores and removed to create the SBLB pharmacophores set.</caption><thead><tr><th>NR</th><th>Pharmacophores</th><th>Generated<br>pharmacophores</th><th>Optimized<br>pharmacophores</th><th>SBLB<br>pharmacophores</th></tr></thead><tbody><tr><td rowspan="2">AR</td><td>LB</td><td>253</td><td>15</td><td rowspan="2">23 (13 LB / 10<br>SB)</td></tr><tr><td>SB</td><td>29</td><td>12</td></tr><tr><td rowspan="2">ERα</td><td>LB</td><td>140</td><td>8</td><td rowspan="2">21 (4 LB / 17 SB)</td></tr><tr><td>SB</td><td>15</td><td>18</td></tr><tr><td rowspan="2">ERβ</td><td>LB</td><td>120</td><td>4</td><td rowspan="2">14 (4 LB / 10 SB)</td></tr><tr><td>SB</td><td>24</td><td>10</td></tr><tr><td rowspan="2">GR</td><td>LB</td><td>232</td><td>8</td><td rowspan="2">9 (8 LB / 1 SB)</td></tr><tr><td>SB</td><td>9</td><td>2</td></tr><tr><td rowspan="2">PPARγ</td><td>LB</td><td>100</td><td>3</td><td rowspan="2">4 (3 LB / 1 SB)</td></tr><tr><td>SB</td><td>37</td><td>1</td></tr><tr><td rowspan="2">TRα</td><td>LB</td><td>34</td><td>3</td><td rowspan="2">3 (3 LB / 0 SB)</td></tr><tr><td>SB</td><td>16</td><td>1</td></tr></tbody></table>
In order to assess the impact of optimization, we compared the values of several metrics associated with the binding compounds prediction performance on the EPA training data sets of the generated pharmacophores on one side and the optimized pharmacophores on the other side (Figure S3 and Table S4). Focusing on TRα, no SB optimized pharmacophore was retained, and the LB optimized pharmacophores similar performance to the initially generated LB pharmacophores, except for MCC values, which were higher in the optimized pharmacophores. For the remaining 5 NR and for both pharmacophore categories (LB and SB), the optimization protocol generally resulted in pharmacophores with higher MCC values. However, the impact on the other performance metrics (Sp and Se) varied between SB and LB pharmacophores. For 4 out of the 5 NR (AR, ERα, ERβ and PPARγ), the optimized LB pharmacophores present higher Sp values but decreased Se values. Conversely, for the optimized SB pharmacophores of 3 NR (AR, ERα, and ERβ), the opposite trend was observed, with higher Se values but decreased Sp values.

For each NR, the LB and SB optimized pharmacophores were merged into an SBLB pharmacophore set with redundant pharmacophores between the LB and SB optimized sets removed. For TRα, since no SB optimized pharmacophores was retained in the SBLB pharmacophores ensembles after redundancy removal, the performance of the SBLB pharmacophore models is identical to that of the optimized LB pharmacophore models (Figure S3 and Table S4). A similar trend is observed for GR for which only one SB optimized pharmacophore was selected. Conversely, for ERα, ERβ, PPARγ, performance of the SBLB pharmacophore models is comparable with those obtained with SB pharmacophores (Figure S3 and Table S4). AR SBLB pharmacophore models performance shows more variations from the AR LB or SB pharmacophore models performance, with increased Se values but similar Sp and
390 MCC values. The SBLB pharmacophore models were selected for combination with the
391 previously identified molecular docking protocols.

392 **2.4 Combinations**

393 Two different protocols, called consensus and hierarchical, were evaluated to combine the
394 selected molecular docking and pharmacophore models for 5 out of the 6 NR. No combination
395 was possible for TRα since no molecular docking protocol had been previously selected. The
396 performance of the two combination protocols was compared to that of the molecular docking
397 models and the pharmacophores models considered separately (Figure S4, Table S5). The
398 combination protocols demonstrated either enhanced or similar performance compared to
399 individual models. For all NR except GR, the consensus protocol resulted in higher Se (Se >
400 0.69), while the hierarchical protocol is associated with increased Sp. Both combination
401 protocols also led to higher MCC values, but no clear trend favoring one protocol over the other
402 was observed.

403 For each NR, a final prediction model was selected from the 4 available options: , i.e.,
404 molecular docking, pharmacophores generation, consensus or hierarchical protocols (Table 5).
405 The Positive and Negative Likelihood Ratio (+LR and -LR, respectively) were computed for the
406 final prediction models (Table 5). The +LR values for all NRs are greater than 1, with the
407 highest +LR obtained for the TRα pharmacophore model and the lowest for the GR docking
408 model. The -LR values for all NR final prediction models are less than 1 and show similar
409 values, except for the TRα pharmacopore model, which present a higher -LR value.

410
<table><caption>Table 5: Performance of the final models selected for each NR on the whole EPA data sets</caption><thead><tr><th>NR</th><th>AR</th><th>ERα</th><th>ERβ</th><th>GR</th><th>PPARγ</th><th>TRα</th></tr><tr><th>Model</th><th>Consensus</th><th>Consensus</th><th>Consensus</th><th>Docking</th><th>Consensus</th><th>Pharmaco.</th></tr></thead><tbody><tr><th>Se</th><td>0.705</td><td>0.842</td><td>0.756</td><td>0.812</td><td>0.694</td><td>0.179</td></tr><tr><th>Sp</th><td>0.707</td><td>0.509</td><td>0.623</td><td>0.413</td><td>0.707</td><td>0.959</td></tr><tr><th>MCC</th><td>0.353</td><td>0.202</td><td>0.265</td><td>0.238</td><td>0.175</td><td>0.224</td></tr><tr><th>+LR</th><td>2.410</td><td>1.715</td><td>2.005</td><td>1.383</td><td>2.369</td><td>4.366</td></tr><tr><th>-LR</th><td>0.417</td><td>0.310</td><td>0.392</td><td>0.455</td><td>0.433</td><td>0.856</td></tr></tbody></table>

## 2.5 Final models external validations and comparison with CoMPARA and CERAPP models

The performance of the final models was evaluated using the NR-DIND validation sets (Table 5). For all NR except TRα, screening of the NR-DBIND validation sets with the final models resulted in high Se values but low Sp and MCC values. Conversely, for TRα, high Sp values were obtained whereas the remaining metrics presented low values. Focusing on the +LR and -LR metrics, ERα, ERβ and PPARγ consensus models are associated with +LR greater than 1 and -LR values less than 1, whereas AR, GR and TRα models have +LR values less than 1 and -LR values greater than 1.

The AR NR-DBIND validation set, on one side, and the ERα and ERβ NR-DBIND validation sets, on the other, were also used to evaluate the performance of the CoMPARA and CERAPP binding models, respectively (Table 6). The CoMPARA binding models were able to reach high
Se value but low Sp and MCC values. Similarly, the CERAPP models with both ERα and ERβ

NR-DBIND validation sets demonstrate moderate Sp values and low MCC values. Both

CoMPARA and CERAPP binding models were associated with +LR values less than 1 and -LR

values greater than 1.

429
<table><caption>**Table 6:** Performance of the final models selected for each NR on the NR-DBIND data sets and of the CoMPARA binding model on the AR NR-DBIND data sets and of the CERAPP binding model on the ERα and ERβ. (NaN: Not a Number)</caption><thead><tr><th>NR</th><th>Model</th><th>Se</th><th>Sp</th><th>MCC</th><th>+LR</th><th>-LR</th></tr></thead><tbody><tr><td rowspan="3">AR</td><td>Consensus</td><td>0.935</td><td>0.047</td><td>-0.030</td><td>0.981</td><td>1.383</td></tr><tr><td>Consensus with AD</td><td>0.936</td><td>0.047</td><td>-0.028</td><td>0.982</td><td>1.362</td></tr><tr><td>CoMPARA</td><td>0.698</td><td>0.087</td><td>-0.198</td><td>0.765</td><td>3.471</td></tr><tr><td rowspan="3">ERα</td><td>Consensus</td><td>0.983</td><td>0.024</td><td>0.024</td><td>1.007</td><td>0.708</td></tr><tr><td>Consensus with AD</td><td>1</td><td>0.017</td><td>0.114</td><td>1.017</td><td>0</td></tr><tr><td>CERAPP</td><td>0.434</td><td>0.402</td><td>-0.139</td><td>0.726</td><td>1.408</td></tr><tr><td rowspan="3">ERβ</td><td>Consensus</td><td>0.895</td><td>0.432</td><td>0.291</td><td>1.576</td><td>0.243</td></tr><tr><td>Consensus with AD</td><td>0.925</td><td>0.471</td><td>0.371</td><td>1.749</td><td>0.159</td></tr><tr><td>CERAPP</td><td>0.530</td><td>0.457</td><td>0.045</td><td>0.976</td><td>1.028</td></tr><tr><td rowspan="2">GR</td><td>Docking</td><td>0.998</td><td>0.0</td><td>-0.011</td><td>0.998</td><td>NaN</td></tr><tr><td>Docking with AD</td><td>1</td><td>0</td><td>NaN</td><td>1.000</td><td>NaN</td></tr><tr><td rowspan="2">PPARγ</td><td>Consensus</td><td>0.998</td><td>0.042</td><td>0.153</td><td>1.042</td><td>0.048</td></tr><tr><td>Consensus with AD</td><td>1</td><td>0.012</td><td>0.094</td><td>1.012</td><td>0</td></tr><tr><td rowspan="2">TRα</td><td>Pharmacophore</td><td>0.090</td><td>0.786</td><td>-0.138</td><td>0.421</td><td>1.158</td></tr><tr><td>Pharmacophore with AD</td><td>0.108</td><td>0.75</td><td>-0.098</td><td>0.432</td><td>1.189</td></tr></tbody></table>

## 2.6 Applicability domain

The AD was defined separately for each NR using molecular descriptors computed with the RDkit python package on the EPA data sets compounds on one side and the NR-DBIND data sets compounds on the other side. The performance of the final prediction models was re-
evaluated for each NR, considering only the NR-DBIND compounds that fell within the defined AD (Table 6). The composition of the NR-DBIND subsets falling within the AD is described in the Table 2. Applying the AD led to an improvement in Se and MCC values for all NRs, except for GR and PPARγ in MCC. Meanwhile, Sp increased for ERβ but decreased for ERα, PPARγ, and TRα.

## 3. Discussion

### 3.1 Data

For this study, we decided to use the chemistry and associated toxicity data from the CompTox Chemicals Dashboard<sup>10</sup>, designed to provide data and models to prioritize chemicals for experimental testing and subsequently reduce the use of animals in chemical testing. In particular, we used the ToxCast (EPA's Toxicity Forecaster) campaign outcomes obtained using high-throughput screening assays. This choice was guided by a previous review of the *in silico* studies dedicated to the NR family<sup>5</sup> highlighting a large number of studies published between 2010 and 2020 conducted in a toxicological context using the ToxCast data for NR models prediction generation. Among all the compounds included in this dashboard, we selected and extracted only the compounds associated with binding assays outcomes for one of the 6 NR targeted by this study—to ensure gathering compounds sharing the same NR modulation mechanism of action. Compounds evaluated only in other NR associated assays such as agonism and antagonism evaluation assays may modulate upstream or downstream NR signaling pathways and not acting through direct binding to NR. The selection of compounds tested on binding assays is particularly important for docking methods that are used to predict the conformation and orientation of a compound while binding a protein. Similarly, pharmacophores
generation requires a training data set of compounds presenting a similar mechanism of action to identify the part of the structures responsible for the biological activity.

Prediction models quality is directly linked to the data used to construct the models and several bias can be encountered along the training and test sets construction protocols³⁶. We thus analyzed: (1) the distribution of the physicochemical properties between binding and non-binding compounds, (2) the promiscuity of the binding compounds included in each NR sets and (3) the balance between the number of binding and non-binding compounds.

We conducted a PCA for each of the 6 EPA data sets. PCA is a statistical technique used to reduce the dimensionality of a data set containing various interrelated variables, while preserving the maximum of the variation in the data set ³⁷. It is frequently used for data visualization, in particular to compare the chemical space coverage among different data sets based on structural and physico-chemical descriptors ³⁸. The results obtained for the 6 EPA data sets (Figure 3, Table S2) show that both binding and non-binding compounds present a similar chemical space coverage according to the studied descriptors. Even considering 50 principal components to account for at least 80% of the variability, the clustering analysis showed that binding and non-binding compounds cannot be separated. The physicochemical descriptors used to conduct the PCA analysis are therefore not sufficient to explain the NR binding ability of these compounds. Moreover, this analysis also highlights that the prediction models generated using the EPA data sets will not be biased according to a diversity in the considered physico-chemical properties between binding and non-binding compounds.

The redundancy of binding compounds among the 6 NR data sets (Figure 4) can be explained by the ToxCast protocol in which the same pool of chemicals, associated with toxicological concerns for human and ecosystem exposure, was tested in hundreds of high-throughput
screening assays. Additionally, NR ligands selectivity has been largely studied in drug discovery campaigns ³⁹, especially for different NR isoforms such as ERα and ERβ. These two isoforms display high LBD similarity, and several compounds are able to bind to both isotypes with comparable affinities, including natural ER ligands⁴⁰. Natural selective ligands such as the phytoestrogen genistein ERβ agonist ligand⁴¹ or synthetic ones are also frequent. Surprisingly, the ERα and ERβ EPA data sets share only 68 common binding compounds, including 26 not found in any other EPA data sets. The ERα data set is even the second one with the largest number of unique binding compounds. In the same way, very few compounds are promiscuous, since no binding compound is present in the 6 data sets and less than 15% of the binding compounds are able to bind to more than 2 NR.

In all the EPA data set, except for GR, the number of non-binding compounds is largely superior to the number of binding compounds. An opposite trend is observed for the 6 NR-DBIND validation sets, in which the number of binding compounds outnumbered the number of non-binding compounds. This imbalance must be considered for the analysis of the performance metrics. Similarly, the PPARγ EPA data set contains a relatively low number of binding compounds whereas the corresponding NR-DBIND data set includes a high number of binding compounds. This difference is easily explained by the fact that numerous synthetic ligands have been produced for PPARγ, a well-known therapeutic target for many diseases⁴², but its involvement in some EDC chemotypes binding has been discovered more recently.

## 3.2 Docking protocol

Different approaches, search algorithms, scoring functions and protein structures can be used to conduct protein-ligand molecular docking studies, and retrospective evaluation of molecular docking protocols must be conducted to ensure its performance for each target. In a previous
study conducted on $ER\alpha^{13}$, we identified among 3 molecular docking software (smina, PLANTS and Surflex-dock) smina as the one the most performant to discriminate binding and non-binding compounds of $ER\alpha$. We thus decided to also use smina in this study, that is extended to 5 other NR. For each one of the 6 studied NR, an adequate and specific molecular docking protocol was selected by evaluating the molecular docking results obtained for all structures with 4 scoring functions and 2 approaches called single structure and ensemble structures molecular docking.

No or limited improvement of the best AUC values is associated with the ensemble structures molecular docking approach. Considering the higher computational times associated with this approach, we selected only the single structure molecular docking approach for our prediction models.

We observed that no scoring function is associated with the best overall performance in binding compounds prediction for all NR. For each NR, except for $TR\alpha$, we identified the scoring function and the PDB structure associated with the highest AUC. The scoring function associated with the best AUC was `dkoes_scoring` for 3 NR and `vina` and `vinardo` for 1 NR. In most cases, using the same scoring function, the amplitude of AUC variation observed between the different PDB structures is low. As the NR ligands (and EDC in general) included in the training sets present different pharmacological profiles, no clear trend was observed for the pharmacological profile of the ligand bound in the selected PDB structure.

Regarding $TR\alpha$, all the studied molecular docking protocols failed to discriminate binding and non-binding compounds. For this NR, only 4 PDB structures matching our inclusion criteria (see Material and Methods) were available in the PDB at the time of this study, all in complex with a $TR\beta$ specific ligand. The conformation of the query protein can impact the molecular docking performance$^{43-45}$ and this might explain the very low prediction performance obtained for $TR\alpha$.
We decided to exclude all the studied molecular docking protocols from the TRα prediction models. Additional approaches taking into account the protein flexibility s or newly solved TRα structures may enable filling the gap for this NR in the future.

The molecular docking protocol selected for each NR (except TRα) was completed with a scoring threshold that will be used to predict compounds as binding compounds (molecular docking scores above the threshold) and non-binding compounds (molecular docking scores below the threshold). The threshold was defined for each NR using the Screening Explorer tool by selecting the value maximizing both sensitivity and predictivity values. The thresholds retained for each NR are relatively comparable (from -6 to -7), probably because the proteins all belong to a same family and share some binding site similarities.

Alternative strategy for ensemble molecular docking results analysis ⁴⁶ and other protocols such as molecular docking methods taking into account the flexibility of the protein with limited impact on the computational cost⁴⁷,⁴⁸, or interaction fingerprint approaches⁴⁹ could be evaluated in future studies to improve the docking prediction performance of our models.

### 3.3 Pharmacophore models

In a previous study dedicated to ERα¹³, we presented an optimization protocol designed to reduce the large number of SB and LB pharmacophore models initially generated. This optimization step is important to select the pharmacophores associated with the highest prediction performance but also to reduce the computational time that increases with the number of pharmacophores. In this study, the optimization protocol has been modified (see Material and Methods) to allow its automation. During this optimization protocol, we arbitrarily fixed a binding compounds/non-binding compounds ratio equal or superior to 1 to select optimized
pharmacophores. A few other higher ratio values were tested during the development of the
optimization protocol, but the prediction performance associated with the selected
pharmacophores was decreased. An alternative to this binding compounds /non-binding
compounds ratio called the growth rate was used in the study of Métivier et al.⁵⁰ by computing
the ratio of the pharmacophore frequency of fits within the active molecules and within the
inactive ones. This alternative could be evaluated to improve the optimization process in future
studies.

The final optimization protocol used in this study enables dividing the total number of
pharmacophores included in the prediction models by a factor 13 (Table 4). The impact of this
optimization step on the pharmacophores binding compounds performance varies among SB and
LB pharmacophores and among NR except for the MCC values that are enhanced for all the
optimized pharmacophores. Additionally, improving one of the two metrics, Se or Sp, always
negatively impacts the second.

For each NR, the optimized SB and LB pharmacophores were then combined in an SBLB
pharmacophore models set. A total of 74 pharmacophore models are included in these SBLB
pharmacophore models sets.

## 3.4 Combinations

The models presented in this study were built using a combination of structure-based and ligand-
based methods. Ligand-based methods are associated with an applicability domain⁵¹, i.e., their
usage is often limited to provide prediction for compounds similar to those used in the model
training data set. Conversely, structure-based method, such as the docking method used in this
study, are independent from the structure of the compounds in the training sets. Combining
structure-based and ligand-based approaches can help overcome applicability domain constraints
while leveraging the strengths of each method and mitigating their respective limitations⁻⁵²,⁵³.

However, we decided to compute the applicability domain of our final prediction models by combining a bounding-box approach and a Mahalanobis distance-based approach.

Different strategies have been developed to combine SB and LB methods⁵³. In a previous study¹³, we introduced two protocols to combine the molecular docking and pharmacophore models for ERα, which we named consensus and hierarchical protocols. These two combinations protocols were also applied in the present study, except for TRα, for which no molecular docking protocol was retained due to very low prediction performance. For the 5 remaining NR, neither protocol led to a simultaneous improvement in all performance metrics. For example, we noticed for example that a high Se often comes at the expense of a lower Sp and reciprocally. The predictive models developed in this study aim to be used for prediction of compounds able to bind NR in a toxicological context, that is to flag compounds that should be further investigated for their potential toxicity and EDCs potency. Therefore, models selection was based on their ability to minimize the number of false negative and correctly identify binding compounds ³³,⁵⁴, i.e. prioritizing low -LR and high Se values.

The consensus protocol was selected for AR, ERα, ERβ and PPARγ as it yields the highest Se, and MCC values and the lowest -LR values. For GR, the performance of the hierarchical protocol was similar-to those of the molecular docking models. We chose to select the molecular docking protocol over the hierarchical one for GR since the SBLB pharmacophore models exhibit a high permissiveness by retrieving 371 hits from the 398 initial compounds (binding and non-binding compounds) including in the EPA data set.

The -LR values obtained for each NR using the EPA datasets highlight the models' ability to correctly identify non-binding compounds, aligning with our objective in the toxicological
context to minimize false negatives. +LR values indicate a moderate informative power of the NR models in predicting a positive result for all NRs, except for TRα, which is associated with the highest +LR value.

Comparison of the distribution of the pharmacological profile among the binding compounds and the final prediction models hitlists shows that the final prediction models are not biased towards a specific pharmacological profile except for GR (Table S6). Specifically, for GR, no agonist compound was predicted as a binding compound using the final prediction model, which relies solely on the docking protocol. This may be explained by the fact that the GR structure used in the docking protocol is in an antagonist bound conformation.

## 3.5 Final models validations and comparison with other NR dedicated prediction models

6 external data sets generated from the NR-DBIND and dedicated to AR, ERα, ERβ, GR, PPARγ and TRα were used to validate the selected protocols for each NR. In comparison to the performance obtained with the EPA data sets, screening of the NR-DBIND data sets with the selected protocols yielded higher Se values but drastically decreased Sp values except for TRα. The +LR and -LR values suggest that the AR, GR, and TRα predictions may not be informative on the NR-DBIND datasets. For ERα, ERβ, and PPARγ, the low -LR values but moderate +LR values indicate that these models are useful for identifying non-binding compounds within the NR-DBIND datasets but less effective for predicting binding compounds. An explanation in the difference of performance obtained with both data sets can be found in the distinct chemical space covered by the compounds of the EPA and NR-DBIND sets (Table S2). This was expected as the NR-DBIND was constructed by gathering drug-like molecules while the EPA database
contains structurally diverse environmental compounds originating from different sectors such as food and agriculture or cosmetic industries.

Various NR prediction models have been described in the literature. Recently, we published a review summarizing 89 *in silico* studies on NR ligands binding prediction, spanning from 2010 to 2020. In this review, we described several webservers and software, such as OpenVirtualToxLab⁵⁵ or ProTox-II⁵⁶ designed to provide predictions on ensemble of NR. However, the use of these tools is often constrained by their single compound input limitation.

We also highlighted that LB methods, particularly QSAR methods, employing machine learning and deep learning algorithms, demonstrate high efficiency and can predict several endpoints related to NR and EDC. This is notably exemplified in the case of AR and ER, for which the CoMPARA⁵⁷ and CERAPP⁵⁸ international consortium, respectively, have provided consensus models for ligand binding prediction, including QSAR models. It is to note that both CoMPARA and CERAPP binding models integrate a combination of QSAR and docking approaches. In the CoMPARA binding model, the docking methodology relies on the use of the GOLD software and the 2PNU AR structure. Interestingly, the 2PNU structure was also associated with the best outcomes in our study using the smina software. The CoMPARA binding models exhibit very high Se and Sp on their training set (>0.95). However, on the evaluation set, the Se is lower (0.65), while the Sp remains high (0.90).⁵⁷ The CERAPP binding models include two different docking protocols: one using the GLIDE software and the 1L2I ERα structure, and another using Autodock with two ERα structures, 1GWR and 3ERT. While the CERAPP binding models achieve very high Se and Sp values on their training dataset (>0.93), they exhibit a lower Se on the literature evaluation set.⁵⁸
Comparing the performance of our models with previously described models is challenging due to variations in training data sets (in terms of size, composition and representativeness), differences in the level of complexity of the models and difficulty gathering an external validation sets with compounds not included in one of the training data sets. To provide a preliminary comparison between our final models and reference models, we screened the AR NR-DBIND set with the CoMPARA binding models⁵⁷ and the ERα and ERβ NR-DBIND sets with the CERAPP binding models⁵⁸. The binding prediction performance obtained from the CoMPARA binding model and our AR final consensus model are very similar, with the AR final consensus models slightly outperforming the CoMPARA models across all evaluated performance metrics. When comparing the CERAPP binding models to our ERα and ERβ consensus models, the former presented higher Sp values (dramatically for ERα) but was outperformed in all other metrics, particularly the likelihood ratios. However, this comparison is inherently biased, since our models were specifically designed for ERα or ERβ ligands binding prediction, whereas CERAPP binding models were constructed to assess the overall ER pathway by merging both isoforms. This preliminary and limited comparative study shows that our consensus model dedicated to AR ligand binding prediction performed at least as well as the CoMPARA binding model using the NR-DBIND sets.

To complement the external validation of the models presented in this study, we conducted a case study demonstrating their application. For this, we combined two publicly available lists of suspected EDCs but without precise biological activity annotations or confirmed molecular targets, the EDlist (https://edlists.org/the-ed-lists) and The Endocrine Disruption eXchange list (TEDX list, https://endocrinedisruption.org/interactive-tools/tedx-list-of-potential-endocrine-disruptors/search-the-tedx-list). After removing compounds already present in the EPA data sets,
we obtained 532 compounds which we screened using our final prediction models (consensus protocol) for AR, ERα, ERβ and PPARγ. On the 532 compounds, 467 were flagged as binding compounds for at least one of the 4 studied NR, and only 65 were predicted as non-binding compounds for all 4 NRs (Table S7).

Finally, the herein presented *in silico* protocols only include one SB and one LB method. However, recent advancements in machine learning and deep learning algorithms have significantly improved predictive performance in QSAR-based approaches <sup>59</sup>, including for NR ligand binding prediction<sup>5</sup>. As a next step, future studies should compare the results obtained using our combination of molecular docking and pharmacophore models with those from QSAR-based methods. Additionally, integrating all three approaches could further enhance predictive accuracy and robustness.

## 4. Conclusion

Rapidly and efficiently identifying harmful compounds to human health and ecosystems equilibrium has become a major challenge in recent years. Endocrine disrupting chemicals (EDC), a class of compounds associated with numerous adverse health effects, are of particular concern within the human exposome. While experimental assays remain essential, in silico methods, especially molecular docking and pharmacophore modeling, offer powerful tools to predict compound interactions with nuclear receptors (NR), modeling their “direct” mechanism of action. In this study, we developed and optimized *in silico* protocols combining molecular docking and pharmacophore modelling approaches for 6 NR that have been linked to the “direct” mechanism action of EDC. These protocols can be used, in complement to experimental assays, to prioritize compounds that should be further evaluated to detect any toxicological concern.
688
689 The authors declare that they have no known competing financial interests or personal
690 relationships that could have appeared to influence the work reported in this paper.

691
692 Acknowledgment

693 We would like to thank Prof. T. Langer and Inte:Ligand for the LigandScout 4.5 software
694 license. We would like to thank Taoufik Labib for his assistance for cluster computing.

695
696 A.S. is recipient of a MESRI fellowship

697
698 Abbreviations

699 AR, androgen receptor; AUC, area under the ROC curve; clogP, octanol-water partition
700 coefficient; clogS, solubility; DNA, desoxyribonucleic acid; DBD, desoxyribonucleic acid
701 binding domain; EDCs, endocrine disrupting chemicals; EPA, environmental protection agency;
702 ERα, estrogen receptor alpha; ERβ, estrogen receptor beta; FN, false negative; FP, false positive;
703 GR, glucocorticoid receptor; HBA, hydrogen bond acceptors; HBD, hydrogen bond donors; LB,
704 ligand-based; LBD, ligand binding domain; MCC, Matthew correlation coefficient; MW,
705 molecular weight; NR, Nuclear Receptors; OMF, omitted features; PCA, principal component
706 analysis; PDB, Protein Data Bank; PPARγ, peroxisome proliferator-activated receptor gamma;
707 PSA, polar surface area; SB, structure-based; Se, sensitivity; Sp, specificity; TN, true negative;
708 ToxCast, EPA's Toxicity Forecaster; TP, true positive; TRα, thyroid hormone receptor alpha;
709 TSA, total surface area; t-SNE, t-distributed stochastic neighbor embedding; vinardo, vina radii
710 optimized.

711

712 **Supporting information**

713 pharmacological profile distribution of binding compounds of the EPA data sets (Figure S1); plot
714 of metrics values according as a function of molecular docking scores (Figure S2); radar charts
715 representation of the metric values of generated and optimized LB, SB and SBLB
716 pharmacophores (Figure S3); radar charts representation of the metric values of docking, SBLB,
717 consensus and hierarchical protocols (Figure S4); mathematical equations of the evaluation
718 metrics (Table S1); binding and non binding compounds distribution across the two clusters
719 computed from the first 50 principal components of the PCA analysis (Table S2); AUC values of
720 all evaluated docking approaches (Table S3); prediction performance of the generated and
721 optimized pharmacophores (Table S4); prediction performance of the docking, SBLB
722 pharmacophores, consensus and hierarchical protocols (Table S5); pharmacological profile of the
723 predicted binding compounds (Table S6); screening results of the compounds extract from the
724 EP list and TEDX (Table S7).

725

726 **Data statement**

727 The data set supporting the conclusions of this article are available in the zenodo repository
728 https://doi.org/10.5281/zenodo.10472286 (the password to unzip the repository is
729 Gt_S<'3x;n5Ykhv6)
The ScreeNR.py script and the prediction models supporting the conclusions of this article are available in the zenodo repository https://doi.org/10.5281/zenodo.14974415 (the password to unzip the repository is VWNy707nFiaEijC)

(1) Imai, Y.; Youn, M.-Y.; Inoue, K.; Takada, I.; Kouzmenko, A.; Kato, S. Nuclear Receptors in Bone Physiology and Diseases. *Physiol. Rev.* **2013**, *93* (2), 481–523.
https://doi.org/10.1152/physrev.00008.2012.

(2) Weikum, E. R.; Liu, X.; Ortlund, E. A. The Nuclear Receptor Superfamily: A Structural Perspective. *Protein Science* **2018**, *27* (11), 1876–1892. https://doi.org/10.1002/pro.3496.

(3) Overington, J. P.; Al-Lazikani, B.; Hopkins, A. L. How Many Drug Targets Are There? *Nat Rev Drug Discov* **2006**, *5* (12), 993–996. https://doi.org/10.1038/nrd2199.

(4) Yang, Z.; Gimple, R. C.; Zhou, N.; Zhao, L.; Gustafsson, J.-Å.; Zhou, S. Targeting Nuclear Receptors for Cancer Therapy: Premises, Promises, and Challenges. *Trends in Cancer* **2021**, *7* (6), 541–556. https://doi.org/10.1016/j.trecan.2020.11.007.

(5) Sellami, A.; Réau, M.; Montes, M.; Lagarde, N. Review of in Silico Studies Dedicated to the Nuclear Receptor Family: Therapeutic Prospects and Toxicological Concerns. *Frontiers in Endocrinology* **2022**, *13*.

(6) Toporova, L.; Balaguer, P. Nuclear Receptors Are the Major Targets of Endocrine Disrupting Chemicals. *Molecular and Cellular Endocrinology* **2020**, *502*, 110665.
https://doi.org/10.1016/j.mce.2019.110665.

(7) International Programme on Chemical Safety. *Global Assessment on the State of the Science of Endocrine Disruptors*; WHO/PCS/EDC/02.2; World Health Organization, 2002.
https://apps.who.int/iris/handle/10665/67357 (accessed 2023-03-27).
753 (8) Lee, H.-R.; Jeung, E.-B.; Cho, M.-H.; Kim, T.-H.; Leung, P. C. K.; Choi, K.-C. Molecular Mechanism(s) of Endocrine-Disrupting Chemicals and Their Potent Oestrogenicity in Diverse Cells and Tissues That Express Oestrogen Receptors. *Journal of Cellular and Molecular Medicine* **2013**, *17* (1), 1–11. https://doi.org/10.1111/j.1582-4934.2012.01649.x.

757 (9) *Endocrine disruptor assessment list* - ECHA. https://echa.europa.eu/da/ed-assessment
758 (accessed 2023-03-27).

759 (10) Williams, A. J.; Grulke, C. M.; Edwards, J.; McEachran, A. D.; Mansouri, K.; Baker, N. C.; Patlewicz, G.; Shah, I.; Wambaugh, J. F.; Judson, R. S.; Richard, A. M. The CompTox Chemistry Dashboard: A Community Data Resource for Environmental Chemistry. *Journal of Cheminformatics* **2017**, *9* (1), 61. https://doi.org/10.1186/s13321-017-0247-6.

763 (11) Casals-Casas, C.; Desvergne, B. Endocrine Disruptors: From Endocrine to Metabolic Disruption. *Annual Review of Physiology* **2011**, *73* (1), 135–162. https://doi.org/10.1146/annurev-physiol-012110-142200.

766 (12) Devillers, J.; Marchand-Geneste, N.; Carpy, A.; Porcher, J. M. SAR and QSAR Modeling of Endocrine Disruptors. *SAR and QSAR in Environmental Research* **2006**, *17* (4), 393–412. https://doi.org/10.1080/10629360600884397.

769 (13) Sellami, A.; Montes, M.; Lagarde, N. Predicting Potential Endocrine Disrupting Chemicals Binding to Estrogen Receptor α (ERα) Using a Pipeline Combining Structure-Based and Ligand-Based in Silico Methods. *Int J Mol Sci* **2021**, *22* (6), 2846. https://doi.org/10.3390/ijms22062846.

773 (14) Richard, A. M.; Judson, R. S.; Houck, K. A.; Grulke, C. M.; Volarath, P.; Thillainadarajah, I.; Yang, C.; Rathman, J.; Martin, M. T.; Wambaugh, J. F.; Knudsen, T. B.; Kancherla, J.; Mansouri, K.; Patlewicz, G.; Williams, A. J.; Little, S. B.; Crofton, K. M.; Thomas, R. S.
776 ToxCast Chemical Landscape: Paving the Road to 21st Century Toxicology. *Chem. Res.*

777 *Toxicol.* **2016**, *29* (8), 1225–1251. https://doi.org/10.1021/acs.chemrestox.6b00135.

778 (15) Réau, M.; Lagarde, N.; Zagury, J.-F.; Montes, M. Nuclear Receptors Database Including

779 Negative Data (NR-DBIND): A Database Dedicated to Nuclear Receptors Binding Data

780 Including Negative Data and Pharmacological Profile. *J Med Chem* **2019**, *62* (6), 2894–2904. https://doi.org/10.1021/acs.jmedchem.8b01105.

782 (16) *ChemAxon* - *Software Solutions and Services for Chemistry & Biology*.

783 https://chemaxon.com/ (accessed 2020-07-27).

784 (17) Wolber, G.; Langer, T. LigandScout: 3-D Pharmacophores Derived from Protein-Bound

785 Ligands and Their Use as Virtual Screening Filters. *J Chem Inf Model* **2005**, *45* (1), 160–

786 169. https://doi.org/10.1021/ci049885e.

787 (18) RDKit: Open-Source Cheminformatics. https://www.rdkit.org.

788 (19) Pedregosa, F.; Varoquaux, G.; Gramfort, A.; Michel, V.; Thirion, B.; Grisel, O.; Blondel,

789 M.; Prettenhofer, P.; Weiss, R.; Dubourg, V.; Vanderplas, J.; Passos, A.; Cournapeau, D.;

790 Brucher, M.; Perrot, M.; Duchesnay, É. Scikit-Learn: Machine Learning in Python. *Journal*

791 *of Machine Learning Research* **2011**, *12* (85), 2825–2830.

792 (20) Pettersen, E. F.; Goddard, T. D.; Huang, C. C.; Couch, G. S.; Greenblatt, D. M.; Meng, E.

793 C.; Ferrin, T. E. UCSF Chimera--a Visualization System for Exploratory Research and

794 Analysis. *J Comput Chem* **2004**, *25* (13), 1605–1612. https://doi.org/10.1002/jcc.20084.

795 (21) Koes, D. R.; Baumgartner, M. P.; Camacho, C. J. Lessons Learned in Empirical Scoring

796 with Smina from the CSAR 2011 Benchmarking Exercise. *J Chem Inf Model* **2013**, *53* (8),

797 1893–1904. https://doi.org/10.1021/ci300604z.
798 (22) Trott, O.; Olson, A. J. AutoDock Vina: Improving the Speed and Accuracy of Docking with
799 a New Scoring Function, Efficient Optimization and Multithreading. *J Comput Chem* **2010**,
800 31 (2), 455–461. https://doi.org/10.1002/jcc.21334.

801 (23) Quiroga, R.; Villarreal, M. A. Vinardo: A Scoring Function Based on Autodock Vina
802 Improves Scoring, Docking, and Virtual Screening. *PloS One* **2016**, *11* (5), e0155183.
803 https://doi.org/10.1371/journal.pone.0155183.

804 (24) Réau, M.; Lagarde, N.; Zagury, J.-F.; Montes, M. Hits Discovery on the Androgen
805 Receptor: In Silico Approaches to Identify Agonist Compounds. *Cells* **2019**, *8* (11), 1431.
806 https://doi.org/10.3390/cells8111431.

807 (25) Empereur-mot, C.; Guillemain, H.; Latouche, A.; Zagury, J.-F.; Viallon, V.; Montes, M.
808 Predictiveness Curves in Virtual Screening. *Journal of Cheminformatics* **2015**, *7* (1), 52.
809 https://doi.org/10.1186/s13321-015-0100-8.

810 (26) Empereur-Mot, C.; Zagury, J.-F.; Montes, M. Screening Explorer-An Interactive Tool for
811 the Analysis of Screening Results. *J Chem Inf Model* **2016**, *56* (12), 2281–2286.
812 https://doi.org/10.1021/acs.jcim.6b00283.

813 (27) Wermuth, C. G.; Ganellin, C. R.; Lindberg, P.; Mitscher, L. A. Glossary of Terms Used in
814 Medicinal Chemistry (IUPAC Recommendations 1998). *Pure and Applied Chemistry* **1998**,
815 70 (5), 1129–1143. https://doi.org/10.1351/pac199870051129.

816 (28) Seidel, T.; Wieder, O.; Garon, A.; Langer, T. Applications of the Pharmacophore Concept
817 in Natural Product Inspired Drug Design. *Molecular Informatics* **2020**, *39* (11), 2000059.
818 https://doi.org/10.1002/minf.202000059.
819 (29) Seidel, T.; Ibis, G.; Bendix, F.; Wolber, G. Strategies for 3D Pharmacophore-Based Virtual
820 Screening. *Drug Discovery Today: Technologies* **2010**, *7* (4), e221-e228.
821 https://doi.org/10.1016/j.ddtec.2010.11.004.

822 (30) Vuorinen, A.; Schuster, D. Methods for Generating and Applying Pharmacophore Models
823 as Virtual Screening Filters and for Bioactivity Profiling. *Methods* **2015**, *71*, 113-134.
824 https://doi.org/10.1016/j.ymeth.2014.10.013.

825 (31) Lagarde, N.; Delahaye, S.; Zagury, J.-F.; Montes, M. Discriminating Agonist and
826 Antagonist Ligands of the Nuclear Receptors Using 3D-Pharmacophores. *J Cheminform*
827 **2016**, *8* (1), 43. https://doi.org/10.1186/s13321-016-0154-2.

828 (32) OECD. *Guidance Document on the Validation of (Quantitative) Structure-Activity
829 Relationship [(Q)SAR] Models*; OECD Series on Testing and Assessment; OECD, 2014.
830 https://doi.org/10.1787/9789264085442-en.

831 (33) Trisciuzzi, D.; Alberga, D.; Mansouri, K.; Judson, R.; Novellino, E.; Mangiatordi, G. F.;
832 Nicolotti, O. Predictive Structure-Based Toxicology Approaches To Assess the Androgenic
833 Potential of Chemicals. *J Chem Inf Model* **2017**, *57* (11), 2874-2884.
834 https://doi.org/10.1021/acs.jcim.7b00420.

835 (34) Sahigara, F.; Mansouri, K.; Ballabio, D.; Mauri, A.; Consonni, V.; Todeschini, R.
836 Comparison of Different Approaches to Define the Applicability Domain of QSAR Models.
837 *Molecules* **2012**, *17* (5), 4791-4810. https://doi.org/10.3390/molecules17054791.

838 (35) Sahigara, F.; Mansouri, K.; Ballabio, D.; Mauri, A.; Consonni, V.; Todeschini, R.
839 Comparison of Different Approaches to Define the Applicability Domain of QSAR Models.
840 *Molecules* **2012**, *17* (5), 4791-4810. https://doi.org/10.3390/molecules17054791.
841 (36) Sieg, J.; Flachsenberg, F.; Rarey, M. In Need of Bias Control: Evaluating Chemical Data for Machine Learning in Structure-Based Virtual Screening. *J. Chem. Inf. Model.* **2019**, *59*
842 (3), 947–961. https://doi.org/10.1021/acs.jcim.8b00712.

844 (37) Introduction. In *Principal Component Analysis*; Jolliffe, I. T., Ed.; Springer Series in Statistics; Springer: New York, NY, 2002; pp 1–9. https://doi.org/10.1007/0-387-22440-8_1.

847 (38) Akella, L. B.; DeCaprio, D. Cheminformatics Approaches to Analyze Diversity in Compound Screening Libraries. *Current Opinion in Chemical Biology* **2010**, *14* (3), 325–330. https://doi.org/10.1016/j.cbpa.2010.03.017.

850 (39) Burris, T. P.; Solt, L. A.; Wang, Y.; Crumbley, C.; Banerjee, S.; Griffett, K.; Lundasen, T.; Hughes, T.; Kojetin, D. J. Nuclear Receptors and Their Selective Pharmacologic Modulators. *Pharmacol Rev* **2013**, *65* (2), 710–778. https://doi.org/10.1124/pr.112.006833.

853 (40) Mal, R.; Magner, A.; David, J.; Datta, J.; Vallabhaneni, M.; Kassem, M.; Manouchehri, J.; Willingham, N.; Stover, D.; Vandeusen, J.; Sardesai, S.; Williams, N.; Wesolowski, R.; Lustberg, M.; Ganju, R. K.; Ramaswamy, B.; Cherian, M. A. Estrogen Receptor Beta (ERβ): A Ligand Activated Tumor Suppressor. *Frontiers in Oncology* **2020**, *10*.

857 (41) Shanle, E. K.; Xu, W. Endocrine Disrupting Chemicals Targeting Estrogen Receptor Signaling: Identification and Mechanisms of Action. *Chem. Res. Toxicol.* **2011**, *24* (1), 6–19. https://doi.org/10.1021/tx100231n.

860 (42) Mirza, A. Z.; Althagafi, I. I.; Shamshad, H. Role of PPAR Receptor in Different Diseases and Their Ligands: Physiological Importance and Clinical Implications. *Eur J Med Chem* **2019**, *166*, 502–513. https://doi.org/10.1016/j.ejmech.2019.01.067.
863 (43) Hawkins, P. C. D.; Warren, G. L.; Skillman, A. G.; Nicholls, A. How to Do an Evaluation: Pitfalls and Traps. *J Comput Aided Mol Des* **2008**, *22* (3), 179–190.
864
865 https://doi.org/10.1007/s10822-007-9166-3.

866 (44) Ben Nasr, N.; Guillemain, H.; Lagarde, N.; Zagury, J.-F.; Montes, M. Multiple Structures for Virtual Ligand Screening: Defining Binding Site Properties-Based Criteria to Optimize the Selection of the Query. *J Chem Inf Model* **2013**, *53* (2), 293–311.
867
868
869 https://doi.org/10.1021/ci3004557.

870 (45) Wang, N.; Wang, L.; Xie, X.-Q. ProSelection: A Novel Algorithm to Select Proper Protein Structure Subsets for in Silico Target Identification and Drug Discovery Research. *J. Chem. Inf. Model.* **2017**, *57* (11), 2686–2698. https://doi.org/10.1021/acs.jcim.7b00277.

873 (46) Ricci-Lopez, J.; Aguila, S. A.; Gilson, M. K.; Brizuela, C. A. Improving Structure-Based Virtual Screening with Ensemble Docking and Machine Learning. *J. Chem. Inf. Model.* **2021**, *61* (11), 5362–5376. https://doi.org/10.1021/acs.jcim.1c00511.

876 (47) Spiriti, J.; Subramanian, S. R.; Palli, R.; Wu, M.; Zuckerman, D. M. Middle-Way Flexible Docking: Pose Prediction Using Mixed-Resolution Monte Carlo in Estrogen Receptor α. *PLOS ONE* **2019**, *14* (4), e0215694. https://doi.org/10.1371/journal.pone.0215694.

879 (48) Wong, K. M.; Tai, H. K.; Siu, S. W. I. GWOVina: A Grey Wolf Optimization Approach to Rigid and Flexible Receptor Docking. *Chemical Biology & Drug Design* **2021**, *97* (1), 97–110. https://doi.org/10.1111/cbdd.13764.

882 (49) Jasper, J. B.; Humbeck, L.; Brinkjost, T.; Koch, O. A Novel Interaction Fingerprint Derived from per Atom Score Contributions: Exhaustive Evaluation of Interaction Fingerprint Performance in Docking Based Virtual Screening. *J Cheminform* **2018**, *10* (1), 15.
883
884
885 https://doi.org/10.1186/s13321-018-0264-0.
886 (50) Métivier, J.-P.; Cuissart, B.; Bureau, R.; Lepailleur, A. The Pharmacophore Network: A
887 Computational Method for Exploring Structure-Activity Relationships from a Large
888 Chemical Data Set. *J. Med. Chem.* **2018**, *61* (8), 3551–3564.
889 https://doi.org/10.1021/acs.jmedchem.7b01890.

890 (51) Sahigara, F.; Mansouri, K.; Ballabio, D.; Mauri, A.; Consonni, V.; Todeschini, R.
891 Comparison of Different Approaches to Define the Applicability Domain of QSAR Models.
892 *Molecules* **2012**, *17* (5), 4791–4810. https://doi.org/10.3390/molecules17054791.

893 (52) Sperandio, O.; Miteva, M. A.; Villoutreix, B. O. Combining Ligand- and Structure-Based
894 Methods in Drug Design Projects. *Current Computer-Aided Drug Design* **2008**, *4* (3), 250–
895 258.

896 (53) Vázquez, J.; López, M.; Gibert, E.; Herrero, E.; Luque, F. J. Merging Ligand-Based and
897 Structure-Based Methods in Drug Discovery: An Overview of Combined Virtual Screening
898 Approaches. *Molecules* **2020**, *25* (20), 4723. https://doi.org/10.3390/molecules25204723.

899 (54) Vuorinen, A.; Odermatt, A.; Schuster, D. In Silico Methods in the Discovery of Endocrine
900 Disrupting Chemicals. *The Journal of Steroid Biochemistry and Molecular Biology* **2013**,
901 137, 18–26. https://doi.org/10.1016/j.jsbmb.2013.04.009.

902 (55) Vedani, A.; Dobler, M.; Hu, Z.; Smieško, M. OpenVirtualToxLab—A Platform for
903 Generating and Exchanging in Silico Toxicity Data. *Toxicology Letters* **2015**, *232* (2), 519–
904 532. https://doi.org/10.1016/j.toxlet.2014.09.004.

905 (56) Banerjee, P.; Eckert, A. O.; Schrey, A. K.; Preissner, R. ProTox-II: A Webserver for the
906 Prediction of Toxicity of Chemicals. *Nucleic Acids Res* **2018**, *46* (W1), W257–W263.
907 https://doi.org/10.1093/nar/gky318.
908 (57) Mansouri, K.; Kleinstreuer, N.; Abdelaziz, A. M.; Alberga, D.; Alves, V. M.; Andersson, P.
909 L.; Andrade, C. H.; Bai, F.; Balabin, I.; Ballabio, D.; Benfenati, E.; Bhhatarai, B.; Boyer,
910 S.; Chen, J.; Consonni, V.; Farag, S.; Fourches, D.; García-Sosa, A. T.; Gramatica, P.;
911 Grisoni, F.; Grulke, C. M.; Hong, H.; Horvath, D.; Hu, X.; Huang, R.; Jeliazkova, N.; Li, J.;
912 Li, X.; Liu, H.; Manganelli, S.; Mangiatordi, G. F.; Maran, U.; Marcou, G.; Martin, T.;
913 Muratov, E.; Nguyen, D.-T.; Nicolotti, O.; Nikolov, N. G.; Norinder, U.; Papa, E.;
914 Petitjean, M.; Piir, G.; Pogodin, P.; Poroikov, V.; Qiao, X.; Richard, A. M.; Roncaglioni,
915 A.; Ruiz, P.; Rupakheti, C.; Sakkiah, S.; Sangion, A.; Schramm, K.-W.; Selvaraj, C.; Shah,
916 I.; Sild, S.; Sun, L.; Taboureau, O.; Tang, Y.; Tetko, I. V.; Todeschini, R.; Tong, W.;
917 Trisciuzzi, D.; Tropsha, A.; Van Den Driessche, G.; Varnek, A.; Wang, Z.; Wedebye, E. B.;
918 Williams, A. J.; Xie, H.; Zakharov, A. V.; Zheng, Z.; Judson, R. S. CoMPARA:
919 Collaborative Modeling Project for Androgen Receptor Activity. *Environ Health Perspect*
920 **2020**, 128 (2). https://doi.org/10.1289/EHP5580.

921 (58) Mansouri, K.; Abdelaziz, A.; Rybacka, A.; Roncaglioni, A.; Tropsha, A.; Varnek, A.;
922 Zakharov, A.; Worth, A.; Richard, A. M.; Grulke, C. M.; Trisciuzzi, D.; Fourches, D.;
923 Horvath, D.; Benfenati, E.; Muratov, E.; Wedebye, E. B.; Grisoni, F.; Mangiatordi, G. F.;
924 Incisivo, G. M.; Hong, H.; Ng, H. W.; Tetko, I. V.; Balabin, I.; Kancherla, J.; Shen, J.;
925 Burton, J.; Nicklaus, M.; Cassotti, M.; Nikolov, N. G.; Nicolotti, O.; Andersson, P. L.;
926 Zang, Q.; Politi, R.; Beger, R. D.; Todeschini, R.; Huang, R.; Farag, S.; Rosenberg, S. A.;
927 Slavov, S.; Hu, X.; Judson, R. S. CERAPP: Collaborative Estrogen Receptor Activity
928 Prediction Project. *Environ Health Perspect* **2016**, *124* (7), 1023–1033.
929 https://doi.org/10.1289/ehp.1510267.
930 (59) Mao, J.; Akhtar, J.; Zhang, X.; Sun, L.; Guan, S.; Li, X.; Chen, G.; Liu, J.; Jeon, H.-N.;

931 Kim, M. S.; No, K. T.; Wang, G. Comprehensive Strategies of Machine-Learning-Based

932 Quantitative Structure-Activity Relationship Models. *iScience* **2021**, *24* (9), 103052.

933 https://doi.org/10.1016/j.isci.2021.103052.

934

935