Kamel Mansouri

Ahmed Abdelaziz

Aleksandra Rybacka

Alessandra Roncaglioni

Alexander Tropsha

Alexandre Varnek

Alexey Zakharov

Andrew Worth

Ann M Richard

Christopher M Grulke

Richard S Judson
email: judson.richard@epa.gov

U S Epa

Daniela Trisciuzzi

Denis Fourches

Dragos Horvath

Emilio Benfenati

Eugene Muratov

Eva Bay Wedebye

Francesca Grisoni

Giuseppe F Mangiatordi

Giuseppina M Incisivo

Huixiao Hong

Hui W Ng

Igor V Tetko

Ilya Balabin

Jayaram Kancherla

Jie Shen

Julien Burton

Marc Nicklaus

Matteo Cassotti

Nikolai G Nikolov

Orazio Nicolotti

Patrik L Andersson

Qingda Zang

Regina Politi

Richard D Beger

Roberto Todeschini

Ruili Huang

Sherif Farag

Sine A Rosenberg

Svetoslav Slavov

Xin Hu

CERAPP: Collaborative Estrogen Receptor Activity Prediction Project

Keywords:

à 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.

Introduction

There are tens of thousands of natural and synthetic chemical substances to which humans and wildlife are exposed [START_REF] Dionisio | Exploring consumer exposure pathways and patterns of use for chemicals in the environment[END_REF][START_REF] Egeghy | The exposure data landscape for manufactured chemicals[END_REF][START_REF] Judson | The toxicity data landscape for environmental chemicals[END_REF]. A subset of these compounds may disrupt normal functioning of the endocrine system and cause health hazards to both humans and ecological species [START_REF] Birnbaum | Cancer and developmental exposure to endocrine disruptors[END_REF][START_REF] Diamanti-Kandarakis | Endocrinedisrupting chemicals: an Endocrine Society scientific statement[END_REF][START_REF] Mahoney | Developmental programming: impact of fetal exposure to endocrine disrupting chemicals on gonadotropinreleasing hormone and estrogen receptor mRNA in sheep hypothalamus[END_REF]UNEP and WHO 2013). Endocrinedisrupting 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 [START_REF] Diamanti-Kandarakis | Endocrinedisrupting chemicals: an Endocrine Society scientific statement[END_REF]). Exposure to EDCs can lead to adverse health effects involving developmental, neurological, reproductive, metabolic, cardiovascular, and immune systems in humans and wildlife [START_REF] Colborn | Developmental effects of endocrine-disrupting chemicals in wildlife and humans[END_REF][START_REF] Davis | Medical hypothesis: xenoestrogens as preventable causes of breast cancer[END_REF][START_REF] Diamanti-Kandarakis | Endocrinedisrupting chemicals: an Endocrine Society scientific statement[END_REF].

The estrogen receptor (ER) is one of the most extensively studied targets related to the effects of EDCs [START_REF] Mueller | Estrogen receptors and endocrine diseases: lessons from estrogen receptor knockout mice[END_REF][START_REF] Shanle | Endocrine disrupting chemicals targeting estrogen receptor signaling: identification and mechanisms of action[END_REF]. This concern about estrogen-like activity of man-made chemicals is because of their potential for negatively affecting reproductive function [START_REF] Hileman | Environmental estrogens linked to reproductive abnormalities, cancer[END_REF][START_REF] Kavlock | Research needs for the risk assessment of health and environmental effects of endocrine disrupters: a report of the U.S. EPA-sponsored workshop[END_REF]. The emergence of concerns about EDCs has resulted in regulations requiring assessment of chemicals for estrogenic activity [START_REF] Adler | Alternative (nonanimal) methods for cosmetics testing: current status and future prospects-2010[END_REF]; U.S. Environmental Protection Agency (EPA) 1996; U.S. Food and Drug Administration (FDA) 1996]. There are numerous in vitro and in vivo protocols to identify potential endocrine pathway-mediated effects of chemicals, including interactions with hormone receptors [START_REF] Jacobs | The use of metabolising systems for in vitro testing of endocrine disruptors[END_REF][START_REF] Rotroff | Using in vitro high throughput screening assays to identify potential endocrine-disrupting chemicals[END_REF] volume 124 | number 7 | July 2016 • Environmental Health Perspectives 2013; [START_REF] Shanle | Endocrine disrupting chemicals targeting estrogen receptor signaling: identification and mechanisms of action[END_REF][START_REF] Sung | Detection of endocrine disruptors-from simple assays to whole genome scanning[END_REF]). 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 [START_REF] Hubal | Advancing exposure characterization for chemical evaluation and risk assessment[END_REF][START_REF] Egeghy | The exposure data landscape for manufactured chemicals[END_REF][START_REF] Judson | The toxicity data landscape for environmental chemicals[END_REF].

The 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 [START_REF] Muster | Computational toxicology in drug development[END_REF][START_REF] Vedani | In silico toxicology in drug discovery-concepts based on three-dimensional models[END_REF]. QSARs have been promoted and their use recognized since the pioneering work of Hansch in the 1960s [START_REF] Fujita | A new substituent constant, π, derived from partition coefficients[END_REF][START_REF] Hansch | Correlation of biological activity of phenoxyacetic acids with Hammett substituent constants and partition coefficients[END_REF][START_REF] Hansch | The structure-activity relationship in amides inhibiting photosynthesis[END_REF]. 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 [START_REF] Dearden | How not to develop a quantitative structure-activity or structure-property relationship (QSAR/QSPR)[END_REF][START_REF] Worth | The Characterisation of (Quantitative) Structure-Activity Relationships: Preliminary Guidance[END_REF].

