Profiling the ToxCast Library With a Pluripotent Human (H9) Stem Cell Line-Based Biomarker Assay for Developmental Toxicity

ToxCast chemical library.

EPA’s ToxCast chemical library has been constructed iteratively using criteria including chemical nomination and procurement, dimethyl sulfoxide (DMSO) solubility, and suitability for testing in automated or semiautomated systems. Drivers for procurement included availability of animal toxicity data, mechanistic knowledge to support model development predicting toxicity, and chemicals of heightened regulatory concern for which data are lacking. The ToxCast chemical inventory file is available at the following link: (ftp://newftp.epa.gov/comptox/Sustainable_Chemistry_Data/Chemistry_Dashboard/2018/September/,last accessed February 6, 2020). For a detailed description of the library, see Richard et al. (2016).

The present evaluation addresses the phase I and II ToxCast library that include pesticides accompanied by guideline animal studies (OCSPP 870 series, and some NTP studies that are guide-linelike), data-poor industrial chemicals, and over a hundred pharmaceutical compounds. This list has 1065 unique structures and 13 duplicates for a total of 1078 samples tested here. Chemical compounds were commercially procured, diluted in DMSO to a stock concentration of up to 100mM (approximately 30% of the chemicals were provided at concentrations lower than 100mM), and plated by Evotec (US), Inc (Watertown, Massachusetts). Aliquots from the stock plates were first diluted with 100% DMSO (Sigma-Aldrich, St Louis, Missouri) to a concentration 1000 times the highest test concentration (HTC) (if necessary) and then diluted 1:1000 in the cell culture media for testing. The final concentration of 0.1% DMSO was a major determinant of the HTC because DMSO itself decreases the ORN:CYSS ratio at concentrations in hESCs above 0.2% (Palmer, unpublished data) and adversely impacts mESCs at concentrations > 0.25% (Adler et al., 2006). We coded stock plates to blind the assay provider to chemical identity. For 7 chemical samples, neat compound was procured from Evotec to enable retesting at concentrations above the range achievable by stock dilutions.

Pluripotent H9 hESC culture.

H9 cells (NIH code WA09, WiCell Research Institute, Inc, Madison, Wisconsin) were used as approved for federally funded research and selected because of their commercial availability, genetic stability (normal female karyotype), and scientific legacy (hundreds of publications). Derivation and characterization of the H9 cell line was originally reported by Thomson et al. (1998). Cells were handled as described (Palmer et al., 2013). Briefly, cells were maintained under feeder-free conditions with mTeSR1 media (StemCell Technologies, Vancouver, Canada) on Matrigel hESC-Qualified Matrix (Corning, Bedford, Massachusetts) coated 6-well plates. Cultures were incubated at 37°C in a humidified atmosphere of <5% CO₂. Differentiated colonies were removed daily through aspiration to maintain the undifferentiated stem cell population. Differentiation was based on visual inspection; there is typically < 5% differentiation in a culture, and only highly pure undifferentiated H9 cell populations were used for these experiments. Cultures were passaged using Versene (Life Technologies, Grand Island, New York) or ReLeSR (StemCell Technologies) at 85%–90% confluency, karyotyped approximately every 10 passages, and the absence of mycoplasma was routinely confirmed with the MycoAlert Mycoplasma Detection Kit (Lonza, Rockland, Maine).

All treatments were carried out in Matrigel-coated 96-well plates. H9 cells were plated with a seeding density of 100000 cells per well in mTeSR1 medium containing 10μM Y27632 Rho-associated kinase inhibitor (ATCC, Manassas, Virginia) to increase plating efficiency. Y27632 was removed prior to compound addition at 24h after plating. The passage number of H9 cells used over the course of this study ranged from 31 to 48; anything above passage 40 was karyotyped within 10 passages prior to use in the assay.

Chemical exposure.

H9 hESCs maintained in a pluripotent state were exposed to test compound for 72h with chemical replenishment with media replacement every 24h. Cell-conditioned media from the final 24-h treatment period was collected for analysis of the targeted biomarker, and cell viability of the corresponding cell layer was assessed as described below. The chemical library was tested in blinded fashion. Plate design and sample workflow is summarized in Figure 1 for the single-concentration screen (no cell viability measures, n=3) and/or 8-point concentration-response series (with cell viability measure, n=3). Each test plate included 1μM Methotrexate (MTX; Selleck Chemicals, Houston, Texas) as a positive reference, 5 nM MTX as a negative reference, 0.1% DMSO as the neutral (vehicle) control, and sample-level media blanks.

Dosing strategy for the protocol initial screen was guided by the “cytotoxicity point” previously reported across 38 different ToxCast assays (Judson et al., 2016). We initially selected 141 chemicals for concentration-response testing; the remaining 924 chemical samples were tested in one concentration. Of those, 252 were retested in concentration-response to confirm a response or adjust the concentration range. The HTC for the initial screen was, for most samples, set to 1, 10, or 100 μM so as not to exceed the chemical-specific median AC50 cytotoxicity point based on Z-score as defined by Judson et al. (2016) unless otherwise limited by compound availability. Other considerations for setting the HTC in concentration-response evaluation included outcome of the single-concentration screen, relevant data from ToxRefDB prenatal developmental toxicity animal studies, and reference compounds used by other studies on alternative test platforms for developmental toxicity (AugustineRauch et al., 2016; Daston et al., 2014; Genschow et al., 2002; West et al., 2010). In all, 379 chemicals were tested in concentration-response and the remaining 686 were negatives at the initial screen. The latter has 14 chemicals remaining where the HTC was an order of magnitude below the ToxCast lower bounds of the median cytotoxicity burst (LBC) (see Supplementary Table S1 for details).

Biosample processing.

H9 cell-conditioned media from the final 24-h treatment period was collected for analysis of the targeted biomarker and cell viability was measured from the corresponding cell layer. Cell viability was measured using the CellTiter-Fluor assay (Promega, Madison, Wisconsin) based on proteolytic cleavage of a substrate to fluorescent signal proportional to the number of living cells (Niles et al., 2007). The cell viability Relative Fluorescence Unit (RFU) was background corrected and normalized to mean RFU of the neutral control (0.1% DMSO). The collected H9-cell-conditioned media samples were processed for targeted biomarker analysis as described (Palmer et al., 2013). Briefly, spent media samples were deproteinized (40% acetonitrile) and processed for UPLC-HRMS. Data acquisition was performed using 4 separate UPLC-HRMS systems, consisting of an Agilent 1290 Infinity LC system (Agilent Technologies) interfaced with an Agilent high-resolution mass spectrometer (models G6520A, G6520B, G6530A, and G6224A). A Waters Acquity UPLC BEH Amide column (2.1mm × 50mm, 1.7lm particle size; Waters, Milford, Massachusetts) maintained at 40C was applied for separation of metabolites using a 6.5min solvent gradient with 0.1% formic acid in water and 0.1% formic acid in acetonitrile (1.0ml/min flow rate). Data were acquired using MassHunter Acquisition software (version B 04.00, Agilent Technologies). The extracted ion chromatogram (EIC) areas for ORN, CYSS, and their respective spike-in C¹³-standards were determined using the Agilent MassHunter Quantitative Analysis software, version B.05.00 or newer (Agilent Technologies), and data were normalized as described in Palmer et al. (2017).

ToxCast annotations.

All raw data and metadata were loaded into a central database for ToxCast using standard nomenclature to pipeline the data as outlined by Filer et al. (2017) into invitrodb (v3 pending release, March 2020) under for the Stemina DevTox^qP assay source identifier (asid) for the ToxCast platform, designated STM_H9 (asid 14); specifically, the 2 assay identifiers (aid) STM_H9_secretome (aid 428) and STM_H9_viability (aid 437)representing data measures for the conditioned media and cell monolayer, respectively. This included identifiers for tracking each chemical sample (spid), plate identifier (apid), well position (row, column), micromolar concentration tested (dose), well quality (0=fail, 1=pass), and well type (wllt). The well type identifiers included media blank “b” lacking H9 cells, neutral control “n” 0.1% DMSO, test compound “t,” negative control “m” 0.005μM MTX, and positive control “p” 1.0μM MTX. Invitrodb_v3 stores raw data values (rval) for the following assay components (acid): peak area for C¹³-cystine (acid 1023) and C¹³-ornithine (acid 1024) tracers, peak area for measured cystine (acid 1025) and cystine standardized to the C¹³ tracer (acid 1026), cystine normalized to the plate median value of the neutral controls (acid 1027), peak area for measured ornithine (acid 1028) and ornithine standardized to the C¹³ tracer (acid 1029), ornithine normalized to the plate median value of the neutral controls (acid 1030), the ORN:CYSS ratio calculated from DMSO-normalized values (acid 1031), the targeted biomarker prediction (acid 1032, which is an empty placeholder), background-corrected RFU from the CellTiter-Fluor assay (acid 1113), and cell viability normalized to mean RFU of the DMSO control (acid 1114). Virtual plate maps reconstructed from the ORN/CYSS ratio (Figure 1) were visually inspected to confirm consistency in each well before entering the data into invitrodb_v3.

The corresponding assay endpoint identifiers (aeids) analyzed in the Results presented here address the 4 main assay component features for the predictive model: the decrease in media ornithine reflecting reduced cellular release (STM_H9_ornithineISnorm_perc_dn, aeid 1689), the increase in media cystine reflecting less utilization (STM_H9_cystineISnorm_perc_up, aeid 1682), the ORN/CYSS ratio reflecting a decrease in the ORN:CYSS ratio as the primary biomarker (STM_H9_OrnCyssISnorm_ratio_dn, aeid 1693), and normalized cell viability (STM_H9_Viability_norm, aeid 1858). The processed STM dataset described and analyzed here comprised 79398 individual data points across 1065 unique chemical structures.

ToxCast Data Pipeline (tcpl).

Individual data points from the 379-chemical concentration-response series were processed through the ToxCast Data Pipeline (tcpl) (Filer et al., 2017). Level 0 of tcpl is the entry point for rval from the 4 key assay component features. All individual data points for the same chemical were concatenated at this level and processed through 6 levels of data processing. Level 1 flagged the DMSO-normalized samples for plate position effects, any missing replicate(s), or samples with poor well quality. The latter was found for individual wells across 30 chemicals where the highest concentration caused significant loss of cell viability and no measurable ORN. Generally, the few invalid samples reflected extreme cytotoxicity at higher concentrations. All are included in the concentration-response assessment, and a teratogenic index could still be computed if lower concentrations remained in the noise belt. If not, then we retested the samples at a lower concentration range and all of the data would be included in the tcpl profile. For plotting purposes, we assigned a minimal recorded ORN value (0.001) that was well below the limits of detection of the metabolite (0.01). This was applied only to the 30 data points having poor well quality. Level 2 transformed the individual replicate data points into an appropriate “response unit” computed on a log2 scale. Level 3, typically a normalization step in tcpl, applied no additional normalization but instead inverted the log2 data to graph response profiles in a manner consistent with other ToxCast assay platforms, ie, responses increase from a baseline of zero. At this level, tcpl calculated a baseline median absolute deviation (bmad) of the response variable using only the lowest 2 test compound concentrations from all test wells, concentrations where we expect no activity for the vast majority of chemicals. Note that the 731 chemicals with STM data calls based solely on single-concentration screen data (eg, that had not been tested to date in a definitive concentration series) are pipelined to level 2, with single-concentration hit calls summarized based on a global threshold of 3*bmad across all tested compounds, but not plate-level controls. For the multiconcentration series, the same 3*bmad was implemented as a noise belt distribution in level 4 that calculated the parameters for automated curve-fitting models. Three curve-fitting models were applied (Constant, Hill, and Gain-Loss) and the winning model, selected by Akaike Information Criterion (AIC), progressed to level 5 for graphing. Level 6 then applied warning flags for curve-fitting issues or data quality concerns such as local plate effects, single point hits, and noisy data. Chemicals with a measured response on the ORN/CYSS ratio greater than 3*bmad were classified as “active” at the concentrations tested here. Because “inactive” results (hitc = 0) are not stored in the database, these were replaced by an arbitrary value of 10⁶ μM post-tcpl for computing purposes, consistent with other assays in the ToxCast database.

Developmental toxicity (in vivo) correlation analysis.

To assess STM model performance against in vivo developmental toxicity, we identified an anchoring set of 42 benchmark compounds from the ToxCast phase I/II inventory. An initial list was compiled having overlap with previous literature studies that aimed to evaluate alternative methods for developmental toxicity (Augustine-Rauch et al., 2016; Daston et al., 2014; Genschow et al., 2002; West et al., 2010; Wise, 2016). Note that 18 of the 42 compounds were in the training set and 2 of the compounds were in the test set from the original devTOX publication (Palmer et al., 2013). From this list, we selected those chemicals having traditional (pre-2015) US FDA (Food and Drug Administration) categories for potential risk to the developing fetus if pregnant women are exposed (or in a few cases, well-defined developmental toxicants from the National Toxicology Program). A final set of positives (n=26) and negatives (n=16) was developed using evidence from teratology cohorts tested under EPA Health Effects Test Guidelines OPPTS 870.3700 (https://ntp.niehs.nih.gov/testing/types/devrepro/index.html, last accessed February 6, 2020) (now OCSPP 870.3700 even though it is has not been corrected on regulations.gov). STM model performance was computed from 2 × 2 binary classification tables of true positive (TP), false positive (FP), false negative (FN), and true negative (TN) conditions (Powers, 2011). Overall (regular) accuracy (Rand ACC) was computed from the TP rates: sensitivity = TP/(TP + FN), specificity = TN/(TN + FP), and Rand ACC = (TP + TN)/(TP + FP + FN + TN). Balanced accuracy (BAC) was computed as the average of PPV = TP/(TP + FP) and negative predictive value (NPV) = TN/(TN + FN); BAC = (PPV + NPV)/ 2. Model performance was assessed with Matthews correlation coefficient (Matthews, 1975). Because the testing protocol was independent of the original training by Palmer et al. (2013), the subset of 20 compounds overlapping in the 42 benchmark set that were used to develop the model (Palmer et al., 2013) are not used here for training chemicals; rather, they confer confidence in the testing strategy because the contractor was blinded to their identity.

To evaluate a broader correlation of STM models to in vivo animal studies, we used prenatal developmental toxicity study type (870.3700) downloaded from the ToxRefDB v1.0 data release (https://www.epa.gov/chemical-research/toxcast-toxrefdb-datarelease, October 16, 2015, accessed July 25, 2017). This included chemical-endpoint-specific dosage for the Lowest Effect Level (LEL) in mg/kg/day if the endpoint was tested (or assumed to have been tested) and a treatment-related effect was observed (Knudsen et al., 2009). Lowest Effect Level values were culled from the “Endpoint_summary” download file. The “study-endpoint-category” calls reflected developmental (dLEL) and maternal (mLEL) endpoints for each available study, else “NULL” meaning no LEL was observed at the highest dosage tested (assigned an LEL of 10⁶ mg/kg/day for computing purposes). For the 1049 rat and rabbit studies, 924 ToxRefDB study records (approximately 400 chemicals) were oral route, 60 records (16 chemicals) inhalation, 48 records (16 chemicals) direct, 42 records (16 chemicals) dermal, and 4 records (2 chemicals) “other.” “Direct administration” in ToxRefDB refers to routes of exposure other than oral, dermal, or inhalation (eg, intraperitoneal, intramuscular, and intravenous). Bioavailability might be questionable for the dermal or inhalation study records with a negative outcome for developmental toxicity in ToxRefDB. This situation applied to 35 of 1049 records, of which 30 records correlated with a negative response (TN) in the STM response and only 5 records correlated with a positive response (FP). There may be overlap wherein some chemicals may have records on alternate routes of exposure. The list of chemicals is provided in Supplementary Table S2. For the few inhalation studies having ppm exposure unit, we used a 1:1 mass conversion rather than assuming a particular breathing rate to get at this volume. In all, 3496 records were obtained for dLEL, mLEL. or NULL entries. These records collapsed into 791 outcomes (424 rat, 331 rabbit, 33 mouse, 2 hamster, and 1 dog). Given the preponderance of rat and rabbit studies, we focused only outcomes from those 2 species to build the endpoint classifier model (400 rat OR rabbit, 389 rat, 323 rabbit, and 203 rat AND rabbit). This removed outcomes reported in mouse/hamster/dog from consideration, even if they may have been reported in rat/rabbit as well.

We first assigned levels of evidence for developmental toxicity for each individual ToxRefDB record as follows: “clear evidence” (dLEL ≤ 200mg/kg/day; dLEL < mLEL), “some evidence” (dLEL ≤ 200mg/kg/day; dLEL ≤ mLEL), “equivocal evidence” (dLEL < 10³ mg/kg/day; dLEL > mLEL), and “no evidence” demonstrated by data from a study with appropriate experimental design but no developmental effects observed (eg, no dLEL or dLEL ≥ 1000mg/kg/day) (see Supplementary Table S2). This cutoff was arbitrary, selected from convenience to match the arbitrary in vitro cutoff of 200μM for the Teratogenicity Index (TI) reflecting the critical change in the ORN/CYSS biomarker. Replicate studies conducted with the same chemical compound and species were collapsed into a single dLEL or mLEL value with the evidence hierarchy: clear > no > some > equivocal. Essentially, “equivocal evidence” indicates that a dLEL was determined for the compound, in a rat or rabbit study, but there is not enough evidence to attribute the observed fetal endpoint to a specific developmental effect as opposed to a maternal effect. It is important that “evidence” not be equated to “importance.” A maternally mediated adverse developmental outcome is important, but perhaps falls outside the biological domain of the assay. The final number of compounds for in vivo anchoring was 401 from ToxRefDB and completed to 432 by adding information from the 11-overlapping compounds in the 42 benchmark compounds identified in the literature (described above) but not compiled in ToxRefDB.

To produce a call TP, FP, FN, and TN for each chemical tested, we looked for rat-rabbit concordance where available. This was done by discordance (OR = hit in either species constituted a positive) and concordance (AND = same calls from both species). Stringency models built onto the BM-42 calls as follows. Base model: defines any chemical with a dLEL as positive in either species, else negative, and the STM data accepts as positive any dTP < 1000; n=432. Low stringency model: defines a positive as CLEAR or SOME evidence, else negative (EQUIVOCAL, NO); the STM-positive data are filtered at an arbitrary cutoff of TI <= 200 μM; n=432. Medium stringency: defines a positive when either species (rat OR rabbit) shows CLEAR evidence, else negative; n=285. High stringency: defines positive as a concordant response (rat AND rabbit) showing CLEAR or NO evidence in both cases; n=127. Note that the BM-42 set was defined as CLEAR and NO; there are 11 ToxRefDB among the 42 tested in rat of which 4 were also tested in rabbit.

Biochemical (in vitro) correlation analysis.

To evaluate the biological perturbations underlying the targeted biomarker in the STM assay, we utilized the HTS results from the ToxCast NVS cell-free biochemical-assay platform (Knudsen et al., 2011; Sipes et al., 2013). The information for this analysis is provided in Supplementary Table S3. AC50 (half-maximal activity concentration) values for NVS aeids were selected from invitrodb_v3 hit-call matrix (August 10, 2018 release, accessed September 18, 2018) (https://www.epa.gov/chemical-research/exploring-toxcast-data-downloadable-data). Note that the public download has AC50s for chemical-assay relationships that are flagged as inactive because they do not meet the assay-specific efficacy threshold for bioactivity; the hit-call file converts them to inactive and inserts log10=3 for all inactive AC50s. The NVS dataset had AC50 values for activation or inhibition across 337 biochemical assays for substrate-product enzymatic activity or ligand-receptor binding, corresponding to 420 total features, due to the possibility of enzymatic activation/inhibition or nuclear receptor (NR) agonist/antagonist binding activity. This produced 2918 hit calls in the 420-feature × 1065 chemical matrix expressed as log10 micromolar AC50.

The NVS AC50 matrix was further consolidated to a “gene-score” feature matrix using chemical-specific Z-scores from the global ToxCast cytotoxicity burst (Judson et al., 2016) and the official gene symbol for the human protein function that is annotated in each particular ToxCast NVS assay. This filtered-out any chemical-assay pair having an AC50 value within 3 standard deviations of the cytotoxicity burst (Leung et al., 2016). The resulting gene-score feature matrix, G(chemical, gene_dir), comprised the mean log₁₀(AC50) for each homologous gene symbol and was assigned an “up” or “down” extension to account for directionality (dir) in terms of enzymatic activator/inhibitor or receptor ligand binding activities. Finally, this matrix was transformed to micromolar concentration potency terms using the following equation:

where G is the logAC50 gene-score matrix and GPS is the transformed gene-potency score matrix (1062 chemicals × 267 NVS gene features) in the analysis. Binary calls from the ToxCast STM dataset (1=active STM response, 0=inactive STM response) were then used as the categorical endpoint to fit a multiple logistic regression model (Scikit-learn, v. 0.18.1; Python, v.2.7.13):

where w_i represents the log-odds coefficient term for each gene. These weights can be positive or negative to represent how strongly the gene association correlated to an active or inactive STM response, respectively.

To characterize the phenotypic importance for each gene in the dataset, we used the Human-Mouse: Disease Connection (HMDC) database (http://www.informatics.jax.org/humanDiseases.shtml, accessed July 31, 2018). Although rat is the preferred species in tiered approach for developmental toxicity testing, performing evaluations in the rat and mouse or mouse and rabbit detected 80% of the 105 teratogens examined in the veterinary pharmaceutical products study by Hurtt et al. (2003). As such, the HMDC database provides a reasonable surrogate for modeling rat-human disease connections. The phenotypic weighting was conducted independently of the logistic regression described above and serves as a way to determine the biological relevance of each gene present in the dataset. HMDC integrates data on phenotype and disease model data from the MGI Mouse Phenotype Ontology (MP) with human gene-to-disease relationships from the National Center for Biotechnology Information, Online Mendelian Inheritance in Man (OMIM), and the Human Phenotype Ontology. We mapped all NVS gene symbols to the 28 available HDMC phenotype systems (extracted from the HDMC database) and assigned a bin score (b) based on evidence weighting levels from the number of HMDC entries (n) for any (n=1, b=1), some (n=2–6, b=2), or many (n>6, b=3) annotations, or no HMDC data (n=0, b=0). We then calculated a gene-specific MP weight by summing the bin score across all 28 MP categories for the 233 HMDC gene symbols that could be linked to a NVS feature. This list was standardized from 0 to the maximum value and weighted each NVS target gene by relevance to a human-mouse disease system. It should be noted that each phenotype system has subordinates and not all of which may directly apply to developmental toxicity (the deeper level of ontology was not considered here). Also note that a few NVS gene entries lacked an HMDC record: the normalized MP weighting schema does not filter out those genes but scales them at the low end of the evidence bar for a phenotypic system. The top and bottom 40 from the NVS gene targets list linked to the STM-positive and STM-negative responses, respectively, were annotated separately using the Functional Annotation Tool from the DAVID Bioinformatics Resources 6.8, NIAID/NIH (https://david.ncifcrf.gov/). Categories were selected for a minimum of 4 genes and maximum 30 genes from GO Direct (molecular function, cellular compartment, and biological process), OMIM, KEGG, BioCarta, Reactome, and INTERPRO at a Bonferroni-adjusted p-value ≤ .05 per category. Nothing from OMIM or BioCarta passed. All other redundancies were resolved manually by the lowest false discovery rate value (as these have lesser uncertainties across the specific annotation) or in some cases the most informative annotation record. The resulting annotation records were visualized in a Spearman correlation matrix based on summed weights of represented gene targets from the logistic regression model, colored by STM on one axis and developmental toxicity on the other.

Determination of STM Positivity

The bmad determined in tcpl is calculated as the median absolute deviation of the response variable from the lowest 2 test compound concentrations in the assay between tcpl levels 3 and 4 (Filer et al., 2017). Because of the way tcpl computes the noise belt, a point of departure represents the state of the platform in ToxCast. This threshold is dynamic and could change as more and more compounds are tested in the future. As such, the ORN/CYSS ratio = 0.76 represents the biomarker’s point of departure threshold (eg, teratogenicity Index, TI) in the current state of 1065 cases and is not a recommendation for general applicability. Figure 2 plots the distribution of tcpl data at level 3 for plate-level neutral controls (n=1176) and Methotrexate (MTX) reference for positive (1μM, n=589) and negative (5nM, n=590) response; and at level 4 for the lowest 2 test concentrations of ToxCast samples (n=2169). A threshold drawn at 3*bmad for each feature distinguished all MTX-positive and MTX-negative plate references. CYSS values showed wider variability than ORN values due to a minor change in the culture medium formulation that increased CYSS utilization over the course of the contracting period. Aside from a few outliers that were not removed when calculating bmad, the lowest 2 concentrations from all ToxCast test wells mimicked the neutral control and MTX-negative reference groups.

Feature plots for each chemical tested in concentration-response (n=379) are given in Supplementary Figures S1a–S1f. As an example, Figure 3 displays results for Methotrexate as part of the ToxCast library. The gray hatched zone indicates the global noise belt computed as a 3*bmad for each main feature. Due to the way tcpl calculates bmad, the threshold ORN/CYSS ratio (0.88) originally described (Palmer et al., 2013) fell within the noise belt. Consistent with other ToxCast assays, we used a statistical threshold for positivity that reflects the test concentration at which the activity departs from the noise belt (acb, activity concentration at baseline). This is indicated by the vertical yellow line in Figure 3 and for the ORN/CYSS ratio gives a TI that reflects the concentration threshold predicted for human developmental toxicity. Because several outliers widened the global noise belt this, in turn, resulted in a right-shift in the TI relative to the default ORN/CYSS model from the original assay description (Palmer et al., 2013). For example, 3*bmad (log2) was 0.403, rendering the acb (log2^−0.403) = 0.756 versus 0.88 in the default model. For normalized cell viability (percent control), the 3*bmad = 0.161 rendering the acb (log2^−0.161) = 0.894. This threshold corresponds to a global decrease of 11% cell viability, which is taken as a positive effect for this feature but may or may not reflect overt cytotoxicity. Although minor effects on cell viability could account for changes in cellular ORN release and/or CYSS uptake measured in the ORN/CYSS ratio, altered cell growth and/or survival are potential modes of action in developmental toxicity. As such, compounds that impact the ORN/CYSS biomarker due to minimal effects on cell viability should not be discounted because of it. The blinded ToxCast Methotrexate sample registered a TI = 0.059μM and 11% loss of cell viability at 0.062μM, consistent with the MTX plate references.

To examine the impact of the right-shift in the tcpl-derived critical concentration for positivity, we next plotted the difference in hit calls between the original value (0.88) (Palmer et al., 2013). We binned each chemical by measured maximum response into a distance from bmad. For example, chemicals in the first bin had a maximal response within 1 bmad, chemicals in the second bin had a maximal response within 2 bmads, and so forth. Quantitatively, the tcpl model shifted active calls to a slightly higher concentration versus the published default that fell within the noise belt. This shift, however, had little impact resulting in the loss of only 2 chemicals from the active list (Maneb and 2-methoxy-5-nitroaniline) that displayed marginal efficacy.

Profiling the ToxCast Phase I/II Inventory

Level 2 and level 6 ToxCast STM data are contained in Supplementary Table S1. Because of the way ToxCast assays were contracted, multiwell stock plates (100μM where possible) permitted testing most chemicals at HTC = 100μM (due to potential interference with DMSO). Variations in the HTC was in some cases limited to < 100μM or < 10μM based on a low “median cytotoxicity burst” across dozen cytotoxicity and cell stress assays in the ToxCast portfolio (Judson et al., 2016). On the other hand, tcpl’s automated curve-fitting extrapolated TI above the HTC of 100μM in several cases. Across the 1065 chemical samples tested, we observed a critical decrease in the ORN/CYSS ratio for 208; however, if we consider the response positive only if TI < 200μM as reasonable cutoff, then this condition was met by 205 (19%) of the ToxCast phase I/II library. Concentration-response profiles are shown for ORN, CYSS, the ORN/CYSS ratio, and cell viability plots in Supplementary Figures S1a–S1f. We retested all compounds with activity in the single-concentration screen in concentration-response aside from tetracycline and 3,3’-dimethylbenzidine, which produced a weak effect on the ORN/CYSS ratio at 100μM but were included among the positives.

Concentration-dependent drops in the ORN/CYSS ratio were driven by ORN alone for 15 chemicals (7.3% of the positives, including thalidomide), by CYSS alone for 36 chemicals (17.6% of the positives), and both metabolites in 147 cases (71.7% of the positives). A few chemicals (n=7) produced an overall effect on the ORN/CYSS ratio without scoring a hit on either metabolite alone due to noisy data flagged for individual metabolites. The reduction in CYSS utilization seemed to drive the ORN/CYSS response for most of the positive chemicals, minor variations were observed for a few compounds (see Supplementary Figures S1a and S1d). All-trans-Retinoic acid, for example, increased ORN production (an effect shown for retinoids as a class through inhibition of ornithine decarboxylase [Palmer et al., 2017]), but the active response was driven by a stronger reduction in CYSS utilization, indicating that different mechanisms may lead to changes in the ORN/CYSS ratio. For more complex concentration-response behavior, octyl gallate produced a transient spikein CYSS utilization (acb = 0.39μM) prior to concentrations that decreased CYSS utilization (acb = 1.83μM) (Supplementary Figs. S1a–f). The nonmonotonic changes in CYSS may foreshadow cytotoxicity but perhaps not as a general phenomenon.

The general relationship between TI (Supplementary Figure S1e) and H9 cell viability (Supplementary Figure S1f) is shown graphically for the 205 positives confirmed in the concentration-response series (Figure 4). Although cell viability is measured, it is not included in the prediction by the assay, which is based solely on the ORN/CYSS ratio response. Cell viability is provided as an additional endpoint to aid in the interpretation of the data. The 11% decline in viable H9 cell number (point of departure from the noise belt) could be due to impaired growth, reduced viability, or both. For example, 2 well-recognized human teratogens, Thalidomide (Figure 4, No. 8) and Methotrexate (Figure 4, No. 154), invoke a biomarker response at concentrations differently accompanied by changes in viable cell number. It should be noted that the highest concentration tested here approached within an order of magnitude or exceeded the lower bounds of the ToxCast cytotoxicity point (LBC) for 204 of 205 STM-positive compounds (99.5%) and 843 of 883 STM-negative compounds (95.3%). Approximately one-third of the STM-positive compounds triggered a critical drop in the ORN/CYSS ratio without altering cell viability at the concentrations tested here (sector A in Figure 4 and Table 1). This is exemplified by all-trans-retinoic acid and thalidomide (TI = 0.003 and 1.27μM, respectively). In the remaining two-thirds of the STM-positive compounds, an effect on the ORN/CYSS ratio cannot be generally dissociated from 11% drop in cell viability. Methotrexate and cytarabine, for example, dropped the ORN/ CYSS ratio (TI = 0.059 and 0.054μM, respectively) as well as cell viability (acb = 0.062 and 0.082μM, respectively).

Benchmarking STM Assay Performance for Teratogenicity

To evaluate ToxCast STM assay performance, we compared the in vitro classification with a set of 42 benchmark compounds that have been used by others to evaluate developmental toxicity alternatives (Augustine-Rauch et al., 2016; Daston et al., 2014; Genschow et al., 2002; West et al., 2010; Wise, 2016) (Table 2). Most of these compounds have information from FDA labels on potential risk to the developing fetus if pregnant women are exposed, ranging from safe for use during pregnancy (category A) to contraindicated during pregnancy (category X). Although the FDA labels are no longer used (since 2015), they have been used in evaluating alternatives to animal testing for developmental toxicity and, as such the category descriptors are given for convenience in parentheses as follows. (FDA category A where adequate and well-controlled studies have failed to demonstrate a risk to the fetus in the first trimester of pregnancy [and there is no evidence of risk in later trimesters]; category B, either animal-reproduction studies have not demonstrated a fetal risk but there are no controlled studies in pregnant women, or animal-reproduction studies have shown an adverse effect [other than a decrease in fertility] that was not confirmed in controlled studies in women in the first trimester [and there is no evidence of a risk in later trimesters]; category C where animal-reproduction studies have shown an adverse effect on the fetus and there are no adequate and well-controlled studies in humans, but potential benefits may warrant use of the drug in pregnant women despite potential risks; category D where there is positive evidence of human fetal risk based on adverse reaction data from investigational or marketing experience or studies in humans, but potential benefits may warrant use of the drug in pregnant women despite potential risks; and category X where studies in animals or human beings have demonstrated fetal abnormalities, or there is evidence of fetal risk based on human experience, or both, and the risk of the use of the drug in pregnant women clearly outweighs any possible benefit. The drug is contraindicated in women who are or may become pregnant.) The ToxCast STM assay based on defining a positive using the tcpl threshold 0.76 performed with an accuracy of Rand ACC = 78.6% (sensitivity 0.65, specificity 1.00, n=42) and BAC = 82.0% (PPV 1.00, NPV 0.64, n=42). This was consistent with the pharmaceutical-trained model from the assay provider (77% accuracy, 0.57 sensitivity, 1.00 specificity, n=23) (Palmer et al., 2013). The targeted biomarker outperformed 11% reduction in H9 cell viability as a predictor of teratogenicity (54% accuracy, 0.28 sensitivity, 0.94 specificity, n=42).

Ten compounds in ToxCast phase I/II are found among 28 compounds on the Daston List (DL) for exposure-based comparisons of developmental toxicity with regards to Cmax for maternal plasma toxicokinetic studies at a dosage producing, or not producing, teratogenicity/embryolethality (Daston et al., 2014). We evaluated concordance of the STM model for these 10 compounds across 7 negative and 7 positive exposure-based dosimetry calls from the Daston et al. (2014) study. In 7 cases, the highest concentration exceeded what could be achieved by diluting from 100mM stock plates at the final DMSO concentration. Those were tested from neat compound in solid or solution form. Results are shown in Table 3. Teratogenicity index correctly called 6 of 7 DL-negatives and 5 of 7 DL-positives, again yielding 78.6% concordance. The 3 misses were 1,2-propylene glycol (FP), caffeine (FN), and ethylene glycol (FN). Evaluation of 1,2-propylene glycol yielded a TI of 246664μM although the result is not considered reliable given the exceedingly high (1M) test concentration required to achieve 850000μM for assessing DL-negativity (Daston et al., 2014). Daston List concordance improves to 84.6% if this FP call is discarded. Caffeine is negative in animal studies when consideration is given to appropriate control data (Wise, 2016). Model performance further improves to 91.7% if this FN call is disregarded. No explanation can be offered for missing ethylene glycol aside from the lack of in vitro bioactivation; however, the ToxRefDB result for this chemical is equivocal with regards to prenatal developmental toxicity in rats and rabbits. As such, the exposure-based DL model performed the same (78.6%) or better (> 84.6%) than the 42-benchmark compound model.

Evaluating STM Assay Performance Against ToxRefDB Prenatal Studies

To evaluate a broader concordance of the ToxCast STM assay to in vivo animal studies, we culled the data from endpoint summary files in ToxRefDB prenatal studies that provide no effect level and LEL (Knudsen et al., 2009) and manually collapsed these data into evidence-based calls for developmental toxicity in rat and/or rabbit studies (see Materials and Methods for details). Each study attempts to achieve a maternal (mLEL) and fetal (dLEL) dosage for developmental parameters observed at term. Note here that some endpoint categories (eg, placenta, resorptions, and postimplantation loss) do not strictly map to dLEL outcomes and are thus not reflected in this analysis. Note that information from Table 2 was carried across the performance models, although only 11 of 42 benchmark compounds had corresponding ToxRefDB data. Table 4 summarizes the results from 2 × 2 contingency models, and Supplementary Table S2 provides the data and models used for this analysis. In all, we identified 432 chemicals with STM data and endpoint effects data.

Fundamentally, ToxCast STM performed weakly against an unfiltered compilation of developmental toxicity studies that included any dLEL call in a rat or rabbit (Rand ACC = 47.2%, 0.33 sensitivity, 0.82 specificity, n=432). Concordance to the in vivo classification increased when the criteria for evidence of developmental toxicity became more stringent. Rand ACC increased to 60.4% (n=432), 71.9% (n=285), and 79.5% (n=127) for low, medium, and high stringency anchors, respectively (Table 4). The latter sets a true condition as dLEL < mLEL and dLEL ≤ 200mg/kg/day, and no developmental toxicity as ≥ 1000mg/kg/ day concordantly in rat and rabbit studies. Therefore, ToxCast STM accuracy reaches 78.6% (Rand ACC) to 82.0% (BAC) when there is high confidence in the call for developmental toxicity.

Biochemical Stratification of the STM Response

We next sought to correlate the STM-positive and STM-negative responses to biochemical profiles from the ToxCast NVS dataset. This dataset was selected because it reflects potential molecular initiating events (MIEs) in chemical-target interaction, and within the different STM-positive and STM-negative compounds there are compounds with known MIEs that could serve as indicators for how well the correlation is picking up MIEs. A binary classification outcome model of the ORN/CYSS ratio was built with a logistic regression strategy using assay-specific AC50s from 420 NVS features. Figure 5 shows the workflow for this operation and Supplementary Table S3 provides the corresponding data.

The NVS chemical-assay matrix started with 2918 chemical-biochemical features in the NVS dataset connected to a ToxCast gene score, whereby assays were aggregated for the gene symbol most closely linked to a biochemical target where the cell-free AC50s fell 3*mad below the ToxCast cytotoxicity burst (Judson et al., 2016; Leung et al., 2016). For assay selection, logistic regression identified 267 molecular targets having more statistical concordance over the different STM-positive or STM-negative compounds. The algorithm tries to maximize the distance between weak and strong potency, resulting in gene-potency scores (GPS) typically in the lower nanomolar range (eg, < 75nM). Supplementary Table S3 lists them in rank order. Perturbed ligand binding to NR3C1 (glucocorticoid receptor) and ESR1 (estrogen receptor) represented the top GPS values for the STM-positive and STM-negative responses, respectively.

For further gene reduction, we next weighted the 267-gene list by evidence for gene-to-disease connections across 28 phenotypic systems in the HMDC database. This rearranged the gene list by evidence for gene-to-disease associations from curated literature. Supplementary Table S3 gives the adjusted rank order, and Figure 5 shows the top 10 genes surfacing to the STM-positive and STM-negative responses, respectively, mapped to canonical pathways and biological processes. For example, receptor tyrosine kinases (RTKs) and nonreceptor tyrosine phosphatases joined NR3C1 glucocorticoid receptor as top STM-positive contenders, whereas several NRs and G-protein coupled linked receptors joined ESR1 as top STM-negative contenders. Further demarcation of potential pathways and processes selected the top 40 phenotype-weighted NVS targets in each domain for functional annotation via NIAID/NIH DAVID Bioinformatics Resources 6.8 (see Materials and Methods). Note that 2 gene products (HDAC6, PTPN9) were present in both up/down sides resulting in 82 annotation records. Best results were obtained when STM-positive and STM-negative lists were annotated separately, utilizing a minimum of 4 genes and maximum 30 genes from GO Direct (molecular function, cellular compartment, and biological process), OMIM, KEGG, BioCarta, Reactome, and INTERPRO annotation systems. Annotation records were selected for Bonferroni-adjusted p-value ≤ .05 per category; redundancies were resolved manually by the lowest false discovery rate value (as these have lesser uncertainties across the specific annotation) or in some cases the most informative annotation record. This resulted in 60 overlapping annotation records (Supplementary Table S3) visualized by Spearman correlation (Figure 6).

The 60 annotation records summed for 68 potential MIEs showed an inconsistent relationship to the medium stringency anchor when colored by STM response on one axis and developmental toxicity on the other (eg, rat or rabbit outcomes, n=285) (Figure 6). We further collapsed the 60 records into several “keystone” pathways that accounted for 75% of the MIEs (51 of 68, Supplementary Table S3). Figure 7 maps 34 MIEs (NVS) that correlated with developmental toxicity and detected (TP) or not (FN) with the STM response. The flow of regulatory pathway information points to AKT/FoxO signaling and focal adhesion as determinants in the applicability domain (RTK signaling), and Ca++ second messenger generation (G-protein coupled receptor [GPCR] signaling) via G(q) pathways and as well NR-mediated gene expression as determinants of developmental toxicity outside the applicability domain. Aside from the glucocorticoid receptor (NR3C1), the pathway model provides the basis for a mechanistic speculation of where the STM response is lacking in sensitivity.