Prediction of human fetal–maternal blood concentration ratio of chemicals

Dataset

The logFM values representing in vivo fetal-maternal blood ratio of 55 chemicals were obtained from a previous work collecting in vivo data from 16 published studies (Takaku et al., 2015). The 55 chemicals were randomly divided into a training dataset and a test dataset with 41 and 14 chemicals, respectively. The training dataset is utilized for feature selection, cross-validation of models, and training the final model for independent tests on the test dataset. For the development of quantitative structure-activity relationship (QSAR) models, 1-dimensional (1D) and 2-dimensional (2D) descriptors including physicochemical properties were calculated for each chemical using the PaDEL-Descriptor v2.21 software (Yap, 2011). PaDEL-Descriptor has been shown to be useful for developing QSAR models for several endpoints (Takaku et al., 2015; Huang et al., 2015; Tseng et al., 2017) and currently it is able to calculate 1,875 descriptors (1,444 1D, 2D descriptors and 431 3D descriptors) and 12 types of fingerprints. In this study, a 1,444-dimensional feature vector for each chemical was generated consisting of 1,444 1D and 2D descriptors for model development. The data tables of training and test datasets are available as Tables S1 and S2, respectively.

Model development and feature selection

To construct an interpretable model for logFM prediction, a four-step feature selection method was developed to identify the most important feature set for developing multiple linear regression models. The first three steps of the feature selection are to remove features with (1) extreme values that are at least 100-fold larger than average, (2) more than or equal to 50% zero values (scarcity), and (3) small variation (less than 8 unique values). In the fourth step, the Lasso method (Tibshirani, 1996) was applied to select informative features from the remaining feature set based on a leave-one-out cross-validation (LOOCV). Subsequently, Z-score (Kreyszig, 1979) was applied to normalize the informative features for training multiple linear regression models as shown in Eq. (1). (1)${\displaystyle z=\left(x-\mathrm{\mu }\right)∕\sigma ,}$ where z, x, µand σ represent the normalized feature value, original feature value, mean of the feature and variance of the feature, respectively. Since features with scarcity and small variation are less informative descriptors that are not useful for modeling chemicals, these features were excluded from subsequent analysis. In this study, the LassoCV and LinearRegression functions of the scikit-learn package (Pedregosa et al., 2011) were utilized to implement the Lasso-based feature selection and model training, respectively.

Applicability domain

To determine the applicability domain (AD) of the developed models, a decision tree-based method was proposed to identify the structural descriptors of chemicals prone to be incorrectly predicted. We conducted a rigorous method (Tung, Lin & Wang, 2019) to adjust the AD of the developed models solely based on the training dataset and evaluate the AD based on the test dataset. A four-step algorithm was described as follows. First, the absolute error for each chemical was calculated based on the LOOCV on the training dataset. Second, the feature vectors and absolute errors of chemicals were analyzed by a classification and regression tree (CART) (Breiman, 2017) to identify the structural rules with the highest absolute errors. The DecisionTreeRegressor function of scikit-learn package was utilized to conduct the decision tree analysis. Third, the rule with the highest absolute error was recursively identified and applied to exclude chemicals out of the AD. Finally, the adjusted AD was applied to identify chemicals in the test dataset within the AD for calculating the performance of the developed models.

Model development

To develop a prediction model for logFM, a four-step feature selection method was applied to identify informative features. The 1,444 features were firstly analyzed following the first three feature selection steps. A total of 18, 507 and 22 features with extreme values, scarcity and small variation, respectively, were excluded from the subsequent study. The Lasso feature selection was then applied to identify two informative features of AATSC1c and ZMIC1 representing the average centered Broto-Moreau autocorrelation (lag 1, weighted by charges) and Z-modified information content index (neighborhood symmetry of 1-order), respectively, from the remaining 897 features based on the training dataset.

The two features were utilized to train and cross-validate a model based on the training dataset. The fitting results of logFM for chemicals in the training dataset is shown in Table 1 and the leave-one-out cross-validation (LOOCV) results are shown in Fig. 1. Please note that the feature values shown in Table 1 are z-score normalized values. Their corresponding mean and variance are 1.32E-04 and 8.35E-06 for AATSC1c and 4.65E+01 and 3.33E+02 for ZMIC1, respectively. The developed model is shown in the following Eq. (2). (2)${\displaystyle \mathrm{logFM}=-0.0882\times \mathrm{AATSC1c}-0.2139\times \mathrm{ZMIC1}-0.3161.}$ The negative coefficients in Eq. (2) indicate a negative correlation between the two features and the logFM, i.e., an increase of the two features of AATSC1c and ZMIC1 result in a decrease of the logFM value. To further clarify the importance of the features, standard partial regression coefficients were calculated based on a piecewiseSEM package (Lefcheck, 2016). The partial regression coefficients for AATSC1c and ZMIC1 are −0.2586 and −0.6297 suggesting that a full shift in AATSC1c and ZMIC1 would result in a shift of 26% and 63% along the range of logFM. The correlation coefficient values of model fitting and LOOCV on the training dataset were 0.796 and 0.759, respectively. The small difference between the correlation coefficients of model fitting and LOOCV indicates a small chance of overfitting problems that is consistent with a similar correlation coefficient value of 0.808 obtained from the independent test on the test dataset. Detailed prediction results on the test dataset are shown in Table 2 and Fig. 2 presents the plot of observed and predicted logFM values. Please note that the feature values shown in Table 2 are z-score normalized values.