Recently, in vitro high-throughput screening (HTS) assays have emerged and become a viable tool for large-scale chemical testing [START_REF] Judson | Estimating toxicity-related biological pathway altering doses for high-throughput chemical risk assessment[END_REF][START_REF] Kavlock | Computational toxicology as implemented by the U.S. EPA: providing high throughput decision support tools for screening and assessing chemical exposure, hazard and risk[END_REF][START_REF] Wetmore | Integration of dosimetry, exposure, and highthroughput screening data in chemical toxicity assessment[END_REF]. 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.

With 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 (U.S. EPA 2015). 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 [START_REF] Dix | The ToxCast program for prioritizing toxicity testing of environmental chemicals[END_REF][START_REF] Judson | In vitro screening of environmental chemicals for targeted testing prioritization: the ToxCast project[END_REF][START_REF] Martin | Impact of environmental chemicals on key transcription regulators and correlation to toxicity end points within EPA's ToxCast program[END_REF]). Relevant to this paper, these two projects have in common ~ 1,800 chemicals tested in a battery of 18 ER-related assays [START_REF] Huang | Profiling of the Tox21 10K compound library for agonists and antagonists of the estrogen receptor alpha signaling pathway[END_REF][START_REF] Judson | Integrated model of chemical perturbations of a biological pathway using 18 in vitro high throughput screening assays for the estrogen receptor[END_REF].

This 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 [START_REF] Judson | Integrated model of chemical perturbations of a biological pathway using 18 in vitro high throughput screening assays for the estrogen receptor[END_REF]. 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.

Materials and Methods

Participants and Project Planning

The 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.

Data Sets

Provided 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 [START_REF] Dix | The ToxCast program for prioritizing toxicity testing of environmental chemicals[END_REF][START_REF] Huang | Profiling of the Tox21 10K compound library for agonists and antagonists of the estrogen receptor alpha signaling pathway[END_REF][START_REF] Judson | In vitro screening of environmental chemicals for targeted testing prioritization: the ToxCast project[END_REF]). 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 [START_REF] Judson | Integrated model of chemical perturbations of a biological pathway using 18 in vitro high throughput screening assays for the estrogen receptor[END_REF]U.S. EPA 2014c). 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 [START_REF] Judson | Integrated model of chemical perturbations of a biological pathway using 18 in vitro high throughput screening assays for the estrogen receptor[END_REF]. 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 [START_REF] Judson | Integrated model of chemical perturbations of a biological pathway using 18 in vitro high throughput screening assays for the estrogen receptor[END_REF]. 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.

Prediction 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: • Chemicals 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 [START_REF] Dionisio | Exploring consumer exposure pathways and patterns of use for chemicals in the environment[END_REF][START_REF] Judson | ACToR-Aggregated Computational Toxicology Resource[END_REF][START_REF] Judson | Aggregating data for computational toxicology applications: the U.S. Environmental Protection Agency (EPA) Aggregated Computational Toxicology Resource (ACToR) system[END_REF]; U.S. EPA 2014a). [START_REF] Gaulton | ChEMBL: a large-scale bioactivity database for drug discovery[END_REF]).

The 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) (U.S. EPA 2016). 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 AC 50 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 AC 50 value if the compound was reported as active. All inactive compounds were arbitrarily assigned an AC 50 value of 1 M.

Chemical Structure Curation

Chemical 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 QSARready 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 (Fourches et al. 2010;[START_REF] Wedebye | Use of QSAR to Identify Potential CMR Substances of Relevance under the REACH Regulation[END_REF]. The free and opensource 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 [START_REF] Berthold | KNIME: the Konstanz Information Miner[END_REF]. The workflow performed a series of curation steps: 1) The 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 (NIH 2015), ChemSpider (Royal Society of Chemistry 2015)] or removed if ambiguous.

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

3) The structures were desalted, and inorganic counterions were removed. 4) The 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. 5) Valid 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 (ChemAxon 2014;[START_REF] Reusch | Reaction Examples. Examples of Organic Reactions[END_REF][START_REF] Sitzmann | Tautomerism in large databases[END_REF]). 6) These transformations were followed by neutralizing the charged structures, when possible, and removing the stereochemistry information. 7) Explicit hydrogen atoms were added, and structures were aromatized according to

Hückel's rules implemented in KNIME [START_REF] Berthold | KNIME: the Konstanz Information Miner[END_REF]). 8) The duplicates were removed using the IUPAC (International Union of Pure and Applied Chemistry) InChI (International Chemical Identifier) codes because these are unequivocal identifiers. 9) The final filter was applied to remove chemicals containing metals that often cause problems in molecular descriptor calculations.

Both 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 setremained 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) (U.S. EPA 2016). The identity of these chemicals (name, CASRN) was not provided to the participating modeling groups during the modeling process.

Modeling Approaches

