Predictive Models for Compound Binding to Androgen and Estrogen Receptors Based on Counter-Propagation Artificial Neural Networks

Endocrine-disrupting chemicals (EDCs) are exogenous substances that interfere with the normal function of the human endocrine system. These chemicals can affect specific nuclear receptors, such as androgen receptors (ARs) or estrogen receptors (ER) α and β, which play a crucial role in regulating complex physiological processes in humans. It is now more crucial than ever to identify EDCs and reduce exposure to them. For screening and prioritizing chemicals for further experimentation, the use of artificial neural networks (ANN), which allow the modeling of complicated, nonlinear relationships, is most appropriate. We developed six models that predict the binding of a compound to ARs, ERα, or ERβ as agonists or antagonists, using counter-propagation artificial neural networks (CPANN). Models were trained on a dataset of structurally diverse compounds, and activity data were obtained from the CompTox Chemicals Dashboard. Leave-one-out (LOO) tests were performed to validate the models. The results showed that the models had excellent performance with prediction accuracy ranging from 94% to 100%. Therefore, the models can predict the binding affinity of an unknown compound to the selected nuclear receptor based solely on its chemical structure. As such, they represent important alternatives for the safety prioritization of chemicals.


Introduction
Endocrine-disrupting chemicals (EDCs) are exogenous chemicals that interfere with the normal function of the human endocrine system through a variety of mechanisms and, as such, represent a global human health concern [1,2]. Certain chemicals can act through specific nuclear receptors, such as estrogen receptors (ERs) or androgen receptors (ARs), which are involved in the regulation of many complex physiological processes in humans [3]. All of the nuclear receptors mentioned (ARs, ERs) have an amino terminal group responsible for a ligand-independent activation of transcription, a central DNAbinding domain, and a ligand-binding domain at the carboxyl end. After binding to the ligand, the receptors translocate to the nucleus and bind to specific transcription elements to eventually trigger the transcription of the target genes [4]. It should also be noted that ARs and ERs can also act through very rapid non-genomic biological responses by triggering signaling circuits initiated outside the nucleus. This rapid action occurs through the interaction of nuclear receptors with various signal or scaffold molecules and has been shown to be involved in cell cycle control, proliferation, migration, and exclusion of steroid receptors from the cell nucleus [5][6][7].
Several physiological processes are regulated by estrogens through two estrogen receptors, ERα and ERβ. ERα is mainly expressed in the uterus, ovaries, breast, liver, kidney, bone, and white adipose tissue, whereas ERβ expression is found in the ovaries, male reproductive organs, prostate, central nervous system, cardiovascular system, lung,

