Application of Bayesian networks in determining nanoparticle-induced cellular outcomes using transcriptomics

Abstract

Inroads have been made in our understanding of the risks posed to human health and the environment by nanoparticles (NPs) but this area requires continuous research and monitoring. Machine learning techniques have been applied to nanotoxicology with very encouraging results. This study deals with bridging physicochemical properties of NPs, experimental exposure conditions and in vitro characteristics with biological effects of NPs on a molecular cellular level from transcriptomics studies. The bridging is done by developing and implementing Bayesian Networks (BNs) with or without data preprocessing. The BN structures are derived either automatically or methodologically and compared. Early stage nanotoxicity measurements represent a challenge, not least when attempting to predict adverse outcomes and modeling is critical to understanding the biological effects of exposure to NPs.

The preprocessed data-driven BN showed improved performance over automatically structured BN and the BN with unprocessed datasets. The prestructured BN captures inter relationships between NP properties, exposure condition and in vitro characteristics and links those with cellular effects based on statistic correlation findings. Information gain analysis showed that exposure dose, NP and cell line variables were the most influential attributes in predicting the biological effects. The BN methodology proposed in this study successfully predicts a number of toxicologically relevant cellular disrupted biological processes such as cell cycle and proliferation pathways, cell adhesion and extracellular matrix responses, DNA damage and repair mechanisms etc., with a success rate >80%. The model validation from independent data shows a robust and promising methodology for incorporating transcriptomics outcomes in a hazard and, by extension, risk assessment modeling framework by predicting affected cellular functions from experimental conditions.

Keywords:

• Bayesian networks
• machine learning
• in vitro
• transcriptomics
• nanoparticles
• information gain

Introduction

The term nanoparticles (NPs) covers particulates with one or more dimensions in the sub-100 nm size range which can differ widely in composition, size, shape and solubility as well as other physicochemical properties. They are not a single entity and the diversity in properties is only matched by the range of possible biological interactions associated with such NP properties. This creates a significant challenge for pro-active risk characterization for the adequate prevention of human and/or environmental health effects whilst still reaping the societal benefits of nanotechnology. This is further compounded by the rate of development in nanotechnologies and their increasingly widespread use in commercial products such as electronic devices, food packaging and health care.

Human and environmental exposure to NPs is a current issue and can be both incidental (e.g. occupational exposure) or deliberate (e.g. cosmetics). NPs can enter the human body via a number of routes including inhalation, ingestion or skin absorption where they themselves can become modified and may also cause biological perturbations. The toxicity of different NPs has already been documented by a number of studies (Ajdary et al. 2018; Cao 2018; Gerber et al. 2013; Stensberg et al. 2011) utilizing both in vitro and in vivo models showing a broad range of property dependent effects from relatively minor (e.g. oxidative stress) to the more serious (e.g. mutagenicity). This variation in toxicity demonstrates that not all NPs are toxic and not all NPs are nontoxic, thereby stressing the need for efficient and predictive approaches to hazard characterization of different NPs.

Systems toxicology (Figure 1) extrapolates short-term observations to long-term outcomes, integrating information from experimental systems (molecular, cellular and macroscopic) with potential risks to humans (phenotypic and population level studies) (Belcastro et al. 2018; Boué et al. 2017). Omics technologies (Figure 1 left) study the changes in epigenome, transcriptome, proteome and metabolome induced by a substance. Integration of omics with the toxicity of a substance is called toxicogenomics and bridges molecular and cellular effects (NRC 2005; Sahu et al. 2015; George et al. 2010). Capturing dysfunctions at a molecular level allows omics technologies to trace adverse effects at low doses including early and sub-clinical effects that traditional in vitro and in vivo studies may overlook (Kawata, Osawa and Okabe 2009; Bai et al. 2017). Toxicokinetics models (Figure 1, center) such as Physiologically Based Pharmaco-Kinetic models (PBPK) and Physiologically Based Dose Response (PBDR) models enable in vitro-to-in vivo and in vivo-to-human extrapolation of observations (Judson et al. 2011; Li and Reineke 2012; Raies and Bajic 2016; Carlander et al. 2016). Findings from human exposure and epidemiological studies (Figure 1 right), supplemented by experimental studies, provide the information needed for substance safety regulation.

Figure 1. A systems toxicology paradigm. The pathway aims to extrapolate all available toxicological information to humans using a combination of in silico methods and toxicokinetics modeling. Molecular level analysis provides toxicological information preceding cellular and physiological manifestation of effects. In vitro methods investigate cellular effects of direct exposure to a substance. In vivo studies provide an examination of a variety of tissues in controlled environments. Human studies require biomonitoring, sampling and tissues/fluids analysis. Population level effects is the culmination of this extrapolation with supporting evidence typically requiring longitudinal observations and large-scale data gathering.

Adverse outcome pathways (AOPs) systematically address the organization and gathering of systems of toxicology information for specific substances. AOPs have been developed as identifiable event chains that link the effects at molecular level to human level (OECD 2017). Although recent trends in toxicology focus on toxicity pathways, mode of actions and AOPs (Raies and Bajic 2016; NRC 2007; Shukla et al. 2010; Nel et al. 2013), definite AOPs have not yet been established for specific NPs or classes of NPs (e.g. poorly soluble, low toxicity NPs) due to toxicity data gaps (Gerloff et al. 2017), although there are areas of advancement (Vietti, Lison and van den Brule 2016). Increasingly, omics data have become more widely available given the increase in high-throughput and high content screening methods. However, the development of state-of-the-art knowledge that integrates such data from these various methods with more traditional animal and epidemiological information, is still a developing process (Stagljar 2016; Labib et al. 2016; Buesen et al. 2017; Sauer et al. 2017).