The 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 wellknown and innovative methods including partial least-squares (PLS) [START_REF] Ståhle | Partial least squares analysis with cross-validation for the two-class problem: a Monte Carlo study[END_REF][START_REF] Wold | PLS-regression: a basic tool of chemometrics[END_REF], partial least-squares discriminant analysis (PLS-DA) [START_REF] Frank | A statistical view of some chemometrics regression tools[END_REF][START_REF] Nouwen | Classification of environmentally occurring chemicals using structural fragments and PLS discriminant analysis[END_REF], decision forest (DF) [START_REF] Hong | An in silico ensemble method for lead discovery: decision forest[END_REF][START_REF] Hong | Multiclass Decision Forest-a novel pattern recognition method for multiclass classification in microarray data analysis[END_REF][START_REF] Tong | Decision forest: combining the predictions of multiple independent decision tree models[END_REF][START_REF] Xie | Decision forest analysis of 61 single nucleotide polymorphisms in a casecontrol study of esophageal cancer: a novel method[END_REF], three-dimensional (3D) quantitative spectral data-activity relationship (QSDAR) (Beger et al. 2001;Beger and Wilkes 2001;[START_REF] Slavov | Complementary PLS and KNN algorithms for improved 3D-QSDAR consensus modeling of AhR binding[END_REF], support vector machines (SVM) [START_REF] Cristianini | An Introduction to Support Vector Machines and Other Kernelbased Learning Methods[END_REF], k nearest neighbors (kNN) [START_REF] Cover | Nearest neighbor pattern classification[END_REF][START_REF] Kowalski | The K-Nearest Neighbor Classification Rule (pattern recognition) applied to nuclear magnetic resonance spectral interpretation[END_REF], associative artificial neural networks (ASNN) (Tetko 2002a(Tetko , 2002b)), PASS algorithm derived from Naïve Bayes classifier [START_REF] Poroikov | Robustness of biological activity spectra predicting by computer program PASS for noncongeneric sets of chemical compounds[END_REF], self-consistent regression with radial basis function interpolation (RBF-SCR) [START_REF] Zakharov | A new approach to radial basis function approximation and its application to QSAR[END_REF], OCHEM machine learning methods [START_REF] Tetko | How accurately can we predict the melting points of druglike compounds?[END_REF], docking and consensus of different approaches [START_REF] Horvath | An evolutionary optimizer of libsvm models[END_REF][START_REF] Ng | Competitive molecular docking approach for predicting estrogen receptor subtype α agonists and antagonists[END_REF][START_REF] Sushko | Online Chemical Modeling Environment (OCHEM): web platform for data storage, model development and publishing of chemical information[END_REF]. 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 (Judson et al. 2015), but some pre-existing models were also used that had been trained using other data sets from the literature such as METI (2002). 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 [START_REF] Roberts | LeadScope: software for exploring large sets of screening data[END_REF], PaDEL-Descriptor (Yap 2011), QikProp (version 3.4, http:// www.schrodinger.com/QikProp/), multilevel and quantitative neighborhoods of atoms (MNA, QNA) used by GUSAR and PASS [START_REF] Filimonov | QNA-based 'Star Track' QSAR approach[END_REF][START_REF] Poroikov | Robustness of biological activity spectra predicting by computer program PASS for noncongeneric sets of chemical compounds[END_REF], DRAGON (Talete srl 2012), Mold2 [START_REF] Hong | Mold 2 , molecular descriptors from 2D structures for chemoinformatics and toxicoinformatics[END_REF][START_REF] Hong | Mold 2 molecular descriptors for QSAR[END_REF], GLIDE (version 6.5, http:// www.schrodinger.com/Glide), AutoDock [START_REF] Goodsell | Automated docking of flexible ligands: applications of AutoDock[END_REF], ISIDA [START_REF] Varnek | ISIDA -Platform for virtual screening based on fragment and pharmaco phoric descriptors[END_REF], and other fingerprint generators. Some of the participants applied feature selection techniques, such as genetic algorithms (GAs) [START_REF] Davi | Handbook of Genetic Algorithm[END_REF]) and random forest (RF) [START_REF] Breiman | Random forests[END_REF]). These techniques were applied after calculating descriptors to reduce collinearity and variable dimensionality to keep only the most informative descriptors in the models.

volume 124 | number 7 | July 2016 • Environmental Health Perspectives

Evaluation Procedure for the Categorical and Continuous Models

All 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) (U.S. EPA 2016).

Standard 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).

Evaluation 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: 1) If, 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. 2) If a chemical was an active agonist or antagonist, it also was considered as an active binder if the information was not available. 3) If a chemical was an inactive agonist and inactive antagonist, it was considered also as nonbinder if the information was not available. This 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) (U.S. EPA 2016).

Evaluation set for continuous models. For active chemicals with available quantitative information from concentration-response assays, the log 10 -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 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 [START_REF] Judson | Integrated model of chemical perturbations of a biological pathway using 18 in vitro high throughput screening assays for the estrogen receptor[END_REF]. 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: • Strong: Activity concentration below 0.09 μM. • Moderate: Activity concentration between 0.09 and 0.18 μM. • Weak: Activity concentration between 0.18 and 20 μM. • Very Weak: Activity concentration between 20 and 800 μM. • Inactive: Activity concentration higher than 800 μM. The 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. The number of entries in each class for binding, agonist, and antagonist are summarized in Table 3.

Evaluation 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).

• Chemicals 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). • All chemicals in the full literature data (all literature sources combined). • All chemicals with at least two literature sources. • All chemicals from the literature data excluding the very weak actives. • Chemicals within the applicability domain (AD) of each model (if provided). • Chemicals 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).

To 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.

Score_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.

