Once data was processed, curated, and featurized, the next step involved the selection of appropriate ML models to be trained. Here, due to their excellent performance on materials data, relative simplicity, and interpretability, we limited our search to support vector machines (SVMs), random forests (RFs), and gradient boosting machines (GBMs) discussed in detail in Supporting Information, Section S1 and Section S6. SVM-, RF-and GBM-based classifiers were trained on identical datasets, their performances being evaluated across four metrics; i) accuracy, ii) balanced accuracy, iii) macro-F1 score, and iv) Mattheus Correlation Coefficient (MCC)discussed in the Supporting Information, Section S1. In case it is not clear as to why more than one evaluation metric has been employed, it is because, in several cases, in addition to accuracy, we are concerned about where the model tends to fail. For instance, in the case of organic linker toxicity, apart from the accuracy of the model, we would also want to minimize potentially fatal linkers being falsely classified as safe. In such scenarios, metrics such as macro-F1 tend to shed more insights. As discussed in Section 2, to optimize model performance, it is crucial to identify a subset of the feature space that is correlated to the target variable. Here, we utilized univariate statistical tests based on the ANOVA F-value to select the K-best features, K  {10, 20 ,…, 180}. The three models implemented on the scikitlearn library in Python were trained on all subsets to find the most optimum set of features relative to the model. Using a single five-fold stratified cross-validation approach, the SVM showed the worst performance across all metrics, with an accuracy of 71.620  0.003 %, balanced accuracy of 71.587  0.003%, macro-F1 score of 0.715  0.003, and MCC of 0.574  0.005 (Figure ). As described in Scheme 1, this model was trained on a coarse grid of hyperparameters pertaining to the regularization parameter, kernel type, and shape of the decision function. The GBM was optimized on a coarse grid of hyperparameters pertaining to the learning rate, maximum iterations, maximum leaf nodes, minimum samples required to be a leaf node, and regularization parameter; it showed the best performance on the training set with an accuracy of 83.219  0.009 %, balanced accuracy of 83.131  0.009%, macro-F1 score of 0.831  0.010, and MCC of 0.748  0.014 (Figure ). In turn, the RF was optimized on a coarse grid of hyperparameters pertaining to the number of estimators and maximum depth (refer to Supporting Information, S6); it showed comparable performance to GBM with an accuracy of 81.909  0.002 %, balanced accuracy of 81.814  0.002%, macro-F1 score of 0.8117  0.009, and MCC of 0.729  0.003 (Figure ). Due to better performance on test data and overall better interpretability, we chose the RF-based classifier for subsequent investigations. The RF-based classifier showed the best performance across all metrics on a subset of 110 features selected using the F-value test. The general trend observed was an increase in performance from 20 features to 110 features, followed by a decrease in performance up to 170 features, with a slight improvement at 180 features (Figure ). Figure visualizes the hyperparameter optimization process, suggesting that at maximum depth > 20, with number of estimators > 400, the model generally performed well, showing the best performance at a maximum depth of 40 with 700 estimators. Following this, pursuant to Scheme 1, we constructed a finer grid of hyperparameters around those reported above, to finetune the model performance. Fine grid optimization as shown in Figure suggests that the model performs best at a maximum depth of 32, with 725 estimators, although these improvements were marginal at best. As pointed out earlier, apart from accuracy, a major concern for toxicity-related models is correctly identifying molecules that are fatal (specificity) and correctly identifying molecules that are not safe (sensitivity). To investigate this, we calculated the receiver-operator characteristic (ROC) curves discussed in Section S5. In the ROC curve (Figure ), the y-axis represents the true positive rate (TPR, analogous to sensitivity), while the x-axis represents the false positive rate (FPR, analogous to 1specificity). A perfect model would have a TPR of 1 and FPR of 0, while a random guess would follow the curve denoted by the dashed lines in As highlighted in the previous section, the induced toxicity largely depends on the route of administration. To investigate this, we chose a subset of the same database, reporting the LD50 of 22,449 molecules, when administered to mice through the oral route. This data was put through the same pipeline with regards to cleaning and featurization. Data sampling was conducted to produce 3000 datapoints of each category. For an equivalent comparison, a 110-feature subset was chosen using the F-test, and a RF-based model was trained, with hyperparameters optimized using the coarse-grid methodology. The model reported an accuracy of 83.18 %, balanced accuracy of 83.06 %, macro-F1 score of 0.833, and MCC of 0.743. Notably, as shown in Figure , not only are the global feature rankings for the oral route different to the i.p. route with regards to the order of the features, but there are also 26 features that are unique to the route of administration. In the case of the oral route, most of these unique features are fragment-based, suggesting that oral toxicity has a large dependence on the presence/absence of specific molecular fragments.