Although toxicity has been observed with a variety of NPs to variable extents and outcomes, how the specific physicochemical properties of different NPs and cellular properties of exposed systems relate to adverse effects is not comprehensively understood. The availability of data surrounding the interaction between biological systems and NPs is still relatively limited as compared to other chemical compounds (Bai et al. 2017). Dealing with the variety of NPs physicochemical properties, routes of exposures, dosimetry, and their multiple effects on biological systems means that Machine Learning (ML) tools are particularly well suited toward the prediction of biological effects based on NPs properties (Winkler et al. 2013; Marvin et al. 2017; Sizochenko et al. 2014). Furthermore, computational models are fast and cheap and as they rely on information inputs rather than physical test materials, they can be used to predict the impact of materials not yet synthesized, thereby contributing to the development process, particularly in relation to safe-by-design approaches. This efficient approach assists in the development of safer Nanomaterials (NMs) without stifling their commercial development (Winkler et al. 2014) in the way that iterative processes of physical material development and in vitro/in vivo safety testing can through constraints of cost and time. Most modeling efforts have thus far been based on Quantitative Structure-Activity Relationships models (QSAR or nano-QSAR) using nano-descriptors and physicochemical properties (Napierska et al. 2010; Choi et al. 2018; Chen et al. 2017). QSAR models do not provide insights into the toxicity pathways, but instead predict endpoints from in vivo or in vitro studies (e.g. carcinogenic or noncarcinogenic) (Jagiello et al. 2016). Recently, novel models called RASARs (Read-Across Structure Activity Relationships) are being used to define chemical similarity (Luechtefeld et al. 2018). Another emerging recent approach, Quantitative Structure-Toxicity Relationships (QSTR) perturbation models, have been applied to estimate toxicity and ecotoxicity of NPs for different endpoints and experimental conditions with remarkable performance (Kleandrova et al. 2014a; Luan et al. 2014; Kleandrovaet al. 2014b, González-Durruthy et al. 2017). QSTR-perturbation models apply moving average analysis to reconstruct the case descriptors and perturbation theory to get the final optimization function based on differences of case pair combinations. Whether based on data fusion or approximate methods such approaches are promising in the field of nanotoxicology.

Bayesian networks (BNs) are ML probabilistic graphical models that merge automated probabilistic analysis and human insight for a wide range of problem solving. BNs have been used in a plethora of studies for forecasting complex system outputs, including risk evaluation, medical diagnosis, skin sensitization and clinical decision-making (Morgan 2005; Jaworska, Gabbert and Aldenberg, 2010; Jaworska et al. 2013). They are increasingly relied upon on the modern, data rich society we find ourselves in and unsurprisingly are finding a significant role in nanosafety and have been used to support risk assessment and prioritize NMs hazard assessment (Marvin et al. 2017; Murphy et al. 2016; Sheehan et al. 2018).

The aim of this study is to develop and apply a BN classifier that predicts NP-induced cellular effects using data inputs addressing NP physicochemical properties, experimental exposure conditions in vitro characteristics and transcriptomics outcomes. The study data gathered for output use functional analysis of the altered genes; investigating which biological pathways are most disrupted after exposure to a test substance (Römer et al. 2014; McDermottet al. 2013; Lin and Lane 2017). This represents a growing field within toxicology because of its capacity to improve risk assessment (Shukla et al. 2010; Buesen et al. 2017; Tsiliki, Nymark, Kohonen, Grafström and Sarimveis 2017; Franceschi et al. 2017). Using this approach, we seek to predict disruption of biological processes owing to genetic perturbations generated by physicochemical properties of NPs. The focus of this study is to determine the predictive capability of a modeling approach focused on biological pathways in vitro which align with AOPs relevant to nanotoxicity.

Materials and methods

Approach

BNs were constructed and evaluated as outlined in Figure 2 with the input data derived through a review of the peer reviewed, scientific literature. The extracted data were then applied either in its unprocessed, raw form (Raw Dataset, RD) or processed to address imbalance, data gaps, and discretization of values to form the Processed Dataset (PD). Imbalance issues in the dataset occurs when one of the outputs has more counts than the other which in turn can limit classifiers performance (Choi et al. 2018; Wang et al. 2018).

Figure 2. Workflow of BN construction and evaluation.

Data were then analyzed using a Correlation Feature Selector (CFS) algorithm, which provides some initial insight into the relationships among the attributes (inputs and outputs). Based on these relationships, the BNs were manually crafted into prestructured networks with the two datasets, RD and PD, randomly split into training (60%) and validation (40%) sets. The training PD set was enriched by applying a Synthetic Minority Technique (SMOTE) to address imbalance issues. For each dataset, different BNs were constructed; a BN based on the RD (BN-RD) and another based on the PD (BN-PD). A third BN with a structure learned from the data using the K2 algorithm was also applied for the PD (BN-K2). The K2 algorithm uses a greedy search with no restriction on the number of parents a node can have (Lerner and Malka 2011) while BN-RD and BN-PD have maximum four parents in each node. All BNs are trained and applied as binary classifiers individually for each of the nine outcomes (Furxhi et al. 2018). At the final stage, all BNs were validated against the training and validation sets with an additional validation, using a reliability set, finally used to evaluate the BN performance when data from new studies are introduced.

Data collection and curation

Toxicity outcomes of the studies gathered

Data from in vitro differential gene expression analysis studies were gathered from peer-reviewed scientific literature. Omics methods in nanotoxicology is advocated due their ability to capture different outcomes. Whole-genome transcriptomics studies were selected since the use of existing standardized tests might miss, in principle, novel endpoints Marchese Robinson et al. (2016). In targeted analysis, such as targeted genes sequencing using polymerase chain reaction (PCR) only the regions of select set of genes are analyzed. Unexpected changes outside the examining genes are not observed such as in whole-genome analysis (Boon and Faas 2013; Hrdlickova, Toloue and Tian, 2017). We collected non-targeted analyses carried out for the complete genome (whole-genome analysis) in human cell lines using a systematic search strategy using Boolean logic operators (AND, OR and NOT) combined with key search terms (e.g. nanoparticle, nanomaterial, transcriptomics, gene expression, whole-genome) to form defined search strings. These search strings were applied to various publicly available search engines and databases (Google Scholar, ScienceDirect, PubMed, NCBI GEO Datasets, NanoMiner and ArrayExpress). A summary of the retrieved literature is available at (Furxhi 2019).

Data completeness is a measure of the extent to which the data which serve to address a specific need are available (Marchese Robinson et al. 2016). In our case, the focus was studies that extract enriched Gene Ontology analysis terms from the most significant dysregulated expressed genes that are associated with biological processes. Almost all studies selected had, as a first minimum requirement, to have validated quality and integrity of RNA extraction prior to microarray analysis. In addition, most of the studies further validate the results by applying real-time quantitative PCR analysis. Furthermore, in order to maintain high relevance for human health effects, all studies gathered here analyze human cell lines after exposure to NPs (Marchese Robinson et al. 2016).

The studies gathered had a considerably variety of different human cell models including both primary cells and immortalized cell lines which exhibit varying responses. In addition, the cellular models were restricted to human cells to limit uncertainties related to species-specific biological mechanisms (Judson et al. 2011; Shukla et al. 2010; Jaworska, Gabbert and Aldenberg 2010). The use of a broad variety of cell models was to breadth of biological interactions arising from a variety of possible NP exposure routes. For example, the use of type-II alveolar epithelial (A549) cells and macrophage cells (e.g. differentiated monocyte derived macrophages) to represent inhalation exposure, intestinal cells (e.g. Caco-2) to represent ingestion route and gastrointestinal exposure, endothelial cells to reflect cardiovascular exposure, either deliberate (intravenous delivery) or incidental (translocation).

Toxicity determinants in the studies gathered