N BA N pred N total N pred N total BA N pred _ 1 3 1 ToxCast ToxCast _ ToxCast _ _ N _ i i _ i i 1 N filters filters # # + + = 1 score = c p | [1]
where 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, N filters represents the number of five filters applied to the evaluation set chemicals and i the steps 2, 3, 4, 5, and 6.

Score_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).

score_2 = 1 ⁄ 2 (BA ToxCast + BA all filters ) [2]
The 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.

Consensus Modeling

The 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.

Considering 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.

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

w score score 1 i i j j n = = / | [3]
where n is the total number of models that provided predictions for the chemical i, and score j is the score of the jth model predicting chemical i.

Next, the consensus potency level C i of each chemical was determined using the predicted potency classes P j of the n available models and their corresponding weights w as follows (Equation 4): The classification into actives and inactives is based on a consensus between the literature data sources that were in agreement.

C w P 1 i j j n j # = = | [4]

Results and Discussion

Models and Evaluation

A 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 andS5) 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.

The 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) (U.S. EPA 2016).

The 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.

Most 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.

Most 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 [START_REF] Steinmetz | Methods for assigning confidence to toxicity data with multiple values-identifying experimental outliers[END_REF]). 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.

Finally, 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 [START_REF] Sahigara | Comparison of different approaches to define the applicability domain of QSAR models[END_REF]. 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 [START_REF] Kleinstreuer | A curated database of rodent uterotrophic bioactivity[END_REF].

The 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.

The 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.

The 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.

Consensus Model

The 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."

Confusion 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 The literature data with more than six sources represents the most consistent part of the evaluation set.

1029 per chemical was considered, a large increase in the sensitivity was noticed (0.23 to 0.85).

To 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 [START_REF] Steinmetz | Methods for assigning confidence to toxicity data with multiple values-identifying experimental outliers[END_REF].

Corrections to the Consensus Model

The 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: 1) If a chemical i is active in the categorical consensus, then it is also considered active in the quantitative consensus. 2) If 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. 3) If a chemical i is predicted active by less than three categorical models, then it is considered inactive also in quantitative consensus.

These 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. 4) If 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. An 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:

Concordance between models

Very weak Moderate

Potency of active chemicals

Strong Weak volume 124 | number 7 | July 2016 • Environmental Health Perspectives After 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) (U.S. EPA 2016).

The 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 andTable 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. The 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) (U.S. EPA 2014c). 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."

For 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).

Conclusion

The 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. One 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).

We 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 and3 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.

The 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.

In 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.

  • Strong: Average score > 0.75 • Moderate: 0.5 < Average score between ≤ 0.75 • Weak: 0.25 < Average score ≤ 0.5 • Very weak: 0 < Average score ≤ 0.25 • Inactive: Average score = 0

Figure 1 .

 1 Figure 1. ROC curves of the categorical corrected consensus predictions for binding evaluated against different sets of the evaluation set with variable numbers of literature sources. The number of available chemicals in the evaluation set (between brackets) decreased with higher numbers of literature sources. The true and false positive rates are determined based on the number of actives in the different sets of the evaluation set.

Figure 2 .

 2 Figure 2. Box-plot of the positive class potency levels in the corrected quantitative consensus predictions for binding. The concordance between models is the fraction of the number of models that agrees on the prediction of a certain chemical. Boxes extend from the 25th to the 75th percentile, horizontal bars represent the median, whiskers indicate the 10th and 90th percentiles, and outliers are represented as points.

Figure 3 .

 3 Figure 3. Variation of the balanced accuracy of the corrected categorical consensus predictions for binding with positive concordance (agreement between models on predictions for active chemicals) threshold at different numbers of literature sources.

  and https://www3.epa. gov/research/COMPTOX/CERAPP_files. html, EvaluationSet.zip) (U.S. EPA 2016).

Table 1 .

 1 Methods adopted by the participant groups (alphabetic order) in the modeling procedure.

	Model name	Calibration method	Descriptors software/type	Training set (No. of chemicals)	Predictions type
	DTU	PLS/fragments	Leadscope	METI (595,481)/ToxCast ™ (1,422)	Categorical
	EPA_NCCT	GA + PLSDA	PADEL	ToxCast ™ (1,529)	Categorical
	FDA_NCTR_DBB (Ng et al. 2014)	DF	Mold2	ToxCast ™ (1,677)	Categorical
	FDA_NCTR_DSB	PLS	3D-SDAR	ToxCast ™ (1019)	Categorical
	ILS_EPA (Zang et al. 2013)	SVM + RF	Qikprop	ToxCast ™ (1,677)	Categorical
	IRCCS_CART (Roncaglioni et al. 2008)	CART-VEGA	2D descriptors	METI (806)	Categorical
	IRCCS_Ruleset	Ruleset	SMARTS	ToxCast ™ (1,529)	Categorical
	JRC_Ispra (Poroikov et al. 2000)	PASS	MNA	-	Categorical
	Lockheed Martin	kNN	Fingerprints	ToxCast ™ (1,677)	Categorical + continuous
	NIH_NCATS	Docking	AutoDock score	-	Categorical
	NIH_NCI_GUSAR (Filimonov et al. 2009)	RBF-SCR	MNA, QNA	ToxCast ™ (1,677)	Categorical
	NIH_NCI_PASS (Poroikov et al. 2000)	PASS	MNA	ToxCast ™ (1,677)	Categorical
	OCHEM (2015)	Consensus	11 Descriptor types	ToxCast ™ (1,660)	Categorical + continuous
	RIFM	SVM	Fingerprints	ToxCast ™ (1,677)	Categorical
	Umeå (Rybacka et al. 2015)	ASNN	DRAGON	METI + (Kuiper et al. 1997; Taha et al. 2010)	Categorical
	UNC_MML	SVM+RF	DRAGON	ToxCast ™ (120)	Categorical
	UNIBA (Trisciuzzi et al. 2015)	Docking	GLIDE score	ToxCast ™ (1,677)	Categorical
	UNIMIB	kNN	DRAGON + fingerprints	ToxCast ™ (1,677)	Categorical
	UNISTRA (Horvath et al. 2014)	SVM	ISIDA	ToxCast ™ (1,529)	Categorical + continuous

