CERAPP: Collaborative Estrogen Receptor Activity Prediction Project

1There are tens of thousands of natural and synthetic chemical substances to which humans and wildlife are exposed. A subset of these compounds may disrupt normal functioning of the endocrine system and cause health hazards to both humans and ecological species. Endocrine-disrupting chemicals (EDCs) can mimic or interfere with natural hormones and alter their mechanisms of action at the receptor level, as well as interfere with the synthesis, transport, and metabolism of endogenous hormones. Exposure to EDCs can lead to adverse health effects involving developmental, neurological, reproductive, metabolic, cardiovascular, and immune systems in humans and wildlife.

2The estrogen receptor (ER) is one of the most extensively studied targets related to the effects of EDCs. This concern about estrogen-like activity of man-made chemicals is because of their potential for negatively affecting reproductive function. The emergence of concerns about EDCs has resulted in regulations requiring assessment of chemicals for estrogenic activity. There are numerous in vitro and in vivo protocols to identify potential endocrine pathway-mediated effects of chemicals, including interactions with hormone receptors. However, experimental testing of chemicals is expensive and time-consuming and currently impractical for application to the vast number of synthetic chemicals in use. Consequently, toxicological data and especially estrogenic activity data are available only for a limited number of compounds.

3The use of in silico approaches, such as quantitative structure–activity relationships (QSARs), is an alternative to bridge the lack of knowledge about chemicals when little or no experimental data are available. These structure-based methods are particularly appealing for their ability to predict toxicologically relevant end points quickly and at low cost. QSARs have been promoted and their use recognized since the pioneering work of Hansch in the 1960s. The conceptual basis of QSARs is that chemicals with similar structures are hypothesized to exhibit similar behavior in living organisms. Thus, it should be possible to predict biological activity of new chemicals based on published experimental data. Several guidance documents to develop these modeling techniques are available in the literature.

4Recently, in vitro high-throughput screening (HTS) assays have emerged and become a viable tool for large-scale chemical testing. HTS generates substantial amounts of data that can be used as a knowledge base to correlate chemical structures to their biological activities. Thus, QSARs can identify key structural characteristics in active chemicals and can use them to virtually screen large chemical libraries. Although there is concern about the overall accuracy of a QSAR model to predict the true activity of a particular chemical, accuracy can be high enough to use the results for prioritizing chemicals that are worth subjecting to experimental testing.

5With the increasing number of new substances submitted to the U.S. EPA and the European Chemicals Agency for registration (~ 1,500 chemicals every year), there is a need to prioritize chemicals to speed up the process and lower the overall costs of testing. The Toxicology Testing in the 21st Century (Tox21) collaboration and the U.S. EPA’s Toxicity ForeCaster (ToxCast™)projects are screening thousands of chemicals in HTS in vitro assays for a broad range of targets. Relevant to this paper, these two projects have in common ~ 1,800 chemicals tested in a battery of 18 ER-related assays.

6This paper describes the results of the Collaborative Estrogen Receptor Activity Prediction Project (CERAPP), which was organized by the National Center for Computational Toxicology at the U.S. EPA. The aim of the project was to use ToxCast™/Tox21 ER HTS assay data to develop and optimize predictive computational models, and to use their predictions to prioritize a large chemical universe of 32,464 unique chemical structures for further testing. Seventeen research groups from the United States and Europe participated in this project. These groups submitted 40 categorical models and 8 continuous models using different QSAR and structure-based approaches. Most of the newly developed models used a training set consisting of 1,677 chemicals, each assigned a potency score quantifying their ER agonist, antagonist, and binding activities, obtained from a computational network model that integrates data from 18 diverse ER HTS assays. All models were evaluated and weighted based on their prediction accuracy scores (including sensitivity and specificity) using ToxCast™/Tox21 HTS data, as well as an evaluation data set collected from different literature sources. To overcome the limitations of single models, all predictions were combined into a consensus model that classified the chemicals into active/inactive binders, agonists, and antagonists and provided estimates of their potency level relative to known reference chemicals.

7The 17 international research groups that participated in this project are listed in alphabetic order in Table S1. The goals of the project, outlined in Table S2, were achieved in multiple steps, including chemical structure curation, experimental data preparation from the literature, modeling and prediction, model evaluation, consensus strategy development, and consensus modeling. Each step was assigned to a subgroup of participants according to their interests and areas of expertise.

