1. Introduction

Neurodevelopmental disorders such as impairment of learning and memory, and cognitive dysfunction, are of serious concern due to the health risks and consequences on the developing brain resulting from exposure to exogenous chemicals [1], [2], [3]. The assessment of developmental neurotoxicity (DNT) is not a mandatory requirement in the European Union or the United States of America and is not routinely conducted. However, it may be undertaken when data from developmental and/or reproductive toxicity studies on adult animals indicate a possible concern for neurotoxicity [4]. When testing is carried out, it is based on available DNT testing guidelines using in vivo methods [5], [6], [7]. These animal tests are a starting point in deciphering complex endpoints such as DNT. However, with the limitations of in vivo testing for DNT [8], there is an opportunity to consider new approach methodologies (NAMs), such as a battery of in vitro test methods, omics technologies and in silico models as a viable alternative. Whilst NAMs for DNT are not yet standardised or required by regulatory authorities, they can provide valuable mechanistic insights regarding potential developmental neurotoxicants [4], [9], [10], [11].

Frameworks are required to organise information from NAMs to predict complex toxicological endpoints, including neurotoxicity and DNT. A promising approach to organise DNT information and subsequently use the information to develop a predictive model is provided by the adverse outcome pathway (AOP) concept [12]. An AOP represents a formal description of a series of events from the molecular initiating event(s) (MIE(s)), key events (KEs) at corresponding molecular, cellular, and tissue levels to adverse outcomes (AOs) at organ, organism and population levels [12]. Recent progress in the development of qualitative and quantitative AOPs underlines their utility to design appropriate experiments and computational simulations [13], [14]. However, as DNT is a complex process with multiple molecular and cellular paths, no single AOP is able to fully explain this complex event. Thus, a network of AOPs allows for a better depiction of the overall mechanistic understanding of DNT than a single AOP. To illustrate this point, given the limited resources available to quantify biological paths of DNT, identification of common key events (CKEs) that intersect the individual paths can assist in the design of testing strategies and in the computational modelling of data-rich biological events [15]. Such CKEs are characterised by high connectivity located, are centrally within the network of AOPs and are essential to link multiple linear AOPs, i.e., MIE(s) to the AO(s) [16], [17].

Quantitative AOPs (qAOPs) and qAOP networks are increasingly being modelled using probabilistic methods [18], [19]. Bayesian modelling is an approach that fits, and makes inferences from, data using Bayes’ theorem for variables. Bayes’ theorem is a mathematical formula that transforms the prior, i.e., our knowledge about the data before seeing the data, into posterior distributions based on the evidence provided [20], [21]. An advantage of this approach, besides the ease of computing predictions and associated uncertainties relating to variables of interest, lies in the reproducibility of the predictions; once the prior beliefs are defined, similar values will be obtained each time the model is run unless new evidence is added [20], [22].

Given the limited availability and heterogeneity of DNT information, i.e., in vivo and in vitro data, and the potential for Bayesian machine learning to investigate chemical-induced DNT, the aim of this investigation was to quantify a simplified, reduced version of the AOP network for neurotoxicity developed by Spînu et al. [23]. The main objective was to predict the probability that a compound induces individual CKEs, in addition to predicting the probability of inducing the AO. The analysis took into account potential correlations and causal relations informed by the key event relationships (KERs) and additional information, such as physicochemical, in silico and in vitro data. This investigation also aimed to explore whether Bayesian hierarchical modelling is fit-for-purpose in chemical risk assessments informed by qAOP models.

2. Materials and methods

2.2. Data description

Two in vitro studies [26], [27] that tested compounds for their DNT potential were chosen as primary data sources. The compounds tested from the two studies were merged into a single list (e.g., nine pairs of compounds, tested as different forms of salts, were combined) and aligned with information referring to the compound name, Chemical Abstracts Service Registry Number (CAS RN), simplified molecular input line entry system (SMILES) string and the US EPA Comptox Chemical Dashboard substance identifier (DTXSID). In total, 88 compounds served as a starting point to collect additional information to improve the modelling. The list contained different types of compounds: pharmaceuticals, pesticides and industrial chemicals. The two in vitro studies differed in the concentration ranges used, type of viability assay applied, plating densities of the cells, calculation of the effective concentration (EC), and the exposure period (five days vs 12 days). Synaptogenesis was measured by a battery of high content imaging and microplate reader-based assays, while the viability was based on the number of cells per field or determined in sister plates using a luminescent assay [27]. Neural network activity was measured on day 12 by the microelectrode array recordings and the viability by two assays, the total lactate dehydrogenase (LDH) release and Alamar blue assay [26]. Also, valproate was tested in two concentration ranges in both studies: lower and higher. The range of higher concentrations showed response/activity in both studies and, thus, the EC_x values of the range of higher concentrations tested, and associated responses, for valproate were retained. The difference in the concentration range and time of effect shows that impacts on synaptogenesis are observed at lower levels than impacts on neural network formation. This concentration–response concordance fits the expected pattern for a causal relationship between the two key events in an AOP [28].