Numerous factors can determine the relative toxicity of substances compounds, such as the route of exposure (e.g. inhalation, oral, dermal), external exposure, tissue dose, persistence in the biological environment (duration of biological interaction), as well as physicochemical properties of NPs (shape, core size, coating etc.). Furthermore, different tissue types, cell models and experimental parameters (e.g. particle dispersion and cell culture conditions) add further data variables.

A challenge associated with defining information criteria for NPs is that the current understanding of independent variables, such as physicochemical properties and other experimental variables, which contribute most to the outputs of assays is arguably insufficient (Marchese et al. 2016). The lack of harmonization in toxicological analysis, NP characterization and preparation as well analysis and reporting within the nanotoxicological literature, driven in part lack of standard protocols for NP in vitro toxicity studies (Drasler et al. 2017; Warheit 2018) generates significant heterogeneity of the literature data.

In this study, we collected all available information on NP characterization from each study keeping a minimum requirements check list in order to identify in silico the physicochemical properties and other experimental variables that significantly affect the outputs of transcriptomics studies. The studies selected include “priority” properties of interest for NP characterization during toxicological assessment as reported in (Hackley et al. 2013). Besides RNA extraction quality check, the second minimum requirement for a study to be considered was reporting NP size and a traceable origin, i.e., either the synthesis method of the NP, or product name and manufacturer for commercially available NPs. Additional physicochemical parameters (such as surface area) are deemed important for the data to be considered complete (Hackley et al. 2013; Marchese Robinson et al. 2016). Studies missing information for such attributes were also included when this information could be retrieved from the NP manufacturing or synthesis information.