Data
For modeling, we created a dataset of structurally heterogeneous compounds from published data, including compound name, CAS, QSAR-Ready SMILES (SMILES representations of desalted, de-isotoped, stereo-neutral forms of chemical structures associated with specific chemical substances), and activity details (available in the Supplementary File supdata.xlsx).
Activity data were taken from the CompTox Chemicals Dashboard (https://comptox. epa.gov/dashboard/, accessed 10 October 2022) [26]. This is a publicly available, web- based application developed by the U.S. Environmental Protection Agency to provide access to systematically compiled and consolidated chemistry, toxicity, and exposure information for more than 900,000 chemicals. All Tox21 data used are assays redouts obtained by measuring reporter genes via receptor activity and developed using an inducible reporter (β-lactamase induction) detected using the GAL4-β-lactamase reporter gene [27]. The shift in fluorescence emission (from green to blue color) is used to identify the activation of the reporter gene as a result of ligand GR binding to the glucocorticoid response element. The β-lactamase reporter gene assay has two readouts: background (Channel 1, ch1, or green channel) and gene expression (Channel 2, ch2, or blue channel), which are used to calculate a ratio (ch2/ch1) for analysis. Details of assay protocols are described in EPA publication Toxicity Forecaster (Toxcast) In Vitro Assays [28]. Potentially active compounds (positive hit calls) were then determined by curve fitting using the ToxCast Data Pipeline [29].
Due to the large differences between the results of different laboratories, only compounds found to be positive in at least 50% of the laboratories. A similar number of compounds, identified as negative in the CompTox Chemicals Dashboard [26], was chosen. The compounds that were selected were those that were tested as active the least in individual laboratories.

Methods
The DRAGON [30] program package was used to calculate 3690 structural descriptors for all compounds. The descriptors that had negligible variance or were not calculated for all compounds were removed and then analyzed using the principal component analysis (PCA) method. PCA is a popular method for analyzing multidimensional data [24,25,31,32] and is a mathematical transformation of the original variables into new variables or principal components (PCs). The new variables are sorted according to the proportion of variance they contain. Typically, only a handful of new variables is necessary to describe the entire variance hidden in the data. In the reported case, thousands of descriptors were replaced with 22 new variables.
The CPANN method was used for modeling, which is an example of a self-organizing map technique developed for analyzing data in multidimensional space. This technique is based on a nonlinear projection from multidimensional space onto a network of neurons organized as a two-dimensional map. To achieve a topology-preserving projection, a nonlinear algorithm called training is used. A neuron is in fact a vector with a dimension equal to the number of descriptors. The weights are the individual components of the vectors determined by a nonlinear iterative training algorithm. The objects are presented to the networks in this procedure, and the weights are modified to be 'similar' to the descriptor values. The training is repeated until the changes in weights between two successive steps fall below a certain threshold. Since the basic property of the trained network is that similar objects are close to each other, it is expected that chemicals with similar structural features form clusters [23]. An additional output layer is added to the CPANN to the network associated with the property. In our case, the output layer is two-dimensional, with one dimension corresponding to binders and the other to non-binders. When a compound is presented to the model, it is placed in a group of 'similar' objects. The output layer is used to read the property. In our case, the prediction was given as a two-dimensional vector with elements expressing the specific class. The threshold for classifying a substance as a binder or non-binder was set to 0.5. After the leave-one-out (LOO) tests, the technical parameter of the model, namely the dimension of CPANN and number of epochs, was determined. All models were finally trained with 500 epochs.
In our study, we used four common evaluation metrics to assess and compare the performance of the model, namely specificity (SP) (Equation (1)), sensitivity (SE) (Equation (2)), accuracy (ACC) (Equation (3)), and the Matthew correlation coefficient (MCC) (Equation (4)). Calculations were based on true positive (TP), true negative (TN), false positive (FP), and false negative (FN) predictions. These were derived from the statistics of the model's prediction results, i.e., the confusion matrix.

Results and Discussion
In order to create the most accurate model possible, the careful selection of experimental data is paramount. A comparative analysis of ED-relevant data from several sources, including the CERAPP [33], CoMPARA [34,35], and Tox21 projects, as well as the ChEMBL and PubChem databases, revealed poor agreement between experimental values from different sources [36]. For this reason, we used data from only one source, the Tox21 screening program, where the results were obtained using the same assay and under the same conditions. Unfortunately, even in this case, the reproducibility of the assay is low. The determined assay reproducibility is 80.6% for the AR agonist, 67.2% for the AR antagonist, 76.9% for the Erα agonist, and 76.8% for the Erα antagonist [37]. No reproducibility was found for the ERβ assays. On the other hand, the comparison of different in silico models has shown that the size of the dataset is not directly related to the performance of the model [36]. Therefore, to ensure the highest quality of input data, we only included compounds that were identified as positive by at least 50% of the laboratories in the model.
For all selected compounds, 3690 structural descriptors were calculated using the program package DRAGON [30]. The descriptors that had negligible variance or were not calculated for all compounds were removed. In further PCA analysis, 2386, 2495, 2444, 2518, 2404, and 2484 descriptors were used for the AR agonist, AR antagonist, ERα agonist, ERα antagonist, ERβ agonist, and ERβ antagonist a, respectively. As an example, we present the scree plots ( Figure 1) for AR agonists showing that the first PC carries more than 40% of the total variance, the second PC about 8.0-9.5% of the variance and the others even less. In fact, the scree plots are very similar for all cases studied. The common feature is that the first 22 PCs account for about 80% of the total variance.    4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 % Variance explained     4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 % Variance explained  Subsequently, we considered the 22 PCs as input descriptors in the CPANN modeling. The architecture of the CPANN model and the LOO test are described in the above section. The LOO test results are shown in Table 2, considering different dimensions of CPANN.   Subsequently, we considered the 22 PCs as input descriptors in the CPANN modeling. The architecture of the CPANN model and the LOO test are described in the above section. The LOO test results are shown in Table 2, considering different dimensions of CPANN. As a final model, we chose the 18 × 18 dimension for the models, except for the Erβagonists 8 × 8 due to the small number of molecules. In the case of the Erβ-agonists, the model gives perfect separation, which is a consequence of the small number of molecules that have similar structures, especially for non-active molecules. The upper maps are shown in Figure 3.   Light blue represents neuron with one inactive compound, dark blue neuron with two or more inactive compounds, light red represents neuron with one active compound, dark red represent neuron with two or more active compounds, and green represents conflict neuron, where at least one substance is classified differently than others.
The model for the AR agonist receptor shows near-perfect separation with only one misclassified molecule and an accuracy of 0.994. The only conflict (highlighted in green) in the AR agonist 'recall-ability' test presented in the top map is on the neuron (118) occupied by deslorelin and eledoisin, both short peptides with a molecular weight of >1000 g/mol and c-terminal pyrimidine, which are difficult to separate in the model. The top map of the AR antagonist shows nine conflict neurons. On neuron (1,6), there are four compounds (romidepsin, dronedarone, nelfinavir mesylate, and clindamycin palmitate) with molecular weights between 540 and 700 g/mol, which contain sulfur. Only clindamycin palmitate was non-binding in in vitro experiments, while the other three compounds were found to be binding. On the other hand, there are also four molecules on neuron (18,  4 5 6 7 8 9 10 11 12 13 14 15 16 17 4 5 6 7 8 9 10 11 12 13 14 15 16 17   Light blue represents neuron with one inactive compound, dark blue neuron with two or more inactive compounds, light red represents neuron with one active compound, dark red represent neuron with two or more active compounds, and green represents conflict neuron, where at least one substance is classified differently than others. The model for the AR agonist receptor shows near-perfect separation with only one misclassified molecule and an accuracy of 0.994. The only conflict (highlighted in green) in the AR agonist 'recall-ability' test presented in the top map is on the neuron (118) occupied by deslorelin and eledoisin, both short peptides with a molecular weight of >1000 g/mol and c-terminal pyrimidine, which are difficult to separate in the model. The top map of the AR antagonist shows nine conflict neurons. On neuron (1,6), there are four compounds (romidepsin, dronedarone, nelfinavir mesylate, and clindamycin palmitate) with molecular weights between 540 and 700 g/mol, which contain sulfur. Only clindamycin palmitate was non-binding in in vitro experiments, while the other three compounds were found to be binding. On the other hand, there are also four molecules on neuron (18,9) presented in Table 3, all of which contain a heterocycle with nitrogen (1-butyl-4-methylpyridinium, 1-butylpyridinium, 1,8-diazabicyclo [5.4.0]undec-7-ene, and butylmethyl-imidazolium), with 1-butyl-4-methylpyridinium binding as an antagonist to AR in in vitro experiments, whereas the other three compounds were not binding. We would like to point out that 1-butyl-4-methylpyridinium and 1-butylpyridinium differ only by a methyl group attached to the pyridine, while the in vitro experiment shows a different classification. Table 3. Compounds on neuron (18,9) in top map for AR antagonist.

Structure AR Antagonist Classification
In Vitro In Silico

AR Antagonist Classification In Vitro
In Silico Similarly, the compounds 1-butyl-4-methylpyridinium and 1-butyl-2-methylpyridinium at the conflict neuron (12,6) in the ERα-antagonist model are structural isomers (Table 4). They have the same chemical formula but differ in the position of the methyl group. They were classified as non-binding and binding, respectively, in vitro, whereas our model classified both as non-binding.  Table 3, all of which contain a heterocycle with nitrogen (1-butyl-4methylpyridinium, 1-butylpyridinium, 1,8-diazabicyclo [5.4.0]undec-7-ene, and butylmethyl-imidazolium), with 1-butyl-4-methylpyridinium binding as an antagonist to AR in in vitro experiments, whereas the other three compounds were not binding. We would like to point out that 1-butyl-4-methylpyridinium and 1-butylpyridinium differ only by a methyl group attached to the pyridine, while the in vitro experiment shows a different classification. Table 3. Compounds on neuron (18,9) in top map for AR antagonist.

AR Antagonist Classification In Vitro
In Silico Similarly, the compounds 1-butyl-4-methylpyridinium and 1-butyl-2-methylpyridinium at the conflict neuron (12,6) in the ERα-antagonist model are structural isomers (Table 4). They have the same chemical formula but differ in the position of the methyl group. They were classified as non-binding and binding, respectively, in vitro, whereas our model classified both as non-binding. non-binding non-binding 1-Butyl-3methylimidazolium 9) presented in Table 3, all of which contain a heterocycle with nitrogen (1-butyl-4methylpyridinium, 1-butylpyridinium, 1,8-diazabicyclo [5.4.0]undec-7-ene, and butylmethyl-imidazolium), with 1-butyl-4-methylpyridinium binding as an antagonist to AR in in vitro experiments, whereas the other three compounds were not binding. We would like to point out that 1-butyl-4-methylpyridinium and 1-butylpyridinium differ only by a methyl group attached to the pyridine, while the in vitro experiment shows a different classification. Table 3. Compounds on neuron (18,9) in top map for AR antagonist.

AR Antagonist Classification In Vitro
In Silico Similarly, the compounds 1-butyl-4-methylpyridinium and 1-butyl-2-methylpyridinium at the conflict neuron (12,6) in the ERα-antagonist model are structural isomers (Table 4). They have the same chemical formula but differ in the position of the methyl group. They were classified as non-binding and binding, respectively, in vitro, whereas our model classified both as non-binding.   (18,9) in top map for AR antagonist.

AR Antagonist Classification In Vitro
In Silico Similarly, the compounds 1-butyl-4-methylpyridinium and 1-butyl-2-methylpyridinium at the conflict neuron (12,6) in the ERα-antagonist model are structural isomers (Table 4). They have the same chemical formula but differ in the position of the methyl group. They were classified as non-binding and binding, respectively, in vitro, whereas our model classified both as non-binding.

non-binding non-binding
Similarly, the compounds 1-butyl-4-methylpyridinium and 1-butyl-2-methylpyridinium at the conflict neuron (12,6) in the ERα-antagonist model are structural isomers (Table 4). They have the same chemical formula but differ in the position of the methyl group. They were classified as non-binding and binding, respectively, in vitro, whereas our model classified both as non-binding.
In the ERβ-antagonist model, three structurally related compounds (fenticlor, bithionol, and phencapton) are found in the conflict neuron (1, 1) ( Table 5). Although fenticlor was non-binding in the in vivo experiment, the other two were. All three molecules contain chlorine bound to an aromatic ring. However, fenticlor and bithionolate show greater similarity in structure, as both consist of two chlorinated phenols bridged with sulfur to form C1, differing only in the number of chlorines attached to the phenol. These examples clearly demonstrate the robustness of the model, as it only has difficulty classifying molecules that are structurally closely related.  In the ERβ-antagonist model, three structurally related compounds (fenticlor, bithionol, and phencapton) are found in the conflict neuron (1,1) ( Table 5). Although fenticlor was non-binding in the in vivo experiment, the other two were. All three molecules contain chlorine bound to an aromatic ring. However, fenticlor and bithionolate show greater similarity in structure, as both consist of two chlorinated phenols bridged with sulfur to form C1, differing only in the number of chlorines attached to the phenol. These examples clearly demonstrate the robustness of the model, as it only has difficulty classifying molecules that are structurally closely related.  In the ERβ-antagonist model, three structurally related compounds (fenticlor, bithionol, and phencapton) are found in the conflict neuron (1,1) ( Table 5). Although fenticlor was non-binding in the in vivo experiment, the other two were. All three molecules contain chlorine bound to an aromatic ring. However, fenticlor and bithionolate show greater similarity in structure, as both consist of two chlorinated phenols bridged with sulfur to form C1, differing only in the number of chlorines attached to the phenol. These examples clearly demonstrate the robustness of the model, as it only has difficulty classifying molecules that are structurally closely related.  In the ERβ-antagonist model, three structurally related compounds (fenticlor, bithionol, and phencapton) are found in the conflict neuron (1,1) ( Table 5). Although fenticlor was non-binding in the in vivo experiment, the other two were. All three molecules contain chlorine bound to an aromatic ring. However, fenticlor and bithionolate show greater similarity in structure, as both consist of two chlorinated phenols bridged with sulfur to form C1, differing only in the number of chlorines attached to the phenol. These examples clearly demonstrate the robustness of the model, as it only has difficulty classifying molecules that are structurally closely related. In the ERβ-antagonist model, three structurally related compounds (fenticlor, bithionol, and phencapton) are found in the conflict neuron (1,1) ( Table 5). Although fenticlor was non-binding in the in vivo experiment, the other two were. All three molecules contain chlorine bound to an aromatic ring. However, fenticlor and bithionolate show greater similarity in structure, as both consist of two chlorinated phenols bridged with sulfur to form C1, differing only in the number of chlorines attached to the phenol. These examples clearly demonstrate the robustness of the model, as it only has difficulty classifying molecules that are structurally closely related. Unfortunately, ED-relevant data from several sources showed poor agreement between experimental values from different sources [26]. Therefore, we were not able to create a test set for the external validation of the models. However, the LOO validation procedure performed results in a reliable and unbiased estimate of model performance.
Following a previous study in which a comparative analysis of docking models (Endocrine disruptome software and VirtualToxLab) and their interpretation using Tox21 in in vitro data was performed, our models show better accuracy than the selected docking models. The prediction accuracy of the endocrine disruptome model was 0.53 for AR, 0.58 for Erα, and 0.62 for Erβ, while Virtual ToxLab showed slightly better results of 0.57, 0.59, and 0.62, respectively [38]. In this regard, the QSAR models show better accuracy. The comparison of different published QSAR models for screening AR showed that the balanced accuracy ranged from 0.55 to 0.81 [26]. The QSAR model with the best performance was presented by Todorov et al. [39], and was based on the multiparameter formulation binding binding Unfortunately, ED-relevant data from several sources showed poor agreement between experimental values from different sources [26]. Therefore, we were not able to create a test set for the external validation of the models. However, the LOO validation procedure performed results in a reliable and unbiased estimate of model performance.
Following a previous study in which a comparative analysis of docking models (Endocrine disruptome software and VirtualToxLab) and their interpretation using Tox21 in in vitro data was performed, our models show better accuracy than the selected docking models. The prediction accuracy of the endocrine disruptome model was 0.53 for AR, 0.58 for Erα, and 0.62 for Erβ, while Virtual ToxLab showed slightly better results of 0.57, 0.59, and 0.62, respectively [38]. In this regard, the QSAR models show better accuracy. The comparison of different published QSAR models for screening AR showed that the balanced accuracy ranged from 0.55 to 0.81 [26]. The QSAR model with the best performance was presented by Todorov et al. [39], and was based on the multiparameter formulation of the common reactivity pattern (COREPA) approach. Better results were obtained with OASIS, a QSAR platform covering both ER and AR binding [40]. Evaluation with ToxCast™ of in vitro binding data gave an accuracy of ER of 0.88 and of AR of 0.84 [41]. Despite the fact that the above QSAR models show excellent results in compound classification, we managed to further improve the accuracy (0.94-1) of the predictions (presented in Table 2) with CPANN. Moreover, with the exception of the endocrine disruptome, none of the aforementioned models distinguishes between agonists and antagonists, which is an additional added value of our models.

Conclusions
In the interest of avoiding animal testing, in silico modeling focuses on models for categorizing and grouping substances using various chemometric approaches. In this work, we present CPANN models for agonistic and antagonistic binding to the androgen and estrogen (ERα and ERβ) receptors. The structural properties of the compounds were fully described by the set of molecular structural descriptors given by DRAGON. Variable reduction was performed by the PCA-based clustering technique. In all cases, the first PC contributes about 40% to the total variance, and thus, plays a crucial role in the distribution of the objects. Score plots using the first two PCs show a weak separation between active and non-active compounds. Subsequently, we used CPANN models for classification. Notwithstanding the uncertainties in the classification of the compounds in the selected in vitro data, we divided the compounds into binders and non-binders. The models that were created were validated using statistical parameters for LOO, which showed that the models had a high level of balance and accuracy. Comparison with other QSAR and docking models presented in the literature with the CPANN models suggests that CPANN models are superior, although an absolute measure of comparison is not possible due to the different sets of compounds used. As far as we know, these are only QSAR models that distinguish between antagonists and agonists. This is an important detail that further increases their usefulness. Overall, the CPANN models we created are robust and predictive of new chemicals. Based on chemical structure alone, the models can predict the binding affinity of an unknown compound to the chosen nuclear receptor. Thus, these models provide a viable, effective, and practical tool for rapidly screening the endocrine activity of organic compounds and prioritizing them for further testing.

Data Availability Statement:
In vitro activity data were taken from the CompTox Chemicals Dashboard: https://comptox.epa.gov/dashboard/ (accessed on 10 October 2022). Further data supporting the findings of this study are available from the corresponding authors upon reasonable request.