The compounds were assigned as positive or negative for DNT induction based on the in vivo studies summarised and evaluated in a literature review by Mundy et al. [29]. In the latter paper, the authors applied rigorous selection criteria to 408 chemicals, concluding that 97 showed evidence of developmental neurotoxicity. Of these 97 DNT-positive chemicals, 21 had evidence in human studies. The compounds were assigned positive/negative for the reduction of BDNF based on a literature review performed to evaluate the peer-reviewed publications that studied the impact of compounds on this effect (Supplementary material, Table S2). Given that at present there is no standardised protocol for assessing the inhibition of BDNF available and thus, no existing datasets, the main selection criteria of the peer-reviewed publications included any available in vitro and/or in vivo studies that showed a decreased or increased level of BDNF for each chemical included in the test set. Compounds were further classified as active/inactive for the CKEs relating to the decrease of synaptogenesis and neural network formation; this was based on the selectivity and specificity identified from the corresponding in vitro studies taking into account cell viability and the results of toxicity assays.

Other data collected included the calculated logarithm of the distribution coefficient (LogD, pH = 7.4), a measure of lipophilicity for ionizable compounds; blood–brain-barrier (BBB) permeability; and the capability to bind to the P-glycoprotein (P-gp) transporter, i.e., if the compound acts as a P-gp inhibitor or substrate, or is (non-)active against P-gp. Predictions for BBB permeability and P-gp interactions were made using in silico models based on curated SMILES [30], [31], [32]. BBB permeability is essential for understanding whether a compound crosses into, and has the possibility to act on, the central nervous system (CNS) [33]. P-gp is a transmembrane protein belonging to the ATP-binding cassette family of transporters (ABC-transporters), highly associated with the absorption, distribution, metabolism, excretion, and toxicity (ADMET) properties of compounds. P-gp may contribute to a decrease in toxicity by eliminating the compound from cells and preventing its intracellular accumulation [30].

The final data set contained missing information given by the in silico models for e.g., inorganic compounds, and there were compounds with no available evidence in the literature for the reduction of BDNF. In addition, inactive compounds in any of the in vitro assays were treated as missing information. Eight compounds contained complete details for this biological path. The final data set is provided in Table S3 in the supplementary material. The types of data utilised for modelling DNT are illustrated in Fig. 1, while further details are provided in Table 1.

Open in a separate window

Fig. 1

Types of information collected for model development exemplified for bisphenol A. This figure illustrates how different streams of data can be integrated for causal predictions to complement the information on key events. The full data set is provided in Table S3 in the supplementary material. The EC_x values were extracted from the corresponding in vitro studies as published.

Table 1

Types of data and their sources collected for the development of the Bayesian hierarchical model. The table describes all variables, i.e., predictors and outcomes defined as features included in the model for the type of data and performance where applicable. See also Tables S1-S2 in the supplementary material.