Predictions type: 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. Calibration methods: PLS (partial least-squares), PLS-DA (partial least-squares discriminant analysis), SVM (support vector machines), RF (random forest), DF (Decision forest), kNN (k nearest neighbors), ASNN (associative artificial neural networks), PASS (algorithm derived from Naïve Bayes classifier), RBF-SCR (self-consistent regression with radial basis function interpolation).

Table 2 .

 2 Evaluation set for binary categorical models. Distribution of the number of active and inactive chemicals within the three different classes: binding, agonists and antagonists.

	Class/activity	Active	Inactive	Total
	Binding	1,982	5,301	7,283
	Agonist	350	5,969	6,319
	Antagonist	284	6,255	6,539
	Total	2,017	7,024	7,522

Table 3 .

 3 Evaluation set for quantitative models. Distribution of the number of chemicals in the five potency levels within the three different classes (binding, agonists, and antagonists), classifications based on average scores.

	Class/activity	Inactive	Very weak	Weak	Moderate	Strong	Total
	Binding	5,042	685	894	72	77	6,770
	Agonist	5,892	19	179	31	42	6,163
	Antagonist	6,221	76	188	10	10	6,505
	Total	6,892	702	916	81	93	7,253

The classification of the chemicals in the five potency levels is based on the concentration responses from the literature sources that were in agreement. volume 124 | number 7 | July 2016 • Environmental Health Perspectives

Table 4 .

 4 Confusion matrices of categorical consensus predictions for binding.

				Literature evaluation set	Literature evaluation set
		ToxCast ™ data	ToxCast ™ data	(all: 7,283)	(all: 7,283)
	Observed/predicted	predicted actives	predicted inactives	predicted actives	predicted inactives
	Observed actives	76	13	467	1,515
	Observed inactives	25	1,415	268	5,033

Table 5 .

 5 Statistics of categorical consensus predictions for binding on ToxCast™ and literature data.

			Literature evaluation set	Literature evaluation set
	Statistics/used data	ToxCast ™ data	(all: 7,283)	(> 6 sources: 1,257)
	Sensitivity	0.85	0.23	0.85
	Specificity	0.98	0.95	0.97
	Balanced accuracy	0.92	0.59	0.91

Table 6 .

 6 Number of chemicals reclassified after applying each one of the four prediction correction rules.Rule 1: Chemicals that changed from inactive to active in the quantitative consensus based on the categorical consensus. Rule 2: Chemicals that changed from inactive to active in the categorical consensus based on the quantitative consensus. Rule 3: Chemicals that changed from active to inactive in the quantitative consensus based on the predictions of the categorical consensus. Rule 4: Chemicals that changed from inactive to active in the categorical binding consensus based on their agonist and antagonist activity in the categorical consensus.

	Rule used for		Rule 1			Rule 2			Rule 3		Rule 4
	each class	Agonist Antagonist Binding Agonist Antagonist Binding Agonist Antagonist Binding Binding
	Number of	1,288	2,760	1,587	217	14	344	145	161	38	966
	chemicals

Table 7 .

 7 Confusion matrices of the modified categorical consensus predictions for binding.

				Literature evaluation set	Literature evaluation set
		ToxCast ™ data	ToxCast ™ data	(All: 7,283)	(All: 7,283)
	Observed/predicted	predicted actives	predicted inactives	predicted actives	predicted inactives
	Observed actives	83	6	597	1,385
	Observed inactives	40	1,400	463	4,838

Table 8 .

 8 Statistics of the modified categorical consensus for binding predictions on ToxCast™ and literature data.

		ToxCast ™	Literature evaluation set	Literature evaluation set
	Statistics/used data	data	(All: 7,283)	(> 6 Sources: 1,275)
	Sensitivity	0.93	0.30	0.87
	Specificity	0.97	0.91	0.94
	Balanced accuracy	0.95	0.61	0.91


J R Khuntia

Sébastien Proust

K K Khatua

Unsteady open-channel flows over rough bed with and without emergent rigid vegetation: A laboratory experiment

published or not. The documents may come   L'archive ouverte pluridisciplinaire

INTRODUCTION

The study of unsteady open-channel flows is very important for assessing and predicting the flood hazard (flow depth, velocity, time of submersion) and to demarcate the area of inundation in the case of flash floods or floods in rivers of small and moderate size. What is specific to unsteady river flows is a hysteresis in the stage discharge curve, with a river flow conveyance higher during the rising water level than during the falling water level. The estimate of river flow conveyance under unsteady state is expected to be more complex in the presence of emergent vegetation (along floodplains, in the presence of riparian vegetation).

