NanoTox: Development of a Parsimonious In Silico Model for Toxicity Assessment of Metal-Oxide Nanoparticles Using Physicochemical Features.

Metal-oxide nanoparticles find widespread applications in mundane life today, and cost-effective evaluation of their cytotoxicity and ecotoxicity is essential for sustainable progress. Machine learning models use existing experimental data and learn quantitative feature-toxicity relationships to yield predictive models. In this work, we adopted a principled approach to this problem by formulating a novel feature space based on intrinsic and extrinsic physicochemical properties, including periodic table properties but exclusive of in vitro characteristics such as cell line, cell type, and assay method. An optimal hypothesis space was developed by applying variance inflation analysis to the correlation structure of the features. Consequent to a stratified train-test split, the training dataset was balanced for the toxic outcomes and a mapping was then achieved from the normalized feature space to the toxicity class using various hyperparameter-tuned machine learning models, namely, logistic regression, random forest, support vector machines, and neural networks. Evaluation on an unseen test set yielded >96% balanced accuracy for the random forest, and neural network with one-hidden-layer models. The obtained cytotoxicity models are parsimonious, with intelligible inputs, and an embedded applicability check. Interpretability investigations of the models identified the key predictor variables of metal-oxide nanoparticle cytotoxicity. Our models could be applied on new, untested oxides, using a majority-voting ensemble classifier, NanoTox, that incorporates the best of the above models. NanoTox is the first open-source nanotoxicology pipeline, freely available under the GNU General Public License (https://github.com/NanoTox).

Introduction

Nanotechnology has delivered the promise of “plenty of room at the bottom” with transformative applications for human welfare.[1] The distinctive properties of nanoscale materials have been indispensable in industrial and medical applications, including the delivery of biologically active molecules and development of biosensors for human health and disease.[2] Engineered metal-oxide nanoparticles are characterized by a concentration of sharp edges and lend themselves to a variety of uses (e.g., ref (3)). However, there is a potential caveat to nanobiotechnology: the differential nanoscale behavior of nanomaterials might also result in emergent toxic side effects in the biological domain and ecological realm.[4−7] These hazards are related to the capacity of nanomaterials to engender free radicals in the cellular milieu, which inflict damaging oxidative stress. Such events could trigger inflammatory responses, which could balloon out of control, leading to apoptosis and cytotoxicity[8−11] as well as genotoxicity.[12]

The mundane use of nanoparticles has necessitated vigorous safety assessment of toxicity, in the interests of sustainable progress.[13−16] Such methods could also help discern safe-by-design principles that could guide adjustments to the nanoparticle formulation and thereby mitigate adverse effects at the source. Intelligent and alternative testing strategies could accelerate rational design of nanoparticles for optimal functionality and minimal toxicity.[17−20] Various computational methods have been applied to predicting toxicity of engineered nanomaterials,[21−31] but with the accumulation of high-quality data, machine learning methods have shown the most promise.[32] Such techniques provide a noninvasive “instantaneous” readout of nanoparticle toxicity[33−35] and originate from the evolution of quantitative structure–activity relationship (QSAR) models.[36] Machine learning models of nanoparticle toxicity have tended to be either generalized[37] or tissue-specific[38,39] and are built from experimental toxicity data that have been scored, standardized, and curated into databases like the safe and sustainable nanotechnology db (S2NANO).[40−42]

Earlier studies have tended to neglect systematic multicollinearity among the predictor variables, which would lead to confounding and data snooping. Second, gross imbalance between the numbers of nontoxic and toxic instances usually exists, which could lead to overfitting to the “nontoxic” class.[43] Third, we were motivated to develop a model that would be agnostic of in vitro characteristics, such as cell line, cell type, and assay method. A truly general model of nanoparticle cytotoxicity, independent of in vitro factors, would lead to significantly broader interpretability and wider applicability.[44] Our study departs also from the notion that tissue-specific models are superior to generalized models[39] and demonstrates that model interpretability is best achieved using a minimal nonredundant feature space, consistent with Occam’s parsimony. We have deployed insights from our study into a majority-voting ensemble classifier, with a view to increasing reliability. Finally, the end-to-end pipeline of our work, including the ensemble classifier, is made freely available as a user-friendly open-source nanosafety prediction system, NanoTox, under GNU GPL (https://github.com/NanoTox). All implementations were carried out in R (www.r-project.org).

Methods

Problem and Dataset

In vitro parameters such as cell type, cell line, cell origin, cell species, and type of assay could be extraneous to modeling the intrinsic hazard posed by a nanoparticle to cellular viability and the environment. This motivated us to formulate the problem in a feature space devoid of biological predictors. The machine learning task is stated as: given a certain nanoparticle at a certain dose for a certain duration, would its administration prove cytotoxic? To address this problem, we used a hybrid dataset building on the physicochemical descriptors and toxicity data found in Choi et al’s study.[36] All in vitro features were removed from the dataset, as noted above. Extrinsic physicochemical properties, namely, dosage and exposure duration, were retained.[45] The periodic table properties of metal-oxide nanoparticles published in Kar et al.[46] were used to augment the dataset. Only complete cases were considered in the process of matching the two datasets. This process yielded a final dataset of 19 features of five metal-oxide nanoparticles: Al2O3, CuO, Fe2O3, TiO2, and ZnO (Table ). Cytotoxicity was used as the outcome variable, encoded as “1” (true) if measured cell viability was <50% with respect to the control, and “0” (false) otherwise. The novel dataset is available on NanoTox.

Table 1

Physicochemical Features of MeO Nanoparticles Considered in Our Study

s no	type of feature	feature	shorthand
1	intrinsic physicochemical properties	core size	CoreSize
2		hydrodynamic size	HydroSize
3		surface charge	SurfCharge
4		surface area	SurfArea
5		conduction band energy	Ec
6		valence band energy	Ev
7		standard enthalpy of formation	Hsf
8		Mulliken electronegativity	MeO
9		enthalpy of formation of cation	enthalpy
10		polarization ratio	ratio
11	periodic table properties	pauling electronegativity	Eneg
12		summation of electronegativity	esum
13		molecular weight	MW
14		number of oxygen atoms	NOxygen
15		number of metal atoms	NMetal
16		ratio of esum to Noxygen	esumbyo
17		oxidation state	ox
18	extrinsic physicochemical properties	exposure time	Time
19		dosage	Dose

Elimination of Multicollinearity

A nonredundant feature space would translate into an optimal hypothesis space. A simple inspection of the properties in Table suggested the existence of correlated features. Correlated features would adversely impact model performance as well as complicate model interpretation. Multicollinearity is an even deeper problem in the pursuit of a nonredundant feature space.[47] The training set alone was used for the feature selection process, to prevent any data leakage from the test set. The dataset was randomly split into a 70:30 train/test ratio stratified on the outcome variable.[48] The existence of highly correlated (Pearson’s ρ ≥ 0.9) variables was ascertained. To address multicollinearity, we used a systematic variance inflation factor (vif) analysis. Each independent variable was regressed on all of the other independent variables in turn, and the goodness of fit of each of these models (fraction of variance explained; R2) was estimated. The vif score for each independent variable was then calculated using eq . In each iteration of the vif analysis, the variable in the current set that had the largest vif score when regressed on all of the other variables was eliminated. This process was continued until a set of variables all of whose vif scores <5.0 was obtained. Note that a vif score of 1.0 is possible only when a variable is perfectly independent of all other variables (all pairwise Pearson’s ρ identically zero).

Feature Transformation

The feature space could be vulnerable to heteroscedastic effects, given the varying scales for the variables. It is necessary to preprocess and prevent features with large variances from swamping out the rest. Positively skewed features could be stabilized using the log transformation. Ec values, which are negative, were first offset by +6.17, then log-transformed. Dosage spanned many orders of 10 and was log10-transformed. Exposure time spanned 2 orders of magnitude, so we performed a log2 transformation. Surface charge whose values could be either positive or negative was standardized (i.e., Z-transformed). All of the other features were log-transformed (to the base e).

Class Rebalancing

The cost of missing a toxic instance is manifold higher than the cost of missing a nontoxic instance, and the imbalance between toxic vs nontoxic instances could exacerbate this problem. In such situations, where the essential problem is to learn the minority outcome class effectively, resampling techniques could be useful.[49] We addressed the class skew problem using Synthetic Minority Over-Sampling TEchnique (SMOTE).[50] SMOTE synthesizes new minority samples from the existing ones, without influencing the instances of the majority class, thereby increasing the number of “toxic” instances relative to the number of nontoxic instances. Balancing the dataset thus would normalize the learning bias arising from unequal representation of the outcome classes.

Predictive Modeling

The overall workflow of our approach is summarized in Figure . The normalized training dataset was balanced using SMOTE, and a variety of classification algorithms were tried and tested, namely, logistic regression,[51] random forests,[52] SVMs,[53] and neural networks.[54,55]Table shows the classifiers and their hyperparameters considered in our work. The optimal values of the hyperparameters were found using 10-fold internal cross-validation.[56] The performance of each optimized model was evaluated on the normalized and unseen test set. To penalize false positives and false negatives equally, we used an objective measure of performance

Figure 1

Workflow of the study up to predictive modeling. Preprocessing refers to both normalization and class balancing. Only the training dataset was used for feature selection; the test set was kept invisible during the model development process.

Table 2

Classifiers Used in Our Study and Their Respective Hyperparametersa

no.	classifier	type/Basis	package/function	hyperparameters	optimization
1	logistic regression	algebraic	glm	threshold (= 0.5)	n/a
2	random forest	rule-based	randomForest	1 #trees (= 500)	caret::train
				2 mtry
3	support vector machine	geometric	e1071	1 Kernels (linear, radial, polynomial)	e1071::tune
				2 cost
				3 γ
				4 degree
4	neural networks	connectionist	RSNNS	1 #hidden layers = 1,2	caret::train, caret::mlpML
				2 size of each hidden layer
				3 decay rate

mtry represents the number of features used for each split in the random forest model.

Applicability Domain

The specification of the applicability boundaries of machine learning models would increase their reliability and utility.[44] This would define the perimeter of model generalization to new instances and safeguard against application to atypical data. We used a Euclidean nearest-neighbor approach to define the applicability domain (AD) of the machine learning models.[57] For each instance in the training set, its distances to all of the other training instances were found. The nearest neighbors of each instance are then defined as the k smallest values from this set, where k is an integer parameter set to the square root of the number of instances in the training set. The mean distance of an instance to its k-nearest neighbors is found, and this process is repeated for all instances to yield the sampling distribution of these mean distances. The mean and standard deviation of this sampling distribution were designated as μ and σ, respectively. The applicability domain is then defined as followswhere z is an empirical parameter (related to the z-distribution) that characterizes the width of belief in the model, which is here set to 1.96.

Results

Our dataset consisted of 483 instances of the five metal-oxide nanoparticles with 19 features and one outcome variable. Correlogram plots identified the existence of high correlation among these 19 variables (Figure ) and especially among the periodic table properties (Figure S1). Three clusters of high correlation were revealed: one cluster of enthalpy, Hsf, ratio, ox, Noxygen, and esumbyo; a second cluster of Ec and Ev; and a third cluster of esum, NMetal and MW. Based on the vif analysis, we were able to obtain a feature space of just nine uncorrelated nonredundant variables (Table ). The highest vif of any variable in this feature space was <2.02, indicating little residual multicollinearity (Figure ). This optimal feature space included two periodic table properties (Eneg, NOxygen), five other intrinsic physicochemical properties (CoreSize, HydroSize, SurfArea, SurfCharge, Ec), and both the extrinsic physicochemical properties (Dose, Time). This final dataset of 483 instances with nine features and one outcome variable is available at NanoTox.

Figure 2

Correlogram of the 19 features. The correlation between a row feature and a column feature is shown by a dot in the corresponding cell. The size of the dot represents the magnitude of the correlation, and color represents the sign of the correlation—blue: positive; red: negative. White indicates a value near 0, i.e., independence.

Table 3

Vif Scores for the Features in the Final Reduced Seta

s. no.	feature	variance inflation factor
1	CoreSize	1.65
2	HydroSize	1.24
3	SurfCharge	1.85
4	SurfArea	1.58
5	Ec	1.50
6	time	1.19
7	dose	1.21
8	Eneg	2.02
9	NOxygen	1.60

The maximum vif score is ∼2.0, corresponding to maximum R2 ∼ 0.5 (cf. eq ).

Figure 3

Optimal hypothesis space. A correlogram of the optimized feature space shows that no subset of variables in this set would yield multicollinearity.

The nine features were normalized, producing acceptable skew values for HydroSize, SurfArea, Ec, and Time (Table ). The normalized dataset was partitioned using a random 70:30 split stratified on the outcome variable, providing a training dataset of 339 instances (with 55 toxic instances), and an independent test dataset of 144 instances (with 23 toxic instances). The training dataset (and not the test dataset) was balanced for the minority toxic instances using SMOTE resampling, yielding 165 toxic and 220 nontoxic instances, for a training dataset of 385 instances. This normalized and balanced dataset was used to train the various classifiers. The optimal hyperparameters of each classifier were determined using the R e1071 package for SVMs (Figure S2), and the R caret package for the neural networks, both one layer (Figure ) and two layers (Figure S3). The full set of model-wise optimal hyperparameters could be found in Table S1. The trained, optimized classifiers were then evaluated on the unseen test dataset. All of the models, except the SVM with polynomial kernel, achieved perfect sensitivity to the toxic instances, i.e., all cytotoxic nanoparticles were classified correctly. The models were not perfectly specific to the nontoxic instances, however. On this count, the random forest and neural network one-layer models outperformed all of the others. They were each frustrated by eight false positives, resulting in a balanced accuracy of 96.69%. Bootstrapping the test set 500 times yielded standard errors of ∼0.0189 for both the random forest and neural network one-layer models, indicating performance robustness. All of the classifiers achieved balanced accuracy >90%. Table summarizes the performance of all of the models on the test set. Five nontoxic instances were classified incorrectly by all of the models, representing refractory instances and constituting a challenge to perfect learning. One of these instances was only marginally viable (0.52), indicating the possible source of refractoriness.

Table 4

Dataset Normalizationa

features	skewness before	type of normalization	skewness after	range (min–max)
CoreSize	0.92	log	–0.2	2.01–4.82
HydroSize	1.76	log	0.14	4.30–7.52
SurfCharge	0.45	z-score	0.45	–1.62 to +1.98
SurfArea	2.14	log	–0.23	1.95–5.35
Ec	2.68	log, with offset	–0.23	0.00–1.54
Time	1.36	rescale, log	–0.48	0.00–4.58
Dosage	1.74	log10	–1.5	–5.00 to +2.48
Eneg	1.46	log	1.26	0.43–0.64
Noxygen	0.66	none	0.66	1–3

Log-transformation was performed to the base e. Skewness was controlled, and the range of all predictors was brought into the same order of magnitude.

Figure 4

Hyperparameter tuning, for the neural network—1 layer model. It is seen that the cross-validation accuracy is sensitive to the choice of the set of hyperparameters.

Table 5

Performance of the Various Modelsa

		train set			test set
id	classifier	accuracy	balanced accuracy	cross-valid accuracy	accuracy	balanced accuracy
model_1	logistic regression	0.94	0.94	0.93	0.91	0.95
model_2	random forest	0.98	0.98	0.94	0.94	0.97
model_3a	SVM-Linear	0.94	0.95	1	0.9	0.94
model_3b	SVM-Radial	0.94	0.94	1	0.86	0.92
model_3c	SVM-Poly	0.98	0.98	1	0.84	0.85
model_4a	neural network—1L	0.96	0.96	0.94	0.94	0.97
model_4b	neural network—2L	0.96	0.95	0.95	0.91	0.95

Models with balanced accuracy >94% are highlighted.

Deployment

The applicability domain was calculated with the normalized train data, prior to SMOTE balancing. Substituting k = 19 and z = 1.96 in eq yielded the AD threshold = 2.23. About 95% of the test instances (i.e., 137/144 instances) were located within the AD radius. It must be noted that the misclassified instances did not coincide with these outliers. We have provided a workflow, deployment.R (available at NanoTox), for prediction on new, untested oxides. The prediction is executed by a majority-voting ensemble classifier,[58] since bagging the predictions of the best models on the test set improved the performance to just five false positives (∼98% balanced accuracy). Any new instance for classification supplied by the user is preprocessed (normalized), and its “typicality” is determined by calculating its distances to the instances in the original train data and finding the mean, D, of the 19 closest distances. If D is greater than the AD threshold, then the instance is deemed atypical for requesting the ensemble model. Predictions are obtained using the top two models, the random forest and the neural network one layer, and a consensus prediction is sought. In the absence of a consensus, an ensemble of the top five classifiers, all with balanced accuracy >94% (highlighted in Table ), is used. In the end, the majority prediction of the ensemble classifier is the predicted cytotoxicity of the given instance. Deployment.R automates this pipeline for a batch of new, untested oxides of any size. Furthermore, the RDS images of all of the models trained in our study are provided on NanoTox, for the interested scientist.

Discussion

The results are encouraging since the test set constitutes an independent validation dataset. It is clear that SMOTE balancing made a difference in the ability of the classifiers to detect the under-represented toxic instances. Filtering based on applicability domain and use of an ensemble classification strategy further mitigate model uncertainty given the ‘no free lunch’ theorem.[59] Benchmarking our results with Choi et al.,[37] we see that the best model in each classifier from our work outperformed the corresponding best models of their work (Table ). The overall best models in our work (random forest and neural network one layer) yielded a balanced accuracy of ∼97% compared to 93% for their best overall model (“neural networks”). All of the five models from this work with balanced accuracy >93% are deployed in an ensemble classifier to further mitigate uncertainty in prediction.

Table 6

Benchmarkinga

	balanced accuracy (%)
model	Choi et al.b	present work
logistic regression	92	94.63
random forest	91	96.69
SVM	91	(a) 94.21
		(b) 91.74
		(c) 85.21
neural networks	93	(a) 96.69
		(b) 94.63

SVM (a), (b), and (c) correspond to linear, radial, and polynomial kernels. Neural networks (a) and (b) refer to one and two hidden layer(s), respectively. No information regarding model hyperparameters were available in Choi et al. The best-performing models from our work are highlighted.

Ref (37).

Measures of variable importance are central to mechanistic insights.[60] Variable importance was assessed using the varImp caret function for the logistic regression model (Figure S4a), neural network one-layer model (Figure S4b), and random network model (Figure a). Dose emerges as the consensus key attribute for prediction; however, there are subtle ranking differences among the different models. NOxygen is a key attribute in both the random forest and neural network one-layer models, but not so for the logistic regression model. Time emerges as another consensus key attribute in all of the models. Logistic regression provides us with not only the effect size (coefficients) of the individual variables but also an estimate of their significance, in terms of the p-value of the coefficients (Table S2). The sign of the coefficient of each variable indicates the class outcome to which the respective variable contributes. It is notable that the two periodic table properties (Eneg, Noxygen) and the quantum chemical property, Ec, show large effect sizes but poor significance, while all of the other variables remain highly significant. Relative importance plots of the neural network models add a direction representing the favored binary outcome[61,62] and obtain concurrence to these findings (Figures b and S5). Dose emerges as the key variable determining nanoparticle toxicity, and Time, HydroSize, and Eneg are the other variables influencing the toxic prediction. NOxygen emerges as the key predictor influencing the nontoxic prediction, and SurfArea, Ec, and CoreSize are the other predictors in this category. The numeric variable importance scores are given in Tables S3 and S4.

Figure 5

(a) Normalized variable importance for the Random Forest model computed with caret. Dose is by and far the attribute with the greatest effect on the toxicity in the Random Forest model. (b) Relative importance plot for the NeuralNet-1L. Positive values correspond to the “true” (i.e., toxic) class, and negative values correspond to the nontoxic class. It is seen that Dose and NOxygen exert the maximum importance on the outcome class, though in opposite directions.

NeuralNetTools was used to visualize the best-performing one-layer neural network model, with the individual connections weighted by their importance[63] (Figure ). The two-layer neural network model was also visualized (Figure S6). Consensus among the models is necessary for explainable AI,[64] and in this direction, we performed a Lek sensitivity analysis with the neural network one-layer model.[65] How does the response variable change with changes in a given explanatory variable, given the context of the other explanatory variables? On investigating the effect of one explanatory variable, all of the other explanatory variables are clustered into a specified number of lakes with like members. While the unevaluated explanatory variables are held constant at the centroid of one lake cluster, the explanatory variable of interest is sequenced from minimum to maximum in 100 quantile steps, with the response variable predicted at each step, yielding a sensitivity curve. This process is iterated for each lake of the unevaluated explanatory variables, yielding the sensitivity profile of the response variable with respect to the specific explanatory variable in the context of the unevaluated explanatory variables. We set the number of clusters to 10, to visualize a sufficient number of response curves for each explanatory variable. In this way, the sensitivity profiles of the response variable are obtained for each predictor (Figure ).

Figure 6

Schematic of the trained neural network one-layer model, with the weights of the connections indicated by the linewidth. Black lines indicate positive weights, and gray lines indicate negative weights. Two bias units are seen, one for the hidden layer and the other for the output layer.

Figure 7

(a) Lek sensitivity analysis of attributes with positive effect on the outcome class. The steep effect of Dose is evident, with the location of the tipping point moving slightly with the cluster of the unevaluated variables. Increasing exposure times and HydroSize are also seen to tip to toxicity. (b) Lek sensitivity analysis of attributes with relatively consistent negative effect on the outcome class: CoreSize, Ec, NOxygen, and SurfArea. The number of lakes of the unevaluated variables is set to 10 in both the cases.

The two input variables that decisively differentiate the outcome are Dose and Noxygen. Dose appears to exert a nearly thresholding effect on the toxic class. The consistent sigmoidal effect seen in the “dose–response” curve, independent of the lake of unevaluated explanatory variables, echoes the maxim attributed to Paracelsus, “The dose makes a thing poison.” The attributes influencing toxicity also included: (i) Time, with a pronounced effect depending on the lakes of the unevaluated variables; and (ii) HydroSize, with a steady nonlinear effect on toxicity that is also sensitive to the context of the unevaluated explanatory variables. The response profile for Eneg is almost flat at all lakes, indicating little to no effect in changing the outcome. The interpretation of the response with respect to SurfCharge remained obscure. NOxygen emerged as the attribute with the clearest inverse effect on toxicity, with a response profile displaying a tipping point to nontoxic class at most, but not all, of the centroids. Other attributes seen to dial down the toxicity include SurfArea, CoreSize, and Ec. These observations of effect size may be tempered with significance analysis toward a complete understanding.

In summary, the ML models of our work are represented by a purely numeric feature space of just nine predictors, and it is possible to consider them in their entirety, similar to the interpretability of a classical QSAR model. The models conform to the Findable, Accessible, Interoperable, Reusable (FAIR) principles and are presented in a unified ensemble prediction engine, NanoTox (https://github.com/NanoTox). In the interest of reproducible research, all the scripts necessary to replicate, apply, and extend our analysis are available at NanoTox. Our methods may be extendable to other classes of engineered nanomaterials requiring urgent, sustainable, and rapid hazard estimation prior to induction in practical uses.[66−69]

Conclusions

We have optimized the problem formulation of cytotoxicity modeling of nanoparticles using a principled approach agnostic of in vitro characteristics. The feature space is trimmed for multicollinearity, tunable hyperparameters were optimized, and the training data were corrected for class imbalance. These steps led to an optimal hypothesis space, thereby improving the performance of the generated ML models to >96% balanced accuracy. The benefits of a parsimonious approach to modeling nanoparticle toxicity include enhanced model interpretability and generalizability. We have embedded our models into an unambiguous ensemble classifier that surpasses ∼98% balanced accuracy. Our entire workflow is available as a free open-source resource for use and enhancement by the scientific community toward proactive noninvasive testing and design of nanoparticles for varied applications.