Feature	Description/Relevance	Data Type	Performance	Source
Chemical Name	The names used to define the compounds tested in both in vitro studies.	Not Applicable	Not Applicable	[26], [27]
CAS RN	Chemical Abstracts Service Registry Number associated with the tested compounds used to identify and track them during the modelling.	Not Applicable	Not Applicable	[26], [27]
DNT Classification	Each compound was classified as either positive, known or potential inducing DNT/negative, safe or without evidence for inducing DNT based on in vivo studies.	Binary, i.e., positive (i.e., associated with DNT) or negative	Not Applicable	[26], [27], [29]
LogD	The logarithm distribution coefficient calculated based on the compounds’ SMILES strings.	Continuous, unitless values	Not Applicable	ChemSpider database [34]
BBB	Each compound was classified for its capability to permeate the blood–brain-barrier (BBB) based on curated SMILES. Predicting BBB permeability means indicating whether compounds pass across the BBB. Compounds that cross the BBB have the potential to be CNS-active, whereas compounds that do not cross are expected to be CNS-inactive.	Binary, i.e., positive (BBB permeable) and negative	The in silico model available in admetSAR v2.0 has the area under the curve (AUC) with a range from 0.625 to 0.99.	Literature review Online BBB Predictor v.0.9 admetSAR v.2.0 [32], [35]
Cbrain/Cblood	In vivo blood–brain-barrier penetration represented as BB = [Brain]/[Blood], where [Brain] and [Blood] are the steady-state concentration of radiolabelled compounds in the brain and peripheral blood. High absorption to CNS had a value of more than 2.0, medium absorption: 2.0–0.1, and low absorption: less than 0.1. The predictions are based on in vivo data on rats.	Continuous, unitless values	The QSAR model of Ma et al. [31] had R = 0.955 with s = 0.232, used by PreADMET v.2.0 to model the predictions.	PreADMET v.2.0 [31], [36]
P-glycoprotein Status	Each compound was classified based on curated SMILES as a substrate or not, inhibitor or not, active or inactive for P-glycoprotein (P-gp) transporter using an in silico model.	Binary, i.e., yes or no for a compound acting as a substrate or an inhibitor, or activity against P-gp	The non-error rate and the average precision was 0.70 for the external validation set.	[30]

Open in a separate window

Feature	Description/Relevance	Data Type	Performance	Source
Reduction of BDNF	Each compound was classified as either positive, inducing the reduction of BDNF levels, or negative based on a literature search for in vivo studies, e.g., in rats and mice, and/or in vitro studies, e.g., human neuroblastoma SH-SY5Y cell line, embryonic mouse hypothalamus cell line.	Binary, i.e., positive (evidence showing alterations of BDNF), or negative	Not Applicable	Literature review of historical and peer-reviewed studies that evaluated compounds for their BDNF reduction potential
Decrease of Synaptogenesis	Selectivity and potency of a chemical were kept as classified in the reference based on results for viability and effective concentrations in rat primary cortical cells.	Categorical, i.e., inducing or not alterations of synaptogenesis	The battery assay had a sensitivity of 87% and a specificity of 71%.	[27]
Synaptogenesis Viability	The amount of ATP present in each well was calculated to assess compounds for their viability in rat primary cortical cells.	Continuous	Not Applicable	[27]
Synaptogenesis Activity, EC30 (?M)	30% change compared to control expressed as an effective concentration EC30 (?M) for puncta per total dendrite length (the most sensitive endpoint) measured in rat primary cortical cells for five days using an imaging assay.	Continuous	Not Applicable	[27]
Decrease of Neural Network Formation	Selectivity and potency of a chemical were kept the way the reference classified a compound based on results for viability and effective concentrations in rat primary cortical cells.	Categorical, i.e., inducing or not alterations of neural network formation	The model had a mean accuracy of 80.2%.	[26]
Neural Network Formation Viability	Total lactate dehydrogenase (LDH) release upon cell lysis and Alamar blue assay were used to assess compounds for their viability in rat primary cortical cells.	Continuous	Not Applicable	[26]
Neural Network Formation Activity, EC50min and EC50max (?M)	50% change compared to control expressed as an effective concentration EC50 (?M) with minimum and maximum values of all 17 parameters measured in rat primary cortical cells over 12 days using microelectrode array (MEA) recordings.	Continuous	Not Applicable	[26]

Open in a separate window

3. Results

The exploratory data analysis allowed for a better understanding of the data collected and utilised for the Bayesian analysis (Supplementary material, Figs. S2–S5). The statistical parameters of the model did not show the presence of any divergent transitions to approximate the posterior distributions (Supplementary materials, Figs. S6-S7 and Table S10). The estimates of the hyperpriors and priors are summarised in the supplementary material, Fig. S8.