The importance of the study is to predict the complexity of the flow behavior due to unsteadiness under extreme flood events or flash floods. The combined effect of bed roughness and emergent rigid vegetation makes the flow process complicated in real field applications. So, the present laboratory study simplifies the real field case.

Unsteady open-channel flows over rough bed with and without emergent rigid vegetation: A laboratory experiment

S. Proust

Researcher, Riverly, INRAE, Centre de Lyon-Villeurbanne, 5 rue de la Doua, Villeurbanne, France

K. K. Khatua

Professor, Department of Civil Engineering, National Institute of Technology Rourkela, India ABSTRACT: We report here a laboratory study on the effects of bed roughness and emergent rigid vegetation on the unsteady flows associated with flash floods and extreme flood events. The flow structure of unsteady open-channel flows over rough bed with and without emergent rigid vegetation has been investigated in an 18m long and 3m wide laboratory flume. Steady uniform flows were also studied and served as reference flows. For both steady and unsteady flows, two geometries are tested: (1) uniform bed roughness (uniform dense synthetic grass modelling meadow) along the flume and (2) a uniform staggered distribution of emergent wooden circular cylinders (model of rigid vegetation) set on the bed roughness. Transient flow depths are simultaneously measured at six longitudinal positions using ultra-sonic sensors. Transient velocities are measured at one longitudinal position over the water column using a side looking ADV probe to estimate depth-averaged velocity. In order to compute ensemble averages of the flow parameters, 109 runs of the same hydrograph are carried out repeatedly at the flume entrance. Two consecutive runs are separated by a base flow. The ensemble averages of the measured discharge, flow depths and velocity are found to be converged after 72 runs. For the two geometries: at the peak flow of the hydrograph, the vertical profiles of mean streamwise velocity and horizontal Reynolds shear stress are comparable to the profiles obtained for steady uniform cases of same flow depth; and accelerated and decelerated velocity profiles are obtained during the rising limb and falling limb of the hydrograph, respectively. Reynolds shear stresses are also found to be higher during the rising limb than during the falling limb for a fixed flow depth. The hysteresis in the depth-averaged velocity / flow depth relationships is comparable in shape and size for the two geometries highlighting the weak effect of the rigid vegetation compared to the effects of un-stationarity.

Over a rough bed, the vertical distribution of mean streamwise velocity is often approximated by a logarithmic law profile along a big part of boundary layer flows [START_REF] Monin | Statistical Fluid Mechanics: Mechanics of turbulence[END_REF]. In the lower layer (z/H < 0.3) of shallow flows over rough surfaces found in rivers, bottom roughness may affect the development of the log-shaped velocity profile. [START_REF] Nikora | Turbulence characteristics of New Zealand gravel-bed rivers[END_REF] found deviations from the log-law flow in the roughness layer (z/H = 0.1 to 0.3) which is a 3D inner layer of the flow in gravel-bed rivers. In such flows, a velocity inflection has been observed, resulting in so-called s-shaped profiles [START_REF] Marchand | Velocity profile, water-surface slope, and bedmaterial size for selected streams in Colorado[END_REF][START_REF] Bathurst | Velocity profile in high-gradient, boulder-bed channels[END_REF][START_REF] Ferro | Flow velocity profiles in gravel-bed rivers[END_REF]. [START_REF] Franca | The simultaneous occurrence of logarithmic and S-shaped velocity profiles in gravel-bed river flows[END_REF] observed velocity profiles in two shallow gravel-bed rivers with low relative submergence. The log-shaped profiles and S-shaped profiles were simultaneously observed in these flows indicating the possibility of 3D flow inside the roughness layer (k/H >≈ 0.6, where k is the equivalent sand grain size of the bed roughness and H is the flow depth). [START_REF] Bombar | Calculation of the time-varying mean velocity by different methods and determination of the turbulence intensities[END_REF] computed the mean velocities and turbulence intensities with 15 runs of the same triangular shape hydrograph, based on laboratory experiments. The experiments were conducted in a rectangular flume (18m × 0.7m × 05m) with a bed slope of 0.004. The bed of the flume was composed of nonuniform sediment mixture with d 50 = 0.43 and thickness was 20mm. The mean velocities and turbulence intensities were found larger in the rising limb of the hydrograph as compared to the falling limb for the same flow depth. [START_REF] Alfadhli | Velocity distribution in non-uniform/unsteady flows and the validity of log law[END_REF] investigated the vertical profiles of longitudinal mean velocity in steady and unsteady non-uniform open channel flows. It was found that for steady/unsteady flow the Log law is applicable only in the inner region (z/H < 0.2) when compared with the measured longitudinal velocity. [START_REF] Dupuis | Combined effects of bed friction and emergent cylinder drag in open channel flow[END_REF] experimentally investigated steady open-channel flows subjected to a longitudinal transition in hydraulic roughness, from bed friction to emergent cylinder drag and vice versa. The vertical profiles of longitudinal mean velocity in the cylinder alignment, the velocity bulge and constant velocity region were clearly observed in the near-bed region at z/D = 1.9 and at z/D = 4.1 respectively (where D is the diameter of cylinder). [START_REF] Tognin | Experimental Setup and Measuring System to Study Solitary Wave Interaction with Rigid Emergent Vegetation[END_REF] studied the propagation of a solitary wave in the presence of emerging rigid cylinders mimicking the propagation of tsunamis through mangrove forests. Specifically, it was observed that the vegetation reduces the wave height proportionally to the vegetation density. During base flow conditions, the vertical profile of longitudinal velocity clearly showed a logarithmic profile in bare soil condition, whereas it turned into a rather uniform velocity distribution in the presence of vegetation.