Within nanoscience, less attention has been paid to the question of which additional experimental details used, need to be recorded ((Marchese Robinson et al. 2016). In this study, we gathered experimental and in vitro characteristics, such as exposure dose and duration and cell line, cell type and microarray protocol. A third minimum requirement for a study to be considered was to report exposure dose and duration, with one exception, (Peng et al. 2010) where it was not relevant. Cell line would have been a fourth requirement but it was, in any case, always reported.

In general, information on physicochemical properties was scarce but traceable, whereas, exposure (dose and duration) and in vitro conditions were always reported. Supplementary online material provides a classification of information given in the selected studies (Furxhi 2019). All studies provide information for the majority of the toxicity determinants/methods/outputs classes.

For each of the studies reviewed, the combination of experimental exposure conditions, in vitro characteristics, such as cell line, cell type and tissue, physicochemical properties and outcomes are recorded separately in a datasheet, resulting in many extractions per study. If, for instance, another dose is administered during an experiment with all the other conditions unchanged, this was inserted as new instance within the dataset.

Biological effects comprising a number of cellular responses (see Table 1 outputs) were extracted and used as the outcome of interest to be simulated by the BN model. Effect occurrence was expressed in binary manner (Triggered, No effect). The NP is said to have a toxic effect (Triggered) based on the conclusions of the relevant study from which the data was derived. The studies analyzed dysregulated genes filtered with a cutoff criteria (p-value or log2 fold change) in order to identify the most significant perturbed genes at in vitro systems after exposure to NPs. Gene Ontology analysis was then performed to functionally annotate the candidate genes with the related biological effect. Triggered output is defined as the biological effect that correspond to cellular responses of disturbed gene functions after NP exposure. While ‘No effect’ was determined were cellular processes were not significantly disturbed after exposure to a NP.

Table 1. Processed Dataset (PD) inputs and outputs extracted from literature review and their states (categories)/numerical ranges and abbreviations used in this study.

Nodes	Abbreviations	States (categories)
Inputs
NPs Physicoch-emical Properties
Nanoparticle	NP	10 categories: Ag, Au, Polymeric NPs, CuO, ZnO, TiO2, SiO2, Fe2O3, Polystyrene NPs, CoFe2O4 NPs
Core size (nm)	Core size	4 discretized ranges: 0 to 10, 10 to 50, 50 to 100, >100
Shape	Shape	5 categories: Spherical, Ellipsoid, Unknown, Nanotube, Clusters
Surface coatings	Coatings	21 Categories: None, Silica, Citrate, PAH, DMSA, Madelic acid, Lipid bilayer, Hybrid lipid, Methoxyl, Diethyl glyc, Gluconic acid, Citric acid, folic acid, chromium, starch, surfactant, polymeric, colloidal, polyethylenimine, mercaptopropyl trimethoxysilane, carboxyl
Specific Surface Area (m2/g)	Surf. area	9 discretized ranges: 0 to 18, 18 to 36, 36 to 54, 54 to 72, 72 to 90, 90 to 126, 126 to 198, 198 to 216, >216
Zeta Potential (mV)	Zeta	5 discretized ranges: -30 to -15, -15 to 0, 0 to 15, 15 to 30, >30
Exposure Conditions
Exposure dose (μg/ml)	Dose	8 discretized ranges: 0 to 3, 3 to 8, 8 to 12, 12 to 30, 30 to 50, 50 to 75, 75 to 100, >100
Exposure duration (h)	Time	5 discretized ranges: 0 to 6, 6 to 12, 12 to 24, 24 to 48, >48
In vitro model characteris-tics
Cell line	C. line	22 categories: SH-SY5Y, 293T, A549, CACO-2, HDF, PC3, THP-1, HEPG2, VSMC, HACAT, HMDM, JURKAT-T, MDDC, MCF-7, IMR-90, U251, HELA, HMEC 184, EAHY926, SAE, RKO, SK MEL-28
Cell type	C. type	2 categories: Cancer or Normal
Tissue	Tissue	12 categories: Kidney, Brain, Lung, Intestinal, Skin, Prostate, Liver, Colon, Cardiovascular, Breast, Cervix, Blood
Microarray	Array	4 categories: Genechip_Affymetrix, WHG_Agilent, Beadchip_Illumina, Superarray_Bioscience
Outputs
Cellular Effects
Cell cycle and Proliferation responses	C. cycle	Binary mode (Triggered, No effect)
Cell death and Apoptosis responses	C. death
DNA damage and Repair responses	DNA dam.
Cell adhesion and Extracellular matrix responses	C. adhe
Inflammation and Immune responses	Inflam.
Unfolded protein responses (UPR) and Endoplasmatic  reticulum (ER) stress	Stress
Metal Ion Responses	Metal ion
Angiogenenesis responses	Angiogen.
Cytoskeleton organization responses	Cytoskel.

Data preprocessing

As such, missing values in the compiled dataset were not uncommon, with several physicochemical properties not reported in most of the studies, most notably composition (presence of impurities in bulk and coatings) and surface area.

Filling the missing values from articles ‘Supplementary information’, as well as information gathered directly from the manufacturers’ (e.g. material safety data sheets) which resulted in a second dataset (PD). Density values derived from manufactures information combined with particle size distribution data (assuming spherical shape and smooth surface) can be used to calculate the Specific Surface Area (SSA) (Ha et al. 2018). Furthermore, as well as data gap filling, discretization was also performed during preprocessing. The range of the numeric input attributes in the PD is discretized and expressed as nominal equal frequented bins with some input ranges obtained from literature, such as core size and zeta potential (Hristozov et al. 2014).

Data analysis

Correlation Feature Selection algorithm (CFS) was applied to the data to remove irrelevant and redundant information allowing a classifier to operate more effectively (Kohavi and John 1997; Doshi and Chaturvedi 2014). WEKA (Waikato Environment for Knowledge Analysis, version 3.8.2), an open-source Machine Learning workbench was used for the CFS application and BNs construction. We used the algorithm to investigate associations between the variables individually for the two datasets (Hall and Smith 1998). CFS was applied in 10-fold cross-validation (out-of-sample testing) which randomly selects ten equal size subsamples of the dataset for validation. Of the ten subsamples, a single one is used as validation set to test the remaining nine, combining the results to an overall average metric. We used CFS to estimate correlation between subset of attributes and outcomes, as well as inter-correlations between the attributes (Karegowda, Manjunath and Jayaram 2010). Distinguishing the most important associations provided the basis of constructing the BN graphical models.

Data split and balancing

The datasets are split into training (60%) and validation sets (40%) randomly. Most instances in the training dataset correspond to no effect outcomes for the majority of the endpoints. Imbalanced datasets can limit the performance of most classification algorithms, making the prediction biased to the dominant class value (Ha et al. 2018). To avoid this, we adjusted the relative frequency of triggered/no effect instances by resampling the second dataset by applying SMOTE (Synthetic Minority Oversampling Technique), a supervised instance algorithm that oversamples the minority instances using the k-nearest neighbors algorithm (Chawla et al. 2002).

Bayesian networks

A BN is a directed acyclic graphical model that represents variables as ‘nodes’ and their connections as arrows. Each arrow signifies a conditional dependence of the child node to a parent node. The network as a whole represents the joint probability distribution of included variables and use Bayes’ Rule to update conditional probabilities given new data (Zabinski et al. 2016; Zabinski, Pieper and Gibson, 2017; Friedman et al. 2000). The network is practically defined by its graphical structure and the corresponding trained conditional probabilities of the included variables. The probability distribution of all nodes $\mathrm{P}\left(\mathrm{A}\right)=\mathrm{ }\mathrm{P}\left({\mathrm{A}}_1,\mathrm{ }\dots ,{\mathrm{A}}_{\mathrm{n}}\right)$ is specified by the product of all conditional probability tables (CPTs) in the BN model as shown in Equation (1)(1) $$\mathrm{ }P\left(A\right)=\mathrm{ }\prod _{i=1}^{n}P[\left(A_i|\mathit{pa}(A_i\right)]\mathrm{ }$$(1) , where $\mathit{pa}\left(A_i\right)$ is parent of node $A_i$ and $\mathrm{P}\left(A_i|\mathit{pa}(A_i\right)$ a conditional probability distribution. (1) $$\mathrm{ }P\left(A\right)=\mathrm{ }\prod _{i=1}^{n}P[\left(A_i|\mathit{pa}(A_i\right)]\mathrm{ }$$(1)

CFS results between all variables in the dataset (input and output) included inter-relationships between the outputs. A model that would have outcome-to-outcome links would fail to predict the outcomes only from new inputs and therefore would be biased to the prediction of other effects. Applying CFS for each biological effect of the PD dataset separately reveals the association (predictive ability) of each input attribute (all cases of the final PD dataset) for the specific effect (Figure 3). Based on that, we constructed the links that associate input features as parents to the cellular effects enabling the BN model to predict the effects solely from the inputs.

Figure 3. Results of 10-fold cross-validation Correlation Feature Selection algorithm applied in processed dataset. (A) Inter-relationship of input variables, (B) relationship between the input and output variables.

Another BN structure is derived automatically by WEKA using the K2 algorithm (Gámez, Mateo and Puerta 2011). Bayesian rules, such as the K2 algorithm, build a different BN structure for each output by changing both probabilities and links between the inputs and the outputs. K2 algorithm is a score-based algorithm in BN which finds the structure that maximizes the score (Lerner and Malka 2011).

Model validation

Training sets were used for the internal validation of the BNs. The internal validation was performed using a 10-fold cross-validation process which valuated the predictive classifiers by partitioning the original dataset into a ‘training set’ and a ‘test set’. Cross-validation is a standard way to obtain unbiased estimates of a model's goodness of fit. In our case we used 10 subsamples as training data where a single subsample is retained as the validation data. The external validation set (40% of the whole dataset) was used to evaluate the predictive performance ability of the trained model (Raies and Bajic 2016). In addition to the external validation dataset, a reliability validation set (four studies) was used with data from studies that were not included in the previous training and validation datasets. Matthews Correlation Coefficient metric is also calculated (Matthews 1975).

Important attribute analysis

Variable (attribute) importance can be measured relatively, based on information obtained from the model (Choi et al. 2018). We quantified the importance of a feature for classification using information entropy (Doshi and Chaturvedi 2014; Nguyen, Petrović and Franke 2010; Fan et al. 2011) which is a measure of uncertainty of a random variable and defined as: (2) $$H\left(X\right)=\mathrm{ }-\sum P\left(X_i\right)l\mathit{og}2\left(P\right(X_i))$$(2) where, $H\left(X\right)$ is the entropy of variable X and $P\left(X_i\right)$ is the prior probabilities for the i^th value of X. The entropy of X after observing values of another variable Y is defined as (3) $$H\left(X/Y\right)=\mathrm{ }-\sum P\left(y_j\right)\sum P\left(X_i/y_j\right)l\mathit{og}2P(X_i/y_j))$$(3) where, $P\left(X_i/y_j\right)$ is the posterior probabilities of X when values of Y are given.

The amount by which the entropy of X decreases reflects an addition of information about X provided by Y. The reduction of entropy is called information gain and can be defined as: (4) $$\mathit{IG}\left(X/Y\right)=H\left(X\right)-H(X/Y)$$(4) where, $\mathit{IG}\left(X/Y\right)$ is the information gain of variable X due to entropy decrease provided by variable Y.

In order to evaluate the relevant importance of every input to each output we normalized IG for each effect using min-max normalization (Priddy and Keller 2005). The input with the lowest IG, that is, the one with the least influence on the model outcome prediction is zeroed and the most influential attribute scores 1.

Results

Data pre-processing

Information in raw dataset (RD) on the specific surface area (SSA) was not available in 98% of the 245 cases (rows at the datasheet); the presence of impurities in the bulk form or in the coating was not defined in approximately 88%. In 29% of cases, no information about the zeta potential was given, and in 75% of cases, there were no shape descriptors included. Of the various attributes relevant to particle toxicity, primary particle size was the only value consistently reported. More generally, information on physicochemical properties was scarce whereas, exposure (dose and duration) and in vitro conditions, such as cell line, cell type and microarray protocol, were always reported. In relation to the NPs evaluated in the studies extracted, the most common were: Ag (16.3%), Au (18.2%), TiO₂ (17.1%) and ZnO (47.2%). Other NPs appeared (CuO, SiO2, Fe₂O₃, etc.) in less than 2%.

The PD consisted of 245 instances (rows) of 12 inputs and 9 outputs as columns (Table 1). Processed Dataset (PD) had 8% (RD: 75%) of shape missing, 9% (RD: 29%) of zeta potential missing and 5% (RD: 98%) of SSA missing. Reliability dataset was comprised of 19 instances. The PD and reliability datasets are available at (Furxhi 2019).

Data analysis

The relationships between the inputs in the PD dataset are shown in Figure 3(A) and the relationships of each outcome individually to the inputs are shown in Figure 3(B). In both cases, attributes showing the most robust correlation (association 0.9–1) are selected as parents to the respective variable (node).

These input relationships allow us to manually construct the part of the model that is fixed for all outcomes. We build the prediction model by relating each output to the fixed input structure, excluding inter-relationship of outputs that would interfere with the results, which enables the ability to obtain various toxicological outcomes using one prediction model.

Bayesian structure

Two BN graphs were constructed on the basis of the results of the Correlation Feature Selection (CFS) algorithm which reveals significantly associated variables that as nodes, can be linked in a parent–child relationship. We constructed BNs nodes and links that represent the associations only between the inputs based on the results of CFS (Figure 4 center). Each output is linked with the input nodes (Figure 4 peripheral).

Figure 4. Graphical structure of the BN-PD model. Center: Structure based on correlation between the inputs. Perimeter: BN-PD structures based on output prediction ability of inputs.

Validation of the BNs

The classification performance of the three models, BN-PD, BN-RD and BN-K2, are shown in Figure 5–7 as the classifier accuracy, sensitivity and Matthewss correlation coefficient (MCC). Accuracy shows how many instances are correctly classified by the model compared to the total number of instances. Sensitivity demonstrates the classifiers’ ability to predict the triggered instances. MCC is a correlation coefficient between target and predictions. It varies between −1 and +1, lowest values indicating disagreement between actuals and prediction. It is frequently used to measure quality of binary classifications.

Figure 5. Accuracy of the classifiers. The validation of the models is carried out for each cellular effect separately. (A) Internal validation using a 10-fold cross-validation process with the training 60% datasets. (B) External validation using the validation 40% datasets. (C) Reliability validation using the reliability dataset. BN-PD: BN constructed from pre-processed dataset; BN-K2: BN automated constructed using K2-algorithm; BN-RD: BN constructed from raw dataset.

Figure 6. Sensitivity of the classifiers. The validation of the models is carried out for each cellular effect separately. (A) Internal validation using a 10-fold cross-validation process with the training 60% datasets. (B) External validation using the validation 40% datasets. (C) Reliability validation using the reliability dataset. BN-PD: BN constructed from pre-processed dataset; BN-K2: BN automated constructed using K2-algorithm; BN-RD: BN constructed from raw dataset.

Figure 7. Matthews correlation coefficient of the classifiers. The validation of the models is carried out for each cellular effect separately. (A) Internal validation using a 10-fold cross-validation process with the training 60% datasets. (B) External validation using the validation 40% datasets. (C) Reliability validation using the reliability dataset. BN-PD: BN constructed from pre-processed dataset; BN-K2: BN automated constructed using K2-algorithm; BN-RD: BN constructed from raw dataset.

MCC is a correlation coefficient between target and predictions. It generally varies between −1 and +1. −1 when there is perfect disagreement between actuals and prediction, 1 when there is a perfect agreement. It can be more appropriate when negatives actually mean something and it is frequently used to measure quality of binary classifications.

BN-PD showed the highest accuracy for all but one of the cellular effects in internal, external and reliability validation (Figure 5). Model performance improved when data preprocessing was implemented in the second dataset (BN-PD versus BN-RD). BN-PD also had higher performance compared to the automatically structured BN (BN-K2) in all of the cases. BN-RD had a higher accuracy in the case of inflammation (Inflam.) and metal ion responses in internal validation compared to the other two classifiers (Figure 5(A)). In addition, it had higher accuracy in the external validation in the case of metal ion responses (Figure 5(B)). Regarding the reliability dataset, BN-PD showed higher predictive capability than the other classifiers (Figure 5(C)). In general, data-driven approach structure performed better compare to BN-RD and BN-K2.

BN-PD showed the highest sensitivity for most of the cellular effects in internal, external and reliability validation (Figure 6) with model performance improved when data pre-processing was implemented in the second dataset (BN-PD). BN-PD also had higher performance compared to the automatically structured BN (BN-K2) in almost all of the cases. BN-RD had better sensitivity in cell cycle and inflammatory responses at internal validation, whilst BN-K2 had better sensitivity in DNA damage (DNA dam.), unfolded protein and endoplasmatic reticulum (Stress) responses (Figure 6(A)). BN-RD sensitivity scored higher compared to BN-K2 in predicting cell death (C. death), metal ion and cytoskeleton organization (Cytoskel.) responses in the external validation process (Figure 6(B)). Regarding the reliability dataset, BN-PD showed higher predictive capability than the other classifiers (Figure 6(C)). However, as shown in the sensitivity reliability graph (C), BN-RD predicted the triggered instances for DNA damage (DNA dam.) responses, while the other BNs failed to predict any.

BN-PD data-driven approach structure also had higher MCC performance compared to the BN-K2 and BN-RD (Figure 7). BN-PD outperformed the other classifiers in MCC metric in internal, external and reliability validation. BN-RD model had higher MCC score in the case of inflammation (Inflam.) and metal ion responses in internal validation (Figure 7(A)) compared to the other two classifiers and higher MCC in metal ion responses in external validation (Figure 7(B)). Regarding the reliability dataset, BN-PD showed higher predictive capability than the other classifiers (Figure 7(C)). Even for DNA damage (DNA dam.) responses BN-PD is zero compared to negative (worse) performance of BN-RD and BN-K2.

In summary, BN-PD that fully accounts for the CFS correlations had an aggregate predictive performance of >80% in most of the cases and outperformed the BNs that do not used either predefined structure, gap filling or data balancing.

Important attribute analysis

In conducting the important attribute analysis, it was found that exposure dose was the most influential attribute (IG = 1) that affected all outcomes significantly. Unfolded protein and endoplasmatic reticulum stress responses (Stress) were most influenced by NP which, given the classification of some NPs, relates to composition. Attributes such as NP and cell line ranked after dose as important attributes (0.25–0.82) while coating (0.59) strongly effected cytoskeleton organization (Cytoskel.) responses although coating appeared influential in most of the outcomes. Array and tissue occurred in top five determinants in four out of nine outcomes. Exposure duration (time) was influential on cytoskeleton organization, angiogenesis (Angiogen.) and DNA damage and repair (DNA dam.) responses. (0.14–0.31). Zeta potential appeared only in metal ion and unfolded protein (Stress) responses scoring 0.12 and 0.44, respectively. Surface area appeared solely in angiogenesis responses (0.27). Similarly, cell type occurred only in DNA damage and repair responses (0.25). Shape had the smallest information gain to the outcomes compared to the other four most influential input attributes. The input data represented low aspect ratio NPs therefore cannot address fibre-type NP effects nor can be used yet to perform a proper evaluation of how shape differentiates toxicological outputs. Primary particle size showed relatively insignificant influence to the prediction of the outcomes.

In summary, dose, NP composition and cell line had the foremost influences. Array, coating, exposure duration and tissue had notable influences. Zeta, cell type and surface area influenced some outcomes distinctively. Shape and core size in the given set of attributes had the least influence in discriminating the outcome classes. The results of normalized information gain analysis are presented in Table 2.

Table 2. Normalized information gain for each model outcome (bold) by input attributes. Only the top five most influential attributes are shown.

Angiogenesis responses		Metal Ion responses		Inflammation and Immune responses
Input attribute	Information Gain	Input attribute	Information Gain	Input attribute	Information Gain
Exposure dose	1.00	Exposure dose	1.00	Exposure dose	1.00
Nanoparticle	0.61	Nanoparticle	0.82	Cell line	0.25
Microarray	0.37	Zeta potential	0.12	Coating	0.01
Surface area	0.27	Cell line	0.03	Shape	0.00
Exposure duration	0.21	Tissue	0.02	Tissue	0.00
Unfolded protein responses and Endoplasmatic reticulum		Cell death and Apoptosis responses		Cell adhesion and Extracellular matrix responses
Input attribute	Information Gain	Input attribute	Information Gain	Input attribute	Information Gain
Nanoparticle	1.00	Exposure dose	1.00	Exposure dose	1.00
Cell line	0.73	Cell line	0.24	Nanoparticle	0.45
Zeta potential	0.44	Tissue	0.20	Microarray	0.14
Exposure dose	0.42	Coating	0.05	Tissue	0.13
Coating	0.08	Nanoparticle	0.03	Coating	0.01
Cell cycle and Proliferation responses		Cytoskeleton organization responses		DNA damage and Repair responses
Input attribute	Information Gain	Input attribute	Information Gain	Input attribute	Information Gain
Exposure dose	1.00	Exposure dose	1.00	Exposure dose	1.00
Cell line	0.32	Coating	0.59	Cell line	0.59
Nanoparticle	0.28	Nanoparticle	0.42	Nanoparticle	0.43
Microarray	0.20	Microarray	0.23	Exposure duration	0.31
Coating	0.15	Exposure duration	0.14	Cell type	0.25

Discussion

In this study, we describe the application of a BN approach to assess the relationship between in vitro wide-genome transcriptomic outcomes in various cellular models and exposure to different NPs. This requires the integration of biological outcomes (possible markers of toxicity) with physicochemical NP properties, exposure conditions and in vitro model characteristics. The BN model predicts a number of biological effects at cellular level with a predictive capability exceeding 80% based on data derived solely from the literature. Given the vagaries of reporting within the literature, this predictive capability is noteworthy and may be improved further with greater harmonization.

Bayesian network constructs can combine expert knowledge and data-driven modeling (Seixas et al. 2014; Velikova et al. 2013). Using an intermediate feature selection step before constructing and applying the BN model provides insight into the dataset under investigation and guides model development. Feature selection as a classifier prestep analysis has been implemented in several studies (Hall and Smith 1998; Drugan and Wiering 2010; Chormunge and Jena 2018; Cinicioglu and Büyükuğur 2014). In this study, we derive and compare two pre-structure BN based on CFS results (BN-PD and BN-RD) and show that using pre-processed data enhances the classifier accuracy (BN-PD). We compare the CFS-based BN-PD to an automatically derived BN structure using the K2-algorithm (BN-K2); the former performs better suggesting, manual BN structuring based on statistical correlation is optimum. In conclusion, predefined discretized ranges, gap filling and class balancing improves classifier prediction ability.

Previous studies have listed characteristics for aggregated hazard ranking of NMs to investigate the relationship between NMs properties and adverse effects (Marvin et al. 2017; Choi et al. 2018; Sheehan et al. 2018). In this study, artificial nodes of cumulative indices, such as hazard node derived from all the outputs were not included in the BN structure. This allows analysis and validation of the model for each cellular effect individually. In agreement with our findings, Marvin et al. (2017) found that coating plays a significant role in most biological effects such as inflammation and genotoxicity (DNA dam.). From a biological perspective, this would be expected as coating may serve to modify surface properties (e.g. charge, hydrophobicity) and biological interactions with poorly soluble particles occurs at the surface interface. Sheehan et al. (2018) performed a comparative study investigating the efficacy of quantitative Weight of Evidence and Bayesian methodologies in ranking the potential hazard and found particle size, surface coating, administration route and surface charge to have significant influence on the toxicological outcomes depending on the NP. These are all contributors also found to have significant impact on the relative toxicity of particles in a pulmonary inflammation review of in vivo studies (Braakhuis et al. 2014). In our study, zeta potential (surface charge) and tissue (administration route) have a significant influence, but size does not significantly affect the prediction of the molecular biological effects compared to the other attributes. The former is expected as zeta-potential is associated with respiratory toxicity (Cho et al. 2012) and different tissues show different senstivities towards NPs. In a study of 15 different metal oxide NPs, Cho et al. (2012) found that zeta potential was a key physicochemical influencer of respiratory toxicity in vivo. How size variation alone affects the prediction of outcomes can only be seen if the other effects are filtered out by means of a sensitivity analysis (out of the scope of this study). Studies report the effect of particle size on whole-genome gene expression (Bouwmeester et al. 2011; van der Zande et al. 2016; Grzincic et al. 2015). Particle size in terms of agglomeration state and/or primary particle size may change during a toxicity test (biological matrices) or during the life cycle due to dissolution or interaction with test media (Worth et al. 2017). In principle, the approach presented here can correlate different set of properties, quantitative exposure metrics, including in vitro characteristics with the toxicity pathways that NPs may trigger and accounting also for secondary effects.

Including exposure conditions, such as dose and duration, can improve the classification quality as this information is always reported in toxicological studies. Aggregating exposure features not readily comparable in one attribute, for example, exposure doses for different tissues, could reduce the biological accuracy of the model. This is compensated in our study by including tissue as an extra model input attribute. Exposure dose was identified as the most important determinant for most effects in the model, which is in concordance with the cornerstone of toxicology, the dose–response relationship. Such findings were also commensurate with QSAR as Choi et al. (2018) found that dose was the most important attribute when looking at cell viability using neural network QSARs. They reported formation entropy as the second most important attribute which strongly relates to NP in our model. Exposure dose here was clearly an important determinant for toxicity as it is used as a surrogate of delivered (tissue) dose. Tissue dose data from in vitro and in vivo studies could be equivalently combined in the dataset. However, tissue doses are widely overlooked in both in vivo and in vitro studies (Schmid and Cassee 2017), and thus, data from both approaches cannot be readily integrated in one toxicity classification dataset. Reporting and use of delivered dose will refine outputs, particularly where in vitro and in vivo data form base datasets for wider use in risk assessment.

The relative information gain (IG) analysis also showed a variety of other output categories which show agreement in terms of attribute importance with in vitro and in vivo studies as summarized below:

• Angiogenesis regulation responses (e.g. disruption in vascular endothelial growth factor, sprouting or angiopoietins) depend strongly on NP, exposure dose and duration (Bartczak et al. 2013; Mukherjee 2018);

• Cell adhesion and extracellular matrix responses (e.g. disrupted genes involved in cell–cell junction, morphology, movement, migration or structure) are determined by the exposure dose as also described in (Septiadi et al. 2018), NP and tissue (Ahmad Khalili and Ahmad 2015; Engin et al. 2017).

• Cell cycle and proliferation pathway responses (e.g. disrupted genes involved in cell differentiation, division, growth, mitosis progression or cytokinesis perturbation) are influenced by NP, exposure dose and cell line (0.32), in agreement with other studies (Huang, Cambre and Lee, 2017; Hussain et al. 2014).

• Cell death and apoptosis (e.g. mitochondrial fragmentation/apoptosis, phosphorylation of p53, death receptor, cleavage of caspase and PARP pathways) processes are highly determined by exposure dose and in vitro characteristics such as cell line and tissue (Ma and Yang 2016; De Stefano, Carnuccio and Maiuri 2012).

• Cytoskeleton organization (e.g. disrupted genes involved in cell matrix, microtubule, actin polymerization, cytokinesis or signal transduction) responses depend on NP and exposure dose (Kang et al. 2016; Soenen et al. 2010).

• In vitro characteristics (cell line and cell type) influence DNA damage and repair (e.g. disrupted genes involved in oncogenesis, histone or chromatin binding, replication, base excision, mismatch, chromosome damage or segregation) responses that depend strongly on exposure dose (Singh et al. 2017).

In terms of accuracy of prediction, BN-PD had marginally lower accuracy predicting metal ion responses compared to BN-RD in internal and external validation. CFS with the RD revealed that metal ion response is strongly associated with in vitro characteristics (data not shown), whereas it has poor correlation in the PD case (Figure 3). PD CFS and information gain agree with regards to metal ion dependencies where NP and zeta inputs have strong influence whilst BN-RD has higher accuracy with none of these attributes as parent-inputs to the class (Furxhi et al. 2018). Biologically one would anticipate NP to be of primary importance as metal ion responses are dependent on the presence of metal ions (i.e. metal or metal oxide NPs) (Saliani, Jalal and Kafshdare Goharshadi 2016; Poynton et al. 2012). Moreover, cell line/type and tissue would be expected to show strong associations owing to varying expression of, for example thiol rich, binding proteins.

BN-PD had zero MCC score concerning DNA dam. responses and BN-K2, BN-RD had negative scores. Reliability dataset has three triggered cases for DNA dam., corresponding to ZnONPs and A549 cell line. On the contrary, training dataset has no triggered cases, to either ZnONPs or A549 cell lines. In BN-PD (Figure 4) DNA dam. is apparently affected from cell line, cell type, dose and NP. Cell line and NP information gain score was 0.59 and 0.43, respectively (Table 2). There is no link of the values of these two significant attributes to triggered DNA dam. in the training dataset that makes the reliability dataset triggered cases impossible to predict from BN-PD. More data would compensate this issue.

Several microarray protocols were followed in the studies considered yet despite this, intra-laboratory microarray results are comparable (Bammler et al. 2005; Schneider and Orchard 2011). Microarray protocol (Array node) in a real scenario does not contribute to the change of the toxicity. However, the classifier precisely predicts the outcomes of cell adhesion and extracellular matrix responses or angiogenesis responses based on array information gain (0.14 and 0.37, respectively) interrelating dose and NP with array through shape at the Bayesian structure (Figure 4).

A further challenging issue relates to inflammatory and immune responses. Within the analysis, these outcomes depended strongly on exposure dose and in vitro characteristics consequently masking the role of physicochemical properties in inflammation. This is expected when considering that different cell types (e.g. glioblastoma vs. keratinocyte) and cell lines (e.g. primary vs. immortalized) exposed to the same NP can express markedly different cytokine levels and inflammatory profiles.

NPs can cause very overt signs of toxicity (e.g. cell death), especially at higher and often non-physiological doses, that are easy to detect yet they may also cause much more muted disruption of cellular functions at the molecular level. Such small perturbations do not necessarily result in adverse effects in vitro or in vivo (Sauer et al. 2017) or may be part of a wider sequence of events whereby small perturbations can produce large, non-linear effects on a complex system such as an organ or whole organism (Gonzalez-Diaz et al. 2013). By their very nature, small perturbations require targeted forms of analysis to detect, that in itself requires knowledge of the wider adverse outcome pathway(s) for which there are still relatively few in nanotoxicology. Alternatively, high content analysis such as omics can provide the breadth and sensitivity to detect small early changes although again identifying and understanding the link between small perturbations and complex disease endpoints remains a challenge for NP. Future studies expanding BN applicability may incorporate differentiated expression gene values (fold change) to predict the NP hazard in a quantitative manner by linking genes probabilistically with apical endpoint (diseases) such as Alzheimer’s (Alexiou et al. 2017), cancer (Achcar, Brookmeyer and Hunter 1985) or diabetes (Sangi et al. 2015). Inclusion of molecular events, defining states, association assessment and integration in the BN graph are steps to be considered in future model development. To date, limited data exist for proteomics and metabolomics for NPs. Transcriptomics is a well-established technique with high intra-array reproducibility (Fröhlich 2017) but provides only a list of genes associated with a pathway. Future integration of other omics data (multi-omics) will elucidate the potential causative changes that lead to an adverse effect and will lead to an advance in the field of nanotoxicology (Hasin, Seldin and Lusis, 2017, Shin et al. 2018).

The flexibility of BNs to missing or additional generated data allows the increase of predictive ability as the data become available or data quality increases. In addition, BNs are used in data-poor environments which fits in our case of a small set of publicly available microarray gene expression human in vitro studies. Another advantage is that BNs are understandable as a result of their graphical representation, thus they are used both as predictive and descriptive models allowing a great accessibility to the tool (Murphy et al. 2016). A detailed comparison of BN advantages with other modeling approaches in performing classification may be found in (Zabinski 2017). In this study we include exposure conditions and in vitro characteristics as model inputs, integrating quantitative metrics and study characteristics to biological effects. As more experimental data becomes available, the model can be extended to provide additional predictions. Biological phenotypic outcomes, for example based on in vitro or/and in vivo testing could be incorporated to the model to link inputs with perturbations (Key Events) and phenotypic changes (apical effects) (Hasin, Seldin and Lusis, 2017; Bañares et al. 2017). Such prediction can be used to prioritize NPs for further toxicological testing (Chepelev et al. 2015; Labib et al. 2016). Our data focus on human cell lines to minimize the loss of relevant biological information. The use of cell lines does not represent the biological complexity and thereby cannot represent the full response in the whole organism (Sahu et al. 2015). However, the complexity of the in vitro model can be enhanced by using co-cultures reducing the need for in vivo models (Bouwmeester et al. 2011, Snyder-Talkington et al. 2015) or use of organo-typic models such as precision cut lung slices (Hess et al. 2016) which combines the complexity of in vivo systems with the efficiency of in vitro models and reduced animal number. It is worth to mention that in vitro models can efficiently capture the complex molecular functions altered in vivo as demonstrated in (Kinaret et al. 2017; Zhang et al. 2016).

The model proposed here could benefit from studies that report ‘no observed transcriptomics adverse effect level’ (NOTEL). For example, Bajak et al. (2015) found that exposure of 5 and 30 nm Au NPs to Caco-2 cell line for 24 and 72 h respectively did not disturb any gene. Such numeric attributes, used as a benchmark of effects, would provide a valuable output in the established model increasing the model precision and applicability to different stakeholders. Information on physicochemical properties was scarce as seen above, thus future studies should include a detail characterization of the NPs physicochemical properties that are considered relevant to hazard identification such as surface properties, chemical composition, shape, aggregation status etc. Most of the studies gathered tested NPs at an exposure duration of 24 to 48 hours. A longer exposure may show the adaptive mechanisms at the molecular level and allow the determination of adverse versus the reversible effects. Detailed considerations for in vitro testing of NPs which should increase reliability and relevance is mentioned elsewhere (Drasler et al. 2017). The integration of in vitro and in vivo studies in one dataset poses challenges beyond this study. In vivo studies were excluded from this study because a harmonized exposure metric system is needed to transpose in vivo exposure dose to a quantity compatible to the target tissue dose (biologically effective dose) of in vitro studies (Riviere 2009; Li and Huang 2008). A PBPK model can be used to estimate the concentration of NPs in target tissue (Judson et al. 2011; Li and Reineke 2012), and therefore, incorporating a PBPK model to the BN model framework described here could be investigated to combine data from in vivo and in vitro studies (in vitro-to-in vivo extrapolation). Furthermore, using PBPK model can permit the route-to-route and species-to-species (e.g., animal-to-human) extrapolations (Raies and Bajic 2016). A detailed review of available PBPK models can be found at the supplementary material in the final technical report of Nanocomput project (Worth et al. 2017).

QSTR-perturbation models have dealt also with relatively small datasets of cases, while combining multiple endpoints, experimental conditions and toxicities/activities (Concu et al. 2017; Speck-Planche 2015). There are three ways that approach could be fused to our study: The dataset in this study can be preprocessed and modeled using the QSTR-perturbation model of Concu et al. (2017); the Concu et al. (2017) dataset can go through the CFS and BN construction and preprocessing stages to be modeled by BNs; finally, a BN can be used to model the QSTR final classification function, instead of the ANN used (Concu et al. 2017).

The BN model approach, as demonstrated here has enormous potential for the advancement of prescreening of NPs as part of risk characterization as well as the wider materials development process. In a broad context, improved screening can be foreseen as the ability to screen large candidate data sets based on physicochemical data to identify possible concerns. Considering this approach, BN models may be employed to screen combination of attributes to identify which values of the attributes lead to toxicity. This could be in the form of physicochemical data generated through testing of actual materials such as is required during regulatory activities (e.g. REACH registration) or screening of candidate designs prior to the development of prototype materials. The ultimate goal here is to reduce the burden of toxicological testing, especially in higher organisms to align with the 3R’s of toxicological evaluations and improve efficiency (time and money) between concept and commercial product.

Conclusions

NPs have entered mainstream production supply chains, while at the same time the risks of these new materials for humans and for the environment remains to be fully understood, not least when considering the long term effects of chronic exposure. Early-stage nanotoxicity measurements represent a challenge, not least when attempting to predict adverse outcomes and modeling is critical to understanding the biological effects of exposure to NPs. In this study, we demonstrate a novel approach to predict the molecular effects of NPs on human cells by using a variety of inputs such as physicochemical properties, exposure conditions and in vitro characteristics. To this end, we employ a machine learning algorithm, Bayesian networks, to afford a better understanding of how these inputs can predict normal cell function disturbance after exposure to NPs in vitro. The BN successfully predicted a number of disrupted biological pathways associated with NP properties which reflect drivers of toxicity determined also experimentally, such as zeta-potential and surface coating, in vitro and in vivo. Model structure based on statistical associations between input attributes and outputs has considerable biological relevance. The results show that appropriate pre-processing techniques should be applied before implementing a BN. It allows researchers to obtain various toxicological outcomes through one prediction concept and prioritize in vitro experimental set-ups. In summary, BN-PD that fully accounts for the CFS correlations has an aggregate predictive performance of >80%, in most of the cases and outperforms the BNs that do not use either predefined structure, gap filling or data balancing. The biggest challenge in any model development is the lack of sufficient homogenous experimental data. BNs predictability improves with additional data over time. Future integration of data confirming molecular events with phenotypic changes will allow the model presented to quantitatively determine the risk posed from exposure to NPs which can then be extrapolated to humans. The model validation shows a robust and promising methodology for incorporating transcriptomic findings in hazard and risk assessment modeling framework.

Acknowledgements

Not applicable.

Disclosure statement

No potential conflict of interest was reported by the authors

Data availability statement

The data that support the findings of this study are available at (Furxhi 2019).

Additional information

Funding

This work was supported by the European Union’s Horizon 2020 research and innovation program under Grant Number 720851. Craig A. Poland was supported by the Colt Foundation. Funding was also received from H2020 Innovation In SMEs.