8Provided training set. The data that were suggested to be used by the participants as a training set to develop and optimize their models was derived from ToxCast™ and Tox21 programs. Concentration-response data from a collection of 18 in vitro HTS assays exploring multiple sites in the mammalian ER pathway were generated for 1,812 chemicals. This chemical library included 45 reference ER agonists and antagonists (including negatives), as well as a wide array of commercial chemicals with known estrogen-like activity. A mathematical model was developed to integrate the in vitro data and calculate an area under the curve (AUC) score, ranging from 0 to 1, which is roughly proportional to the consensus AC50 value across the active assays. A given chemical was considered active if its agonist or antagonist score was higher than 0.01. In order to reduce the number of potential false positives this threshold can be increased to 0.1.

9Prediction set. We identified > 50,000 chemicals [at the level of Chemical Abstracts Service Registry Number (CASRN)] for use in this project as a virtual screening library to be prioritized for further testing and regulatory purposes. This set was intended to include a large fraction of all man-made chemicals to which humans may be exposed. These chemicals were collected from different sources with significant overlap and cover a variety of classes, including consumer products, food additives, and human and veterinary drugs. The following list includes the sources used in this project:

10Chemicals with documented use, and therefore, with exposure potential (~ 43,000). Available in the U.S. EPA chemical product categories database (CPCat), which is part of the Aggregated Computational Toxicology Resource (ACToR) system.

11The Distributed Structure-Searchable Toxicity (DSSTox). A list of ~ 15,000 curated chemical structures from multiple inventories of environmental interest. In particular, structures for all of the ToxCast™ and Tox21 chemicals are included.

12The Canadian Domestic Substances list (DSL). A compiled list of all substances thought to be in commercial use in Canada (~ 24,000 chemicals). Thus, it includes chemicals with potential human or ecological exposure.

13The Endocrine Disruption Screening Program (EDSP) universe of ~ 10,000 chemicals. The U.S. EPA’s EDSP is required to test certain chemicals for their potential for endocrine disruption.

14A list of ~ 15,000 chemicals used as training and test sets for the different models implemented in the U.S. EPA’s Estimation Program Interface (EPI Suite™) to predict physico-chemical properties.