In the literature, the separated and combined effects of a rough bed and emergent rigid vegetation on unsteady flow parameters (mean velocity and turbulence statistics) has not been studied so far, to the authors' knowledge. The present study therefore focuses on these issues. A particular attention was paid to the convergence of flow parameters based on ensemble averages:

(1) to confirm the mean velocity profiles and Reynolds shear stress for unsteady accelerated, and decelerated flows over a rough bed; and (2) to study the vertical distribution of mean flow and Reynolds shear stress for unsteady accelerated and decelerated over a rough bed with emergent rigid vegetation.

EXPERIMENTAL SETUP AND PROCEDURES

Flume and measuring devices

The experiments have been conducted in a rectangular open-channel flume (18 m long and 3 m wide) in the Hydraulics and Hydro-morphology laboratory of the RiverLy research Unit of INRAE Lyon-Villeurbanne, France (Figs. 1 to 3). The longitudinal bed slope of the flume is S 0 = 1.05×10 -3 m/m. Because of some laboratory constraints, the experiments are performed in the left-hand part of the flume, using 1/3 of the total flume width (working width = 1 m, working length = 18 m). A vertical Plexiglas sidewall is used to isolate a 1 m-wide open-channel on the left-hand side of the flume. Figure 2a shows the side views of the channel with a rough bed only (M = Meadow) and with emergent circular cylinders set on the rough bed (W = Woodland), and Figure 2b shows top views of the array of cylinders and elementary pattern of cylinders with the ADV probe position for velocity measurements. Figure 3 shows the channel view of the two configurations studied (W and M). The 1m wide open-channel is supplied by a 1 m wide inlet tank, the inflow being monitored using a flow meter (Water flux) at an acquisition rate of 50 Hz. The glass bed has been covered by: a) artificial dense grass (M); b) a staggered distribution of wooden circular cylinders set on the artificial grass (W). The Manning's roughness coefficient n of the rough bed surface is 0.0157 s/m 1/3 . A Cartesian right-handed coordinate system is used in which x-, y-, and z-axes are aligned with the longitudinal (parallel to the flume bottom), transverse, and vertical (normal to the flume bottom). The water surface elevations are measured with six ultrasonic sensors (BAUMER UNDK 20I6912/S35A) located at longitudinal positions x = 3 m, 6 m, 10 m, 13 m, 15 m and 17 m. At each measuring point, elevation is recorded at a rate of 50 Hz with an accuracy better than0.3 mm. Velocities are measured over the water column at x = 13m, using a side-looking ADV probe (Vectrino plus, Nortek) with a sampling rate of 100 Hz and an accuracy of 0.5% of the measured velocity. The measurement is taken at one fixed position (x = 13m) for the two configurations M and W. water depth h p and a falling limb until the next base flow (Fig. 4). The associated unsteadiness flow parameter  (Eq. 1, where T r is the duration of the rising limb and S o is the longitudinal slope of the bed) equals to 0.65 (same order of magnitude of hydrographs as those carried out by [START_REF] Lai | Experiments on flood-wave propagation in compound channel[END_REF]. The unsteadiness parameter  represents the ratio of the rising speed of water surface to the wave speed at the peak flow.

 ℎ 𝑝 -ℎ 𝑏 𝑆 𝑜 𝑇 𝑟√ 𝑔ℎ 𝑝                   
Due to software constraints, four sequences of twenty four runs are launched (Fig. 4a), and eventually a sequence of 13 runs. At each measuring point, the 109 runs are carried out without switching off the pump to ensure the repeatability of the process.

Unsteady data treatment

(a) (c)

Second, to perform ensemble averaging of water depth and velocity data, a windowing of instantaneous data is carried out by (1) detecting the peak discharge value of instantaneous flow rate Q (t) and ( 2) considering a window centered on the detected peak value. The window width is chosen to include the base flow (i.e., 7 l/s) before and after the hydrograph, and the rising and falling limbs (Fig. 4b).

Third, an ensemble averaging of the 109 runs is computed (Fig. 4c). If P(t) is the instantaneous value of a flow parameter P, the ensemble average < P(t) > and the associated standard deviation are defined by < 𝑃(𝑡) >= (2)

σ(t) = √ 1 109 ∑ (P i (t)-< 𝑃(𝑡) >) 2 109 i=1 (3)
Note that, a 5 second moving time average is used to smoothen the flow parameters.

3 EXPERIMENTAL RESULTS

Convergence of the flow parameters

We denote Q final , H final and U xfinal as the ensemble averages of discharge Q, flow depth H, and local longitudinal velocity U x , respectively, using the 109 runs. Figure 5 shows the convergence of Q, H and U x as a function of the number of runs. The convergence of the ensemble averages of discharge Q, flow depth H, and local streamwise velocity U x was obtained after 45 runs, 50 runs, and 72 runs, respectively.

Mean streamwise velocity profiles