The y-randomization test, a widely used method for assessing the quality of a developed model by comparing the model performance with random models (Rücker, Rücker & Meringer, 2007), was also utilized to test the model. A total of 100 runs of LOOCV were conducted based on 100 modified training datasets whose corresponding logFM values were randomly shuffled. The mean and standard deviation of the y-randomization test are −0.202 and 0.288, respectively. The y-randomization performance is much lower than the original model (0.759) showing the uniqueness of our model.

Adjustment of applicability domain

While the model is with acceptable performance, it is extremely important to determine the applicability domain (AD) as the number of chemicals in the training dataset is relatively small compared to the huge chemical space. By properly adjusting the AD, the chemicals within the defined AD of the developed model can be identified and the prediction performance is expected to be improved (Tung, Wang & Wang, 2018; Tung, Lin & Wang, 2019). To accurately evaluate the usefulness of AD for predicting unseen chemicals, the AD of the developed model was adjusted by using only the training dataset and independently tested by using the test dataset.

In this study, we proposed a novel AD adjustment method based on the CART algorithm to identify exclusion rules for removing chemicals out of the AD. While the model performance can be continuously improved by appending more exclusion rules, a decreased coverage of chemicals within the AD can limit the practical use of the model. To balance the tradeoff between performance and coverage, we applied a stopping criterion for determining the optimal number of exclusion rules as follows. The exclusion rules were iteratively selected until no significant improvement (<1%) on the correlation coefficient was obtained by an additional exclusion rule. A total of three exclusion rules (Table 3) were selected to exclude chemicals out of AD. The exclusion rules shown in Table 3 indicate that chemicals with a relatively large AATSC1c value (0.782 <x) and medium ZMIC1 values (−0.359 <x <= 0.073) are out of AD. In addition, chemicals with small values of AATSC1c (−1.957 <x <= −0.838) and ZMIC1 (−1.315 <x <= −1.141) are out of AD. Please note that the exclusion rules are based on normalized feature values.

After the adjustment of AD, three corresponding chemicals in the training dataset were identified to be out of AD as shown in Table 1. Mifepristone, didanosine, and diazepam were excluded based on rules of #1, #2, and #3 (Table 3), respectively. The model performance was largely improved with correlation coefficient values of 0.875 and 0.850 for model fitting and LOOCV on the training dataset, respectively. The coverage of chemicals within AD is 92.68% (38/41). The defined AD was subsequently applied to exclude chemicals in the test dataset. As shown in Table 2, two chemicals were identified to be out of AD. Duloxetine and midazolam were excluded based on rules of #1 and #3, respectively. The test performance was substantially improved with a correlation coefficient of 0.847 and coverage of 85.71% (12/14) on the test dataset.

Analysis of informative features

Two features of AATSC1c and ZMIC1 were identified as informative features. AATSC1c belongs to the autocorrelation descriptor representing the average centered Broto-Moreau autocorrelation of lag 1 weighted by charges (Todeschini & Consonni, 2009). The spatial charge descriptor is relevant to lipophilicity that is considered an important factor for placental transfer (Pacifici & Nottoli, 1995). ZMIC1 is an information content descriptor representing the Z-modified information content index (neighborhood symmetry of 1-order) (Todeschini & Consonni, 2009). ZMIC1 is correlated with the molecular branching, size and ring closure and size is another important factor for placental transfer (Pacifici & Nottoli, 1995). The selected features of AATSC1c and ZMIC1 are only moderately correlated with a correlation coefficient of 0.449 showing no multicollinearity problems. To compare with the three features (MW, Hmax and TopoPSA) identified by a previous study (Takaku et al., 2015), the correlation coefficients among the two informative features of this study and previously identified three features were calculated as shown in Table 4. The analysis showed that ZMIC1 positively correlates to MW (molecular weight) with a high correlation coefficient of 0.778. In contrast, AATSC1c negatively correlates to Hmax and TopoPSA with moderate correlation coefficients of −0.611 and −0.645, respectively.

The test performance of the proposed model using the two informative features (correlation coefficient = 0.808) is much better than the previous study using three features (0.714). As a strong correlation (correlation coefficient = 0.613) between Hmax and TopoPSA and a moderate correlation (correlation coefficient = 0.494) between MW and TopoPSA were observed, multicollinearity might be responsible for the worse performance of the previous study. The three features of MW, Hmax and TopoPSA were all excluded by the fourth step of feature selection using the Lasso method and LOOCV. The proposed algorithm coped well with the multicollinearity problem, compared with the previous study based on Akaike Information Criterion (AIC) and the whole training set. After the AD adjustment, the performance was further improved to 0.847. The two informative features are considered more relevant to logFM values.