The posterior distributions, (i.e., likelihood, posterior densities) and posterior predictive distributions were obtained from the Bayesian analysis. The posterior distributions represent the evidence provided by the data combined with the prior that incorporates our knowledge before analysing the data. The posterior probability distributions for each CKE, including the AO, are summarised in Fig. S9. The Bayesian CI was broader for compounds with missing information. Thus, the Bayesian CI quantifies the uncertainty given by the sources of variability and missingness. For instance, fluoxetine with information across all levels had a 95% HDI of 0.83–1.0 with a mean probability of 0.94 for the induction of DNT (Fig. 3). In contrast, sodium fluoride with missing data for the CKEs had a 95% HDI of 0.23–0.99 with a mean probability of 0.65 for the induction of DNT (Fig. 3). Glyphosate, a known negative DNT, had a very low 95% HDI with a maximum of 0.55 with a mean of 0.22 for the induction of DNT (Fig. 3).

Open in a separate window

Fig. 3

Visualisations of the predictions for three compounds: fluoxetine, sodium fluoride and glyphosate to show the relationships between the compounds and CKEs and the AO following the structure of the simplified biological path for DNT. The shaded blue distribution represents the predicted severity of the CKEs and the AO as well as the uncertainty in the prediction for the 95% highest density interval (HDI). For a complete overview of the results of all compounds, the reader is referred to Fig. S9 in the supplementary material. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

The posterior predictive distributions are given by the binary classification of compounds that were analysed against the likelihood. The binary classification was informed by the literature review for the CKE reduction of BDNF and by the in vitro studies for the other two CKEs and expressed as a Bernoulli distribution to describe the outcome in a Boolean way, as explained in the Materials and methods section. Two thresholds derived from the results of the predicted posterior distributions were used to classify the compounds for the low, medium or high probability of inducing a CKE and the AO. Such a classification might be used for screening and prioritisation purposes. All compounds with a probability between 0.59 and 0.60 for posterior predictive distributions were classified overall as having a medium probability of inducing a reduction of BDNF (Supplementary material, Fig. S10). This might be a result of the imputation method chosen to deal with the missing information. The results for the decrease of synaptogenesis showed a probability between 0.36 and 0.76 and as such medium and high levels of inducing this CKE (Supplementary material, Fig. S11). For the decrease of neural network formation, the posterior predictive probability was between 0.33 and 0.77 for all compounds that were classified as having low, medium and high levels of probability to induce this CKE (Supplementary material, Fig. S12). An overview of the distribution of the predicted posterior probabilities is shown in Fig. 4.A. The predictions for the AO, which is of most interest for e.g., regulatory decision making, helped to classify the compounds into the three levels (Fig. 4.B). Relatively few compounds had a low level of probability (5% of compounds), and most of the compounds had a high level of probability of inducing DNT (57% of compounds). This unbalanced classification is explained by the initial list of compounds for which 74% of compounds were known to be positive for DNT in vivo (in humans or animals). The full results for the prediction of CKEs and the AO are provided in Tables S6–S9 in the supplementary material.

Open in a separate window

Fig. 4

The results of the posterior predictive probabilities. A. The distribution of the posterior probabilities of the three CKEs and the AO for positive and negative compounds, and compounds with a level of missingness for the CKE reduction of BDNF. It shows how the mean of the predicted probability of compounds was clustered and distributed based on the two thresholds to describe a low, medium and high probability for the induction of the corresponding CKE and AO. A detailed graphical representation is shown in Supplementary material, Figs. S10–S12. B. Predicted probabilities of compounds for the induction of developmental neurotoxicity. The predicted probabilities are colour-coded based on two thresholds estimated from the results set to group the compounds for their low, medium and high probability. Compounds were listed in the order of increasing probability.

The sensitivity analysis underlined that weakly-informative priors did not have a significant influence on the model, which is instead data-driven (Table S11). The overall performance of the model is summarised in Table S12 with an accuracy of 76% for the prediction of DNT. There were several misclassifications due to a number of reasons including the quality of the data used for modelling and the level of the missing information associated with the CKEs. For example, a posterior predictive probability of 0.71 was obtained for loperamide hydrochloride, which was misclassified with a high level of probability for inducing DNT. The potential reasons for these misclassifications may be in part because of (1) the in silico model used for the prediction of P-gp, which classified the compound as an inhibitor, substrate and active, (2) the absence of data for the CKE of reduction of BDNF, (3) the in vitro assays that identified it as active, (4) the chosen threshold for classification. However, the large CI associated with the predicted mean value of 0.36–0.99 emphasises the caution required in making decisions for this compound.

Notes