The vertical profiles of streamwise velocity U x are found to be superimposed when comparing with the steady uniform flow case and the base flows of same discharge at 7 l/s. This is happening for both the geometries studied (Fig. 6a,b) (i.e., uniform meadow M, uniform woodland W). It is seen in section 3.1 that after 72 repetitions of the same hydrograph, the ensemble average of the velocity become converged. For both the two configurations of the fixed flow (i.e., 15 l/s), and for both the accelerated and decelerated longitudinal velocity profiles are deviating from the uniform flow cases (Fig. 6c,d). A constant velocity region was observed in case of emergent rigid vegetation (Fig. 6d). Higher values are found during the rising cases as compared to that during the falling limb cases. The longitudinal velocity profiles of steady uniform flow case are found to be very similar to the peak flow (i.e., 22 l/s) of the unsteady case for the both configurations studies (Fig. 6e,f). The vertical profiles of Reynolds shear stress are found to be superimposed when comparing with the steady uniform flow case and the base flows of the same discharge i.e. at 7 l/s. This is valid for the both the geometries studied (Fig. 7a,b) (i.e., uniform meadow M, uniform woodland W). For these two configurations, at a fix flow of 15 l/s, both the accelerated and decelerated Reynolds shear stresses are found to be deviating from the results of uniform flow cases. Higher values are found to be occurring during the rising limb cases as compared to those values during the falling limb cases (Fig. 7c,d). The Reynolds shear stresses of the steady uniform flow cases are also found to be very similar to the peak flow of unsteady cases for the both configurations studies (Fig. 7e,f).

Hysteresis loops

Figure 8 shows the clockwise hysteresis loops in the depth-averaged longitudinal velocity / flow depth for both uniform meadow M and uniform woodland W configurations. For both geometries, given a flow depth value, the flow conveyance is higher during the rising limb than during the falling limb. More surprisingly, the extents of both loops are very similar, highlighting that there is very weak effect of the emergent rigid vegetation on the transient features of the flows. flow discharge, local streamwise velocity and flow depth are found to be converged after 72 runs of the same hydrograph.  Over the rough bed with emergent rigid vegetation (configuration W) and without (configuration M), the vertical profiles of local mean velocity and Reynolds shear stress (when ensemble averaged) for the base flows of the hydrographs and the steady uniform flow of same discharge are found to be superimposed, as expected.  At a fixed intermediate flow, accelerated (during the rising limb) and decelerated (during the falling limb) mean velocity profiles and Reynolds shear stress profiles are deviating from the steady-uniform flow case conditions. Higher values are found during the rising limb than during the falling limb for both geometries (W and M) studied.  The velocity profiles and Reynolds shear stress profiles at the peak flow of the hydrograph are found to be very similar to the ones for a steady uniform flow cases for both geometries W and M.  The hysteresis loops in the depth-averaged velocity / flow depth relationships for the two geometries have the same shape and dimension, highlighting that rigid emergent vegetation has no damping effect on the non-stationary features of the studied flows.

Figure 1 .

 1 Figure 1. Schematic top view and cross-sectional view of the open-channel (using 1/3 of the total flume width), INRAE, Lyon

First

  , both instantaneous values of water level and velocity are filtered. Water level measurements are despiked. The ADV velocity data are filtered based on the values of the signal to noise ratio (SNR) and signal correlation (COR). The velocity data with a SNR lower than 22dB and a COR lower than 85% are removed from the time series. To account for a small misalignment of the ADV probe with the longitudinal direction, corrections of transverse velocity data are carried out based on[START_REF] Peltier | Estimation of the error on the mean velocity and on the Reynolds stress due to a misoriented ADV probe in the horizontal plane: Case of experiments in a compound open-channel[END_REF].

Figure 4 :

 4 Figure 4: (a) Times series of inlet flow rate Q during 24 hydrograph runs. (b) Time series of Q during two consecutive runs, highlighting the different time periods. (c) Ensemble average of Q (t) for the 109 runs [Q is normalized by base flow Q b ].

Figure 5 :

 5 Figure 5: Convergence of the ensemble averages of discharge Q, flow depth H, and local velocity U x at elevation of z = 24mm from bed roughness. The dashed lines bound measurement accuracy.

Figure 6 :Figure 7 :

 67 Figure 6: Vertical distribution of streamwise velocity U x for a steady uniform flow and at base flows of 7 l/s for (a) a uniform Meadow (M) and (b) a uniform Woodland placed over the meadow (W), at a fixed flow of 15 l/s for rising limb and falling limb, and steady uniform flow of 15 l/s for (c) a uniform Meadow (M) and (d) a uniform Woodland placed over the meadow (W) at the peak of the hydrograph (22 l/s) and for a steady uniform flow of 22 l/s (e) Uniform Meadow (M) (f) Uniform Wood (W)

Figure 8 :

 8 Figure 8: Hysteresis loops of depth averaged velocity U xd against flow depth H for a channel bed with a uniform meadow (M) and a uniform woodland placed over the meadow (W)



  Authors wish to thank CEFIPRA/IFCPAR for sanctioning Raman-Charpak Fellow 2018-2019, Sanction No: SANCTION NO. FT/092/18-19, Reference No. IFC/RCF/2019/1754 to the first author to conduct experimental research on "unsteady flow in rough open channels" at Hydraulics and Hydro-morphology laboratory (HHLab.) of the RiverLy research unit of INRAE (Formerly IRSTEA), Centre de Lyon-Villeurbanne, France.

  J. R. Khuntia Ph.D. Scholar, Department of Civil Engineering, National Institute of Technology Rourkela, India & Visiting research scholar, Riverly, INRAE, Centre de Lyon-Villeurbanne, Villeurbanne, France