15This virtual chemical library has undergone stringent chemical structure processing and normalization for use in the QSAR modeling study (see “Chemical Structure Curation”) and made available for download on ToxCast™ Data web site under CERAPP data (https://www3.epa.gov/research/COMPTOX/CERAPP_files.html, PredictionSet.zip), is intended to be employed for a large number of other QSAR modeling projects, not just those focused on endocrine-related targets.

16Experimental evaluation set. A large volume of estrogen-related experimental data has accumulated in the literature over the past two decades. The information on the estrogenic activity of chemicals was mined and curated to serve as a validation set for predictions of the different models. For this purpose, in vitro experimental data were collected from different overlapping sources, including the U.S. EPA’s HTS assays, online databases, and other data sets used by participants to train models:

17HTS data from Tox21 project consisting of ~ 8,000 chemicals evaluated in four assays, extending beyond the 1,677 used in the training set.

18The U.S. FDA Estrogenic Activity Database (EADB), which consists of literature derived ER data for ~ 8,000 chemicals.

19Estrogenic data for ~ 2,000 chemicals from the METI (Ministry of Economy, Trade and Industry, Japan) database.

20Estrogenic data for ~ 2,000 chemicals from ChEMBL database.

21The full data set consisted of > 60,000 entries, including binding, agonist, and antagonist information for ~ 15,000 unique chemical structures. For the purpose of this project, this data set was cleaned and made more consistent by removing in vivo data, cytotoxicity information, and all ambiguous entries (missing values, undefined/nonstandard end points, and unclear units). Only 7,547 chemical structures from the experimental evaluation set that overlapped with the CERAPP prediction set, for a total of 44,641 entries, were kept and made available for download on the U.S. EPA ToxCast™ Data web site (https://www3.epa.gov/research/COMPTOX/CERAPP_files.html, EvaluationSet.zip). The non-CERAPP chemicals were excluded from the evaluation set (see “Chemical Structure Curation” section). Then, all data entries were categorized into three assay classes: (a) binding, (b) reporter gene/transactivation, or (c) cell proliferation. The training set end point to model is the ER model AUC that parallels the corresponding individual assay AC50 values, and therefore all units for activities in the experimental data set were converted to μM to have approximately equivalent concentration–response values for the evaluation set. Chemicals with cell proliferation assays were considered as actives if they exceeded an arbitrary threshold of 125% proliferation. For entries where testing concentrations were reported in the assay name field, those values were converted to μM and considered as the AC50 value if the compound was reported as active. All inactive compounds were arbitrarily assigned an AC50 value of 1 M.

22Chemical structures collected from different public sources contained many duplicates, and inconsistencies in the molecular structures. Hence, a structure curation process was carried out to derive a unique set of QSAR-ready structures. All participating groups then used this consistent set of structures for both training and prediction steps. It should be noted that each group likely employed different descriptor calculation software, which could effectively alter structures in some cases. Several different curation approaches were combined into a unique procedure used for this project. The free and open-source data-mining environment KNIME (Konstanz Information Miner) was selected to design a curation workflow to process all structures and provide consistent training and prediction sets. The workflow performed a series of curation steps:

23The original files containing structures in different formats were parsed, checked for valences, and for the integrity of the required structural information to render the molecules. Invalid entries were corrected by retrieving a new structure from online databases using web services [PubChem, ChemSpider] or removed if ambiguous.

24The first filter was applied to check for the presence of carbon atoms and remove inorganic compounds.

25The structures were desalted, and inorganic counterions were removed.

26The second filter, based on molecular weight, was applied and chemicals exceeding a threshold of 1,000 g/mol were removed to speed up molecular descriptor calculations and model calibration.

27Valid QSAR modeling practice requires all chemicals to be structurally consistent by converting tautomers to unique representations. Thus, a series of transformations was applied on the structures to standardize nitro and azide mesomers, keto-enol tautomers, enamine-imine tautomers, ynol-ketene, and other conversions.

28These transformations were followed by neutralizing the charged structures, when possible, and removing the stereochemistry information.

29Explicit hydrogen atoms were added, and structures were aromatized according to Hückel’s rules implemented in KNIME.

30The duplicates were removed using the IUPAC (International Union of Pure and Applied Chemistry) InChI (International Chemical Identifier) codes because these are unequivocal identifiers.

31The final filter was applied to remove chemicals containing metals that often cause problems in molecular descriptor calculations.

32Both training and prediction sets were processed by the same structure curation workflow. At the end of this procedure, 32,464 unique structures—the 32 K set—remained in the prediction set and 1,677 in the training set. These two data sets are made available for download in structure data file (SDF) format on the U.S. EPA ToxCast™ Data web site (https://www3.epa.gov/research/COMPTOX/CERAPP_files.html, TrainingSet.zip and PredictionSet.zip). The identity of these chemicals (name, CASRN) was not provided to the participating modeling groups during the modeling process.

33The participant groups adopted different approaches and used several software programs (proprietary or open-source [commercial or free]) to calibrate categorical and continuous models to the training data (Table 1). A categorical model is one that provides an active/inactive call for each chemical, whereas a continuous model provides a prediction of the potency (in μM) for each active chemical. Models were developed using both well-known and innovative methods including partial least-squares (PLS), partial least-squares discriminant analysis (PLS-DA), decision forest (DF), three-dimensional (3D) quantitative spectral data–activity relationship (QSDAR), support vector machines (SVM), k nearest neighbors (kNN), associative artificial neural networks (ASNN), PASS algorithm derived from Naïve Bayes classifier, self-consistent regression with radial basis function interpolation (RBF-SCR), OCHEM machine learning methods, docking and consensus of different approaches. The set of 1,677 chemicals provided by the U.S. EPA was used by more than 90% of the participating groups as a training set to fit their models, but some pre-existing models were also used that had been trained using other data sets from the literature such as. In addition, each group performed its own analysis to select the appropriate chemicals to be considered as a training set according to their particular modeling procedure. For descriptor calculation and docking procedures, some of the programs used were LeadScope, PaDEL-Descriptor, QikProp (version 3.4, http://www.schrodinger.com/QikProp/), multilevel and quantitative neighborhoods of atoms (MNA, QNA) used by GUSAR and PASS, DRAGON, Mold2, GLIDE (version 6.5, http://www.schrodinger.com/Glide), AutoDock, ISIDA, and other fingerprint generators. Some of the participants applied feature selection techniques, such as genetic algorithms (GAs) and random forest (RF). These techniques were applied after calculating descriptors to reduce collinearity and variable dimensionality to keep only the most informative descriptors in the models.

34All molecular structures of chemicals collected for the evaluation set from the different sources were curated and standardized using the previously described KNIME workflow (Table S2, step 2). All data used as the evaluation set for categorical and continuous models are available on the U.S. EPA ToxCast™ web site (https://www3.epa.gov/research/COMPTOX/CERAPP_files.html, EvaluationSet.zip).

35Standard InChI codes were generated in KNIME and used to identify the chemicals. Data-mining tools available in the KNIME environment were used to concatenate and unify the different information fields from the different sources (CASRN, chemical name, original structure, standardized structure, InChI code, assay name, assay class, protein subtype, species, end point name, end point value, end point unit, and literature reference). Although ToxCast™ chemicals were used in the training sets of many models, they were not removed from the evaluation set to investigate how the predictions will perform on the literature data because there are differences between the AUC values and the literature data and because the sources from which the evaluation set was collected were not fully verified (we cannot assume that all cytotoxicity information was already fully cleaned).

36Evaluation set for categorical models. An important issue with the literature-derived evaluation set was the inconsistency of the results from different sources. To minimize this, the available entries for each chemical structure were grouped into binders, agonists, and antagonists. The results were then categorized into active and inactive classes using all available literature sources by applying three rules:

37If, for a specific chemical within one of the three classes (binding, agonist, and antagonist), the disagreement among the different sources exceeded 20% (e.g., two sources indicating active agonist and three indicating inactive agonist), that chemical was removed from the evaluation data set of that specific class.

38If a chemical was an active agonist or antagonist, it also was considered as an active binder if the information was not available.

39If a chemical was an inactive agonist and inactive antagonist, it was considered also as nonbinder if the information was not available.

40This procedure resulted in a total of 7,522 unique chemical structures with activity data to be used for evaluation of the categorical models (Table 2). It is also available for download on the U.S. EPA ToxCast™ web site (https://www3.epa.gov/research/COMPTOX/CERAPP_files.html, EvaluationSet.zip).

41Evaluation set for continuous models. For active chemicals with available quantitative information from concentration-response assays, the log10-median of the literature values was calculated. Only entries with equivalent end points were considered (e.g., PC50 and EC50). This resulted in 7,253 unique chemicals with quantitative information (Table 3 and https://www3.epa.gov/research/COMPTOX/CERAPP_files.html, EvaluationSet.zip). To reduce the variability that increased with the disparate literature sources, the chemicals with quantitative information were categorized into five potency activity classes: inactive, very weak, weak, moderate, and strong. These five classes were used to evaluate the quantitative predictions. A list of 36 known active and inactive reference chemicals was used for calibrating the mapping from quantitative potency values to the activity potency classes. These same chemicals were used to validate the mathematical model used to generate the AUC values for the training set. The following thresholds were applied to the concentration–response values:

42Strong: Activity concentration below 0.09 μM.

43Moderate: Activity concentration between 0.09 and 0.18 μM.

44Weak: Activity concentration between 0.18 and 20 μM.

45Very Weak: Activity concentration between 20 and 800 μM.

46Inactive: Activity concentration higher than 800 μM.

47The five classes were assigned scores from 0 (inactive) to 1 (strong) with 0.25 increments. Then, for each chemical, the arithmetic mean of the scores of the merged entries from different literature sources was calculated. A new class was assigned to the merged entries according to the following thresholds.

48Strong: Average score > 0.75

49Moderate: 0.5 < Average score between ≤ 0.75

50Weak: 0.25 < Average score ≤ 0.5

51Very weak: 0 < Average score ≤ 0.25

52Inactive: Average score = 0

53The number of entries in each class for binding, agonist, and antagonist are summarized in Table 3.

54Evaluation procedure. This section is focused on the categorical models for their high number compared to the continuous models. The procedure used to evaluate the predictions of the participant groups was based on the categorical and continuous experimental data from ToxCast™ and the evaluation set from the literature. All continuous and categorical models for binding, agonist, and antagonist were evaluated separately on the overlap between their predicted chemicals and the following sets of chemicals (Table S3).

55Chemicals in the U.S. EPA’s ToxCast™ data set (n = 1,529 chemicals after excluding those in the ambiguous AUC range of 0.01–0.1).

56All chemicals in the full literature data (all literature sources combined).

57All chemicals with at least two literature sources.

58All chemicals from the literature data excluding the very weak actives.

59Chemicals within the applicability domain (AD) of each model (if provided).

60Chemicals remaining after applying the previous three filters in steps 3, 4, and 5 to reduce ambiguous predictions (single literature source, very weak actives, and predictions outside the AD).

61To evaluate the models on different criteria, we first determined the sensitivity (fraction of accurately predicted actives out of all actives), specificity (fraction of accurately predicted inactives out of all inactives), and balanced accuracy (BA; average of sensitivity and specificity) for each subgroup of chemicals according to each model. We then used BA values to derive two summary scores for each model, as described below.

62Score_1. Evaluation includes BA of each of the six steps weighted by the fraction of predicted chemicals of the same step, as well as the fraction of the predicted chemicals out of the full prediction set. This score (Equation 1) favors models with a wider AD and those predicting a maximum number of chemicals.

63where BA is balanced accuracy, N_pred is the number of predicted chemicals by a specific model, N_total is the total number of chemicals in the prediction set, Nfilters represents the number of five filters applied to the evaluation set chemicals and i the steps 2, 3, 4, 5, and 6.

64Score_2. Evaluation includes the BA of the model on the ToxCast™ data and the BA on the unambiguous chemicals (i.e., the subgroup of chemicals from the literature that remained after excluding chemicals with only one literature source, very weak chemicals, and chemicals outside of the AD, if provided). It favors models that focused on predicting more accurately but potentially with a narrower AD (Equation 2).

65score_2 = 1/2 (BAToxCast + BAall filters) [2]

66The quantitative predictions were evaluated as categorical models (using the BA) of the five classes after converting the numerical predictions to potency classes as defined earlier (see “Evaluation set for continuous models” section). Scores of the continuous models were calculated using Equation 2.

67The consensus predictions were generated for binders, agonists, and antagonists separately. For each chemical, we derived the average Score 2 value for all categorical models that predicted the chemical as active, and the average Score 2 value for all categorical models that predicted the chemical as inactive; we used the higher of the two averages to classify the chemical as active or inactive. Models that did not provide a prediction for the chemical in question were not included when deriving the average scores. We used Score 2 to derive the consensus classifications because its value for individual models is not penalized for the number of chemicals not predicted by the model. Also, the concordance among models on both active and inactive classes was calculated for each chemical as the fraction of models with positive and negative prediction, respectively.

68Considering only the models that provided predictions, the sum of the concordance among models for actives and inactives is equal to 1. Because most models were associated with comparable scores, the average score used to classify chemicals was mostly in agreement with model concordance (i.e., the average score for actives is high when the concordance among the models with active predictions is high and vice versa). The few exceptions were noticed when model concordance was around 0.5, which means only one or two models were driving the classification.

69For continuous predictions, the weight (w) for each chemical i was calculated from the scores (Equation 3):

70where n is the total number of models that provided predictions for the chemical i, and scorej is the score of the jth model predicting chemical i.

71Next, the consensus potency level Ci of each chemical was determined using the predicted potency classes Pj of the n available models and their corresponding weights w as follows (Equation 4):

72A total of 48 models were received from the 17 participant groups. Each group provided at least 1 categorical model for binding. Only 8 groups built models for agonists, and 6 groups built models for antagonists. The limited number of models for agonists and antagonists was the result of the low number of actives, which caused the training set to be highly unbalanced. The total number of models in each class (Table 1; see also Tables S3 and S5) was a) binding models: 21 categorical and 3 continuous, b) agonist models: 11 categorical and 3 continuous, and c) antagonist models: 8 categorical and 2 continuous.

73The participating groups provided predictions for uneven fractions of the 32 k set. AD information on model predictions was provided by only six groups. All predictions for the individual models are provided on the U.S. EPA ToxCast™ web site (https://www3.epa.gov/research/COMPTOX/CERAPP_files.html, Models.zip).

74The same evaluation procedure was applied to all models following the previously described steps. Note that some models were built using training sets other than what was provided in CERAPP and that these alternative training sets were not all publicly available. Hence, none of the training set chemicals were excluded from the evaluation sets (Table 1). Each model was evaluated on the overlap between the predicted chemicals and the two previously mentioned data sets: ToxCast™ data and the evaluation set collected from the literature. The evaluation results for categorical models are summarized in Table S3. The detailed statistics, including sensitivity and specificity, are provided in Table S4.

75Most compounds were predicted as inactives and the models seemed to be more in agreement in predicting inactives than active compounds. Only 757 chemicals (2.33%) are predicted as actives by more than 75% of binding models. The agreement among the binding models for the 32 k set of the prediction set is illustrated in Figure S1.

76Most categorical models (binding, agonist, and antagonist) are associated with high balanced accuracies on the ToxCast™ data (> 0.8), with no clear difference between models that used it as a training set and those that did not (see Table S3). However, for the evaluation set from the literature, the BA is clearly lower for all models (< 0.7). Nonetheless, the BA increased after removing chemicals with only one source from the literature data. This result could mean that this first filter (i.e., removing chemicals with limited information in the literature for being either positive or negative) reduced the uncertainty in the experimental data from the literature. This is in agreement with related studies showing that the results of QSAR models may change depending on the robustness of the experimental values. The second filter (i.e., removing very weak actives) also increased the BA, which suggests that the literature data may contain a number of false positives. Alternatively, the in vitro assays used by ToxCast™/Tox21 only test chemicals up to 100 μM, so very weak chemicals may not be picked up by these assays and some of the literature reports may have tested chemicals up to much higher concentrations.

77Finally, removing predictions outside the AD did not show improvement of the BA of the categorical models (see Table S3). This is in agreement with literature sources showing that predictions outside the AD are not always less accurate than those within its limits. The performance of most models showed a clear improvement of 0.05 to 0.1 on the BA after applying all the filters on the literature data to keep only the unambiguous chemicals. We believe that this effectively reduced the uncertainty of the literature sources. This step also highlighted differences between ToxCast™ and the literature data and confirmed the existence of uncertainty in the literature data. Uncertainty and data discordance was also reported in literature review of in vivo uterotrophic bioassays.

78The calculated scores for categorical models (see Table S3) take into consideration the whole prediction set (Score_1) and the accuracy of the model on its most reliable predictions (Score_2). The models that provided predictions for the whole or most of the 32 k set of chemicals, and had wide ADs, showed high Score_1 values (Umeå 0.82, OCHEM 0.83). Whereas models with predictions for smaller fractions of the prediction set and narrow AD showed better Score_2 values (UNIMIB_2 0.85, UNIBA 0.80). NIH_NCI_GUSAR (0.87 and 0.84) and FDA_NCTR_DBB (0.88 and 0.84) showed the highest values for both Score_1 and Score_2. Part of the differences among model scores could result from the uncertainty in the literature data.

79The BAs of all antagonist models was low compared with binding and agonist models (see Table S3). This may be due to the highly unbalanced training set with a low number of active antagonist chemicals. Additionally, antagonism activity (in either ToxCast™ or the literature) can be confounded with cytotoxicity because antagonist transactivation assays are loss-of-signal assays.

80The predictions of all continuous models were first converted to five classes using the list of reference chemicals as described in the evaluation set section (see “Evaluation set for continuous models” section). The predictions were then evaluated on the ToxCast™ data and the literature data to calculate the average of BA of the different evaluation steps as the score of each model (see Table S5). All models showed high BA on ToxCast™ data and relatively good BA on the evaluation set.

81The consensus predictions were first evaluated on the ToxCast™ data and then on the evaluation set from the literature. The total number of predicted active binders was 2,661 out of the 32 k set of chemicals (8.2%) based on the method described in the “Materials and Methods” section “Consensus Modeling.”

82Confusion matrices (Table 4) and prediction statistics (Table 5) revealed a clear accuracy difference between the categorical consensus for binding on the ToxCast™ data and on the evaluation set. This difference could result from the fact that the ToxCast™ data, based on a model with inputs from 18 different assays, were used by most of the models as a training set, which we presume reduces the uncertainty. This is in contrast to the literature data, where the number of sources per chemical varied from one to a few hundreds. When only the subset of the evaluation set with more than six literature sources per chemical was considered, a large increase in the sensitivity was noticed (0.23 to 0.85).

83To better understand the effect of the number of sources on the classification accuracy, ROC (receiver operating characteristic) curves were made using the fraction of the binding models in each class as a threshold for the classification predictions and increasing the number of literature sources of the evaluation set. The ROC plot shows an improvement of the classification accuracy of the consensus model as the number of sources increases (Figure 1). Note that the same level of consistency (i.e., 80%) was required to merge the sources regardless of the number of sources (see rule 1 in the “Evaluation set for categorical models” section). This could lead to the conclusion that the low classification accuracy on the full literature data is not because of a lack of accuracy of the consensus predictions, but rather to noise and experimental uncertainty in the literature data. We assume that the high number of false negatives in the confusion matrix of Table 4 is caused by false positives in the full literature data for chemicals tested only a small number of times. Thus, by considering a higher number of sources (i.e., six), the number of false positives is reduced from the evaluation set and so the number of predicted false negatives decreased. This is in agreement with what was observed in the literature.

84The first step of consensus modeling was conducted in an independent way for the categorical and continuous models on binding, agonist, and antagonist predictions. This led to a number of inconsistencies because some chemicals were predicted as active in categorical predictions but inactive in quantitative and vice versa. In addition, some chemicals were predicted as active agonists or antagonists but non-binders. To make all predictions more consistent, a number of corrections were applied on the first consensus predictions. Because the goal of this project was to help in a regulatory prioritization procedure, the modifications aimed to reduce the number of false negatives but without adding an excess of false positives. The rules that were followed to obtain the final consensus predictions are as follows:

85If a chemical i is active in the categorical consensus, then it is also considered active in the quantitative consensus.

86If a chemical i is active in the quantitative consensus and predicted as active by at least three categorical models, then it is also considered active in the categorical consensus.

87If a chemical i is predicted active by less than three categorical models, then it is considered inactive also in quantitative consensus.

88These three rules were applied on the agonist and antagonist consensus models first, then on the binding consensus. A fourth rule was added to establish consistency between agonist and antagonist consensus models and the binding consensus model.

89If a chemical i is an active agonist or active antagonist, then it is considered as active in categorical binding consensus, and its potency level in the quantitative binding consensus is made equal to its potency level as agonist/antagonist.

90An analysis of variance in concordance in each potency level of the active chemicals in the continuous models (very weak, weak, moderate, and strong) is presented as a box-plot in Figure 2. Based on this figure, we noticed a correlation between the concordance of the categorical models and the potency level of active chemicals. This implies that models are more in agreement for strong actives and that the weaker a chemical is the more difficult it is to accurately predict. Therefore, the very weak chemicals are the main source of discordance among the different in silico models and also are the most uncertain experimentally. This relationship between positive concordance (agreement between models on predictions for active chemicals) and potency level for active chemicals can be used to set a quantitative prediction to the newly reclassified active chemicals using the previously mentioned rule 1 of the corrections applied to the consensus predictions. The following thresholds were considered for each potency level:

91Strong: Concordance among models ≥ 0.9

92Moderate: 0.75 ≤ Concordance among models < 0.9

93Weak: 0.6 ≤ Concordance among models < 0.75

94Very weak: Concordance among models < 0.6

95After applying the four correction rules on consensus predictions, the total number of chemicals predicted as actives increased from 2,661 to 4,001, which corresponds to 12.3% of the total number of the prediction set (32,464). Table 6 shows the number of reclassified chemicals based on each one of the four correction rules applied to the consensus predictions. After this step, the predicted activity of several chemicals has changed. The structural information of chemicals and the predictions of the consensus model for the whole 32 k set are provided on the U.S. EPA ToxCast™ web site (see https://www3.epa.gov/research/COMPTOX/CERAPP_files.html, PredictionSet.zip).

96The confusion matrices and statistics for the binding categorical consensus model after modifications evaluated on ToxCast™ data and the literature data are presented in Table 7 and Table 8, respectively. The effect of the number of sources on the classification accuracy of the consensus model is illustrated by a bar plot in Figure S2. This figure shows an improvement of sensitivity with the increase in the number of literature sources in the evaluation set (from ~ 0.3 with at least one source to > 0.6 with six sources and more). This is translated into an increase in BA, whereas specificity is almost constant (~ 0.9) because of the high number of inactives compared to active compounds.

97The results of this project and the ToxCast™ data used as the training set are published online in the EDSP21 dashboard, together with other structural and experimental assay information (see “Consensus CERAPP QSAR ER Model Predictions” under “Chemical Summary” tab on http://actor.epa.gov/edsp21). A comparison of the single classification models to the consensus predictions for the whole 32 k set of chemicals is provided in Table S6. The calculations are done using the categorical consensus predictions as the “observed response.”

98For regulatory or prioritization purposes, one could use a looser definition of active (i.e., allow more disagreement among models) in order to further reduce the chance of false negatives. Figure 3 shows the number of chemicals that can be predicted as potential actives by the categorical consensus for binding using various positive concordance (agreement on actives between the included models) thresholds. When this threshold is set to 0.2, an additional 6,742 more chemicals can be added to the potential positives (this refers to the available binding models). This figure also shows the BA variations at different numbers of literature sources in the literature. Balanced accuracy increases as the concordance threshold increases from 0 to 0.2 because sensitivity increases (false negatives decrease) as the number of chemicals classified as active increases. For chemicals with the highest data quality (seven or more sources), the BA curve reaches a plateau at concordance thresholds of 0.4–0.5, and the number of chemicals classified as active is consistent with the number of active chemicals predicted from our consensus model (n = 4,001.) However, higher concordance thresholds result in declining BA due to increasing numbers of false positive predictions (i.e., decreasing specificity).

99The collaborative efforts of the CERAPP participants resulted in consensus predictions of the ability of chemicals to interact with ER. Up to 48 separately developed categorical and continuous models were received from 17 research groups from the United States and Europe. Separate models were built for agonist, antagonist, and binding activity. The models were applied to a large collection of 32,464 chemical structures that approximate the human exposure universe (chemicals with potential human exposure). A KNIME workflow was developed to carefully curate the large collection of chemical structures to ensure consistency in model development and evaluation. Most of the models were trained using activities derived from a data set combining 18 in vitro assays from ToxCast™ probing various points of the ER pathway. Models were then evaluated using the ToxCast™ data plus a collection of ER in vitro data from the literature. After this process, categorical predictions were combined into a consensus to classify the chemicals into actives and inactives, while continuous predictions were combined to classify the actives into 4 different potency classes: very weak, weak, moderate, and strong.

100One major observation was that most models had comparable performances, independent of the methods used, with a slight improvement for models with narrow ADs. A second and, perhaps, more important observation is that the most concordant predictions come from comparing the consensus of many models with a consensus of many literature sources. For instance, when comparing the consensus of the categorical binding models with the evaluation set from the literature for chemicals with seven or more sources, we achieve a balanced accuracy of about 90% (Table 8).

101We propose several important conclusions from our results. First, there does not appear to be an optimal modeling approach (combination of descriptor set, feature selection, or machine learning algorithm) that will solve the QSAR/docking problem and achieve perfect prediction accuracies. Second, there are inherent limitations to the accuracy of the data being used to train QSAR and docking models. Our analysis of the literature data showed a disagreement in the reported activity of many chemicals. The sources of discrepancy include limits to the concentration ranges tested, true differential activity among tissue sources [e.g., the presence of selective ER modulators, SERMs (selective estrogen receptor modulators)], and a variety of experimental artifacts and errors. Figure 2 shows that the most consistent predictions are achieved for the most potent compounds, whereas weaker compounds are called inactive by some laboratories because these compounds were not tested at a high enough concentration. So chemicals with very weak activity would be more likely to be incorrectly classified as inactive than more potent chemicals. Therefore, 100% accuracy cannot be achieved due to these limitations in the experimental data used for training and evaluation. Figures 1 and 3 help to illustrate this point by showing that higher consistency in the experimental data is associated with an increase in the concordance among model predictions. But this comes at the cost of excluding parts of the experimental data. So, just as every model has limitations, every in vitro assay also has inherent variability in its results.

102The major purpose of this study was to identify potential ER actives out of the large universe of chemicals to which humans potentially are exposed using a consensus of in silico models to overcome the limitations of single models. Most of the chemicals in this collection were predicted to be negatives, with a high agreement among the individual models. The disagreement was the highest for chemicals with weak activity (Figure 2). This disagreement is driven by the difficulties in experimentally assessing the activity of these weak chemicals. In total, the consensus predicted 4,001 chemicals as actives. The testing of these active chemicals will be prioritized from the most potent to the least according to the continuous model consensus predictions. There are 6,742 more chemicals that 20–50% of the models predicted to be positive, which could also be candidates for follow-up analyses. Although this large number of chemicals (~ 10,000 in total) appears to be a daunting set to evaluate experimentally, this is equivalent in size to the current Tox21 library already being tested for activity in ER and many other targets.

103In summary, this project demonstrates the feasibility of screening a large and toxicologically relevant library of chemical structures in an extensive battery of QSAR and docking models to meet important goals in human and environmental health. ER provides a good initial case because of the ready availability of experimental data and pre-existing models. However, through the ToxCast™ and Tox21 programs, and through other large scale data-integration projects, equivalently large data sets will become available for other multiple targets of environmental importance.