In Silico Prediction of Human Intravenous Pharmacokinetic

Subscriber access provided by Nottingham Trent University

Pharmaceutical Modeling

In Silico Prediction of Human Intravenous Pharmacokinetic Parameters with Improved Accuracy Yuchen Wang, Haichun Liu, Yuanrong Fan, Xingye Chen, Yan Yang, Lu Zhu, Junnan Zhao, Yadong Chen, and Yanmin Zhang J. Chem. Inf. Model., Just Accepted Manuscript • DOI: 10.1021/acs.jcim.9b00300 • Publication Date (Web): 12 Aug 2019 Downloaded from pubs.acs.org on August 13, 2019

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 35 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Information and Modeling

In Silico Prediction of Human Intravenous Pharmacokinetic Parameters with Improved Accuracy Yuchen Wanga,1, Haichun Liua,1, Yuanrong Fana, Xingye Chena, Yan Yanga, Lu Zhua, Junnan Zhaoa, Yadong Chena,* and Yanmin Zhanga,* a School

of Science, China Pharmaceutical University, 639 Longmian Avenue, Nanjing, Jiangsu 211198, China. 1 The authors contributed equally to this work and should be considered as co-first authors

ABSTRACT: Human pharmacokinetics is of great significance in the selection of drug candidates, and in silico estimation of pharmacokinetic parameters in the early stage of drug development has become the trend of drug research owing to its timeand cost-saving advantages. Herein, quantitative structure–property relationship studies were carried out to predict four human pharmacokinetic parameters including volume of distribution at steady state (VDss), clearance (CL), terminal half-life (t1/2), and fraction unbound in plasma (fu), using a dataset consisted of 1352 drugs. A series of regression models were built using the most suitable features selected by Boruta algorithm and four machine learning methods including support vector machine (SVM), random forest (RF), gradient boosting machine (GBM), and XGBoost (XGB). For VDss, SVM showed the best performance with R2test = 0.870 and RMSEtest = 0.208. For the other three pharmacokinetic parameters, the RF models produced the superior prediction accuracy (for CL, R2test = 0.875 and RMSEtest = 0.103; for t1/2, R2test = 0.832 and RMSEtest = 0.154; for fu, R2test = 0.818 and RMSEtest = 0.291). Assessed by 10-fold cross validation, leave-one-out cross validation, Y-randomization test and applicability domain evaluation, these models demonstrated excellent stability and predictive ability. Compared with other published models for human pharmacokinetic parameters estimation, it was further confirmed that our models obtained better predictive ability and could be used in the selection of preclinical candidates. INTRODUCTION Pharmacokinetics (PK) is the study and mathematical description of the relationship between the drug dose and its concentration in body fluids and tissues over time. It is one of the most important evaluation indicators for drug candidates.1 The relationship

ACS Paragon Plus Environment

Journal of Chemical Information and Modeling 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

of PK is controlled by the rates of drug absorption, distribution, metabolism and excretion (ADME). The drug dosage adjustments are necessary since variations in these physiological and bio-chemical processes may influence pharmacokinetic parameters, such as volume of distribution at steady state (VDss), clearance (CL), terminal half-life (t1/2), and fraction unbound in plasma (fu).2 For VDss, it is defined as the ratio of its dose in vivo to its plasma concentration at steady state. As an important PK property in predicting the plasma concentration for a given drug in human body and affecting the half-life of the drug, it guides the dosage regimen that clinicians should prescribe to patients.3 CL is typically the most important PK attribute which quantitates the irreversible removal of a drug from the measured matrix (typically blood or plasma) for the reason that it is a determinant of other parameters including terminal half-life, oral bioavailability, and efficacious dose.4 Of all the pharmacokinetic parameters, t1/2, the time required for the plasma concentration to fall by 50% during the terminal phase, is the most frequently reported one and depends not only on drug clearance but also on volume of drug distribution.5 As a basic physicochemical parameter for predicting human VDss,6 fu refers to the fraction unbound in plasma and the subscript “u” has the usual meaning of unbound.2 The prediction of these PK parameters helps to determine the feasibility of clinical drug dosage, and provides reference for the initial dose of the first-in-human trials when dealing with the pharmacokinetics and pharmacodynamics of candidate compounds. Although ADME assessments of candidates in the middle and late stage can improve the success rate of drug candidates in clinical study, the traditional experimental methods are expensive and time-consuming. Therefore, it is essential to develop rapid and effective drug pharmacokinetics evaluation methods in the early stage of drug development. The last two decades have witnessed the flourishing of prediction on human PK parameters for query compounds. Gleeson and his colleagues7 developed the first universal quantitative structure-activity relationship (QSAR) model for the prediction of the human and rat VDss using the combination of three statistical methodologies, Bayesian neural networks (BNN), classification and regression trees (CART), and partial least squares (PLS). The QSAR model based on human data performed well with the R2test value of 0.612 and RMSEtest value of 0.479 for an independent test set. Following by several attempts for VDss prediction,8-12 multiple linear regression (MLR) was found to achieve the best predictive ability. On the basis of a dataset of 642 drugs, a newly linear PLS model was built with the Q2


Page 2 of 35

Page 3 of 35 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


value of 0.70 for the test set.13 As for clearance prediction, small data scale (less than 100 pieces) and several statistical techniques like artificial neural networks (ANN)14, 15,

MLR11, 16-18, and back propagation neural (BPN)19 were applied with R2test ranging

from 0.646 to 0.880. Arnot et al.20 employed MLR to develop and evaluate various QSARs for t1/2 estimation using a much larger dataset consisted of 1105 human t1/2 values. Later, Lu et al.21 attempted to establish a GBM model using the same dataset, yielding better predictive accuracy (R2test = 0.820, RMSEtest = 0.555). Besides, Del et al.13 built a PLS model to estimate the human fu values with a Q2 of 0.54 for the test set. Overall, the previous statistical models employed to predict pharmacokinetic parameters obtained comparable performance to that from much more resourceintensive methods. Meanwhile they were relatively single to some extent, covering only one or two PK parameters. Based on recent release of the human intravenous pharmacokinetic database published by Lombardo and Berellini,2 we have built and validated QSAR models for t1/2, VDss, CL, and fu. The models are based on two- and three-dimensional descriptors and fingerprints using RF, SVG, GBM and XGBoost. New models outperform the previously published ones, further confirming the prediction accuracy for PK parameters. EXPERIMENTAL SECTION Modeling Overview. A workflow for the PK parameters modeling process was shown in Figure 1. Specific details of the various steps would be provided in subsequent sections.

Figure 1. Flowchart of the modeling process for four pharmacokinetic parameters (A: VDss, B: CL, C: t1/2, and D: fu) prediction.



Data Curation and Preparation. The human intravenous pharmacokinetic data on a dataset of 1352 drugs collected by Lombardo et al.2 were used in this study. Since the upper molecular-weight limit for a small molecule is approximately 900 daltons,22 which allows for the possibility to rapidly diffuse across cell membranes so that it can reach intracellular sites of action.23 82 drugs with molecular-weight over the cutoff were removed from the dataset. To ensure the consistency of quality of the rest dataset consisted of 1270 drugs (Table S1), the Simplified Molecular Input Line Entry System (SMILES) notations of these molecules were converted into chemical structures and records with null values of the corresponding parameters VDss, CL, t1/2, and fu were removed. Log-transformation for these parameters was finally utilized for keeping exclusion of leverage points to a minimum, and satisfying the assumption of constant variance of errors which is required to make inferences more reliable from a multiple linear regression model.11 Finally, the remained compounds were randomly divided into training sets and test sets by a ratio of 4:1 and the details for the dataset were listed in Table S2. Principal component analysis (PCA) was further carried out to characterize the chemical space distribution of training sets and test sets for four pharmacokinetic parameters. After reducing the feature dimension of the dataset, three new features were created to summarize the whole descriptors and fingerprints. Feature Calculation and Selection. Molecular descriptors, which are quantitative representations of physical, chemical or topological characteristics of molecules, are of great importance for constructing reliable QSAR models.24 Before calculating descriptors, all the compounds were cleaned up by the molecule wash module in Molecular Operating Environment (MOE) software25 to disconnect group metals in simple salts, keep only largest molecular fragments, deprotonate strong acids, protonate strong bases and add explicit hydrogens. Then, MOE26,27 was used to calculate 206 two-dimensional (2D) structural descriptors of 1270 drug compounds. 2D structures were further subjected to energy minimization by MOE using Merck Molecular Force Field 94 (MMFF 94), where the gradient threshold of potential energy was set to 0.1 kcal/mol/Å2. Followed by conformational search of each energy-minimized structure, the most stable conformer of each structure was selected and saved into the database to generate the common descriptors. Subsequently, 93 three-dimensional (3D) descriptors were computed. The same operation was made to calculate 1001 Mordred two- and three-dimensional descriptors. In order to obtain better insight into compound dataset, 16092 molecular fingerprints were calculated by


Page 4 of 35

Page 5 of 35 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


the PaDEL software.28 Thus, a total of 17392 molecular descriptors and fingerprints were computed. To avoid the phenomenon of over-fitting, the number of variables must be reduced. In this study, training sets selected using random data split method were applied for selecting features by employing three different methods as explained below. Low Variance Filter. Eliminate features with very small variance values (< 0.2), which could clean away variables with the same value. High Correlation Filter. Exclude one of the two descriptors containing the same information, of which the correlation coefficient is higher than 0.85. Boruta Algorithm. Proposed by Kursa et al,29 it uses a wrapper approach built around a random forest classification algorithm, which provides an intrinsic measure of the importance for each feature, known as the Z score. Since the Z score is computed by dividing the average loss by its standard deviation to evaluate the importance of the descriptors, the fluctuations of the average accuracy loss of trees in the forest is taken into account. While this score is not directly a statistical measure for the significance of the feature, we can compare it to random permutations of (a selection of) the variables to test whether it is higher than the scores from random variables, which is the essence of the implementation of Boruta algorithm. By adding randomness to a system and then collecting results from random samples of the bigger system, it is based on a more general idea that one can actually reduce the misleading impact of randomness in the original sample. Model Building. Among a multitude of available modeling machine learning methods, several techniques including support vector machine (SVM),30 random forest (RF),31 gradient boosting machine (GBM),32 and XGBoost (XGB)33 were applied, which are highly effective, robust and have been extensively used in QSAR modeling.9,11,34 On the basis of the structural risk minimization principle from statistical learning theory, SVM regression is one of the widely used machine learning methods in compound ADME properties and protein structural studies.35,36 Random forest, firstly introduced by Breiman,31 is an ensemble learning technique to improve the CART method using samples of the training data and random feature selection in tree induction. Gradient boosting machine is a family of powerful machine-learning techniques whose learning procedure consecutively fits new models to provide more



Page 6 of 35

accurate estimate of the response variable. The rationale behind GBM algorithm is to construct a new underlying learner that is most correlated with the negative gradient of the loss function and associated with the whole ensemble.37 XGBoost stands for “Extreme Gradient Boosting”, which is an advanced implementation of GBM algorithm.33 Also known as “regularized boosting” technique, XGB adds regular items to cost function to control the complexity of the model. Model Evaluation. In this work, the performance of the regression models were assessed by two statistical parameters that measure goodness of fit: coefficient of determination (R2, the square of the Pearson’s correlation coefficient, based on the line of best fit) and root-mean-square error (RMSE) as follows: n

2

R =1―

(yi ― yi)2

∑i

= 1

n

∑i

(yi ― y)2

n

∑i

RMSE =

(1)

= 1

(yi ― yi)2

= 1

n

(2)

where yi, yi and y are the experimental, predicted and mean value of chemicals, respectively. Test set was employed as an external validation to ensure that the derived model from the training set has good generalization ability. The model with the highest R2test and lowest RMSEtest was selected as the best fit for further validation.38 To evaluate the accuracy of the best model for each PK parameter prediction, one commonly evaluation criterion was the fold error (FE), as represented by Eq.(3). If the FE value is smaller than 2, the prediction tends to be successful.39 Besides, the percentage of compounds in the test set predicted within 2-fold, 3-fold, and 4-fold error from experiment value was calculated, respectively. The average fold error (AFE) of the predictive models was also used to provide a measure of bias with equal value to under- and over-prediction (Eq.(4)). FE =

{

yi

yi, yi

if yi > yi

yi, if yi > yi n

∑i

AFE = 10

lg

= 1

(3)

yi

|| yi

n

(4)

Cross Validation. To assess actual predictive power and stability of the PK models, two cross validation procedures, 10-fold and leave-one-out (LOO) cross validation, were applied to the training sets. In 10-fold cross validation procedure, the initial


Page 7 of 35 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


dataset was randomly split into 10 groups, using 9 of the folds for training and the remaining part for testing. The process was repeated 10 times, each time using different training and test sets. Sharing the similar principle, in the LOO cross validation procedure, one compound was excluded from the initial set and the model was derived based on the remaining n-1 compounds. The models were evaluated by cross-validated coefficients Q210-fold and Q2LOO, respectively, according to the same equation as Eq.(1). The higher Q2 values indicates that the model is more stable and obtains stronger internal predictive ability. Y-Randomization Test. Y-randomization, firstly proposed by Klopman and Kalos,40 has been a popular tool for the validation of the performance of QSAR models, like statistical significance, robustness, predictive ability, etc. When constructing models, some selected features are combined together to best fit the given data, there exists an increased risk of chance correlation. By repeating and deliberately breaking the connection between the response variable Y and the independent variable matrix X, Y-randomization attempts to observe the action of chance in fitting given data. For each iteration, the dependent variables are randomly ordered and a new model is constructed using the original independent variable matrix. This procedure is repeated multiple times, 50 random times in this study. The performance for each run is recorded and compared with the square correlation coefficient of 10-fold crossvalidation (Q2) of the original model. It can be expected that the random generated QSAR models have low Y-randomized Q2 values. If all models from Yrandomization have lower Q2 values than that of the original model, it means that for given dataset, it is possible to obtain an acceptable QSAR model with the current modeling method. Applicability Domain Evaluation. According to the Organization for Economic Cooperation and Development (OECD) principles, the QSAR model requires a defined applicability domain (AD) for determining its degree of generalization.41 In other words, a predictive model needs to confirm the model limitations with respect to its structural domain and response space. Generally, QSAR models only cover limited chemical space on the basis of the training set. If the query molecule falls within the AD, it means that the molecule could satisfy the model scope and the prediction of the model is reliable. Otherwise, the prediction may be unconvincing. In this study, the applicability domain was evaluated by the Williams plot, which was defined by one vertical line and two horizontal lines.42 In Williams plot, the leverage values,



Page 8 of 35

calculated by the leverage matrix (H) as Eq. (8), is used to measure the distance from the centroid of the training set. H = X(XTX)

―1 T

(8)

X

where X is the descriptor matrix, XT is its transpose matrix, and (XTX)

―1

is the inverse

of (XTX). The diagonal elements in the H matrix represent the leverage values for the molecules in the dataset. The vertical line represents the warning leverage, h*, calculated by 3*m/n, where m is the number of descriptors and n is the number of training samples.43 The two horizontal lines stand for the prediction errors, calculated as the three times of variances (±3σ).

RESULTS AND DISCUSSION Diversity Distribution. Take VDss modeling as example, 988 and 248 molecules were used for model building and validation, respectively. On the basis of the whole dataset, 17392 molecular descriptors and fingerprints were computed and further taken into PCA analysis. Based on three features created by PCA, Figure 2 showed the spatial distributions of training sets and test sets (represented by red and blue balls, respectively) for four PK parameters. As clearly displayed in Figure 2A and 2B, only a small part of red dots and blue dots were far away from the main group for training sets and test sets for VDss and CL. Compared with the whole dataset, the proportion of outliers could be ignored, meaning that the overall distributions were quite similar. Compared with the fu distributions, the distributions of training set and test set of t1/2 were relatively dispersed, for all the three components concentrated between -15 and 15 (Figure 2C). Figure 2D showed that both the training and test data points of fu were mainly concentrated in the same area, where the first component concentrated between -25 and 5, the second component concentrated between -5 and 15, and the third component concentrated between -15 and 10. It was obvious that for each PK parameter modeling, the chemical space of the test set was roughly within the scope of the training set, which ensured the reasonability of the partitioning method.


Page 9 of 35 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 2. Spatial distributions of training sets (red balls) and test sets (blue balls) for four pharmacokinetic parameters (A: VDss, B: CL, C: t1/2, and D: fu) on the basis of PCA features. Feature Selection. To minimize the regression error, it is necessary to determine the most informative features of the pruned data. In this research, Boruta algorithm was employed to exact important features for modeling of four pharmacokinetic properties. The parameters of this algorithm were set as follows: doTrace = 2 to report the attribute decisions once they are cleared, maxRuns = 200 to confirm the maximal number of importance source runs. All other parameters were remained as default values. After calculation, the features were classified as important or unimportant ones with their names and total number displaying simultaneously. For VDss, CL, t1/2, and fu modeling, 209, 134, 162, and 121 features were marked as important ones and taken into modeling step respectively (Table 1). Illustrated by the case of CL prediction, 134 selected features were grouped into three categories: 59 2D molecular descriptors, 26 3D descriptors, and 49 fingerprints. Their corresponding names were listed in Table 2. For the modeling of the other three parameters, Table S3-S5 summarized the corresponding features, respectively. Table 1. Numbers of Features Selected for Pharmacokinetic Attributes Modeling Parameter

Before

Low

High

selection

variance

correlation


Boruta

After selection


Page 10 of 35

algorithm VDss

17392

14713

774

1696

209

CL

17392

14720

775

1763

134

t1/2

17392

14683

833

1714

162

fu

17392

14757

782

1732

121

Table 2. Statistical Analysis for 134 Important Variables in CL Prediction Class 2D Molecular Descriptors Physical Properties Hueckel Theroy Descriptors Subdivided Surface Areas Atom Counts and Bond Counts

Number

Name

59 5 9

apol, bpol, FCharge, logP(o/w), logS h_ema, h_emd, h_emd_C, h_logS, h_log_dbo, h_log_pbo, h_pKa, h_pKb, h_logD SlogP_VSA0, SlogP_VSA1, SlogP_VSA2, SlogP_VSA3,

12

SlogP_VSA4, SlogP_VSA7, SlogP_VSA9, SMR_VSA0, SMR_VSA1, SMR_VSA2, SMR_VSA3, SMR_VSA5

10

a_aro, a_nN, a_nO, b_1rotN, b_double, chiral, lip_acc, lip_don, opr_brigid, b_max1len

Kier & Hall Connectivity and Kappa Shape

1

KierA3

3

balabanJ, diameter, GCUT_SLOGP_0

5

a_acc, a_acid, a_base, a_don, vsa_other

Indices Adjacency and Distance Matrix Descriptors Pharmacophore Feature Descriptors

PEOE_VSA+0, PEOE_VSA+1, PEOE_VSA+2, Partial Charge Descriptors

PEOE_VSA+4, PEOE_VSA+6, PEOE_VSA-0, 14

PEOE_VSA-1, PEOE_VSA-5, PEOE_VSA-6, PEOE_VSA_NEG, PEOE_VSA_PPOS, Q_VSA_PNEG, PC+, PC-

3D Molecular

26


Page 11 of 35 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Descriptors Potential Energy Descriptors

1

PmiX, pmiY, pmiZ, rgyr, std_dim2, vsurf_CP,

Surface Area, Volume and Shape

E_ele

19

Descriptors

vsurf_CW3, vsurf_D1, vsurf_D8, vsurf_IW6, vsurf_IW7, vsurf_IW8, vsurf_W1, vsurf_W4, vsurf_W6, vsurf_W8, vsurf_Wp2, vsurf_Wp3, vsurf_Wp5

Conformation Dependent Charge

6

ASA+, ASA-, ASA_H, ASA_P, CASA+, DASA

Descriptors

Fingerprints

49

CDK fingerprint

1

FP298

3

ExtFP286, ExtFP467, ExtFP735

1

GraphFP728

MACCS fingerprint

3

MACCSFP82, MACCSFP104, MACCSFP123

Pubchem fingerprint

3

PubchemFP346, PubchemFP366, PubchemFP528

2

SubFP297, SubFP298

CDK extended fingerprint CDK graph only fingerprint

Substructure fingerprint Substructure fingerprint count

SubFPC12, SubFPC84, SubFPC85, SubFPC88, 10

SubFPC136, SubFPC169, SubFPC287, SubFPC296, SubFPC300, SubFPC307 KRFPC297, KRFPC677, KRFPC1147, KRFPC2547,

Klekota-Roth fingerprint count

15

KRFPC2548, KRFPC3224, KRFPC3295, KRFPC3408, KRFPC3410, KRFPC3436, KRFPC3455, KRFPC3591, KRFPC3750, KRFPC3957, KRFPC4080

2D atom pairs

1

AP2D102 APC2D2_N_O, APC2D2_O_O, APC2D3_C_N,

2D atom pairs count

10

APC2D3_N_O, APC2D4_C_C, APC2D4_C_N, APC2D4_O_O, APC2D5_C_O, APC2D9_C_O, APC2D10_C_N



Page 12 of 35

Model Selection. Four machine learning methods including SVM, RF, GBM and XGB were employed for model construction using features selected by Boruta algorithm. For the prediction of four different pharmacokinetic attributes, a series of regression models were built using different training sets. Test sets were then employed to ensure the generalization abilities of the predictive models. The statistical results of the QSAR models in fit and prediction of human PK parameters were given in Table 3. Table 3. Human VDss, CL, t1/2, and fu Statistics of the SVM, RF, GBM and XGB Models on Training sets and Test sets Parameter

VDss

CL

t1/2

fu

Method

𝐑𝟐𝐭𝐫𝐚𝐢𝐧

𝐑𝐌𝐒𝐄𝐭𝐫𝐚𝐢𝐧

𝐑𝟐𝐭𝐞𝐬𝐭

𝐑𝐌𝐒𝐄𝐭𝐞𝐬𝐭

SVM

0.954

0.140

0.870

0.208

RF

0.916

0.188

0.856

0.219

GBM

0.922

0.182

0.845

0.227

XGB

0.934

0.167

0.833

0.236

SVM

0.648

0.411

0.541

0.200

RF

0.882

0.239

0.875

0.103

GBM

0.739

0.402

0.577

0.192

XGB

0.836

0.280

0.724

0.155

SVM

0.854

0.252

0.708

0.203

RF

0.878

0.230

0.832

0.154

GBM

0.867

0.240

0.754

0.184

XGB

0.855

0.251

0.758

0.184

SVM

0.822

0.336

0.772

0.326

RF

0.927

0.214

0.818

0.291

GBM

0.834

0.324

0.778

0.321

XGB

0.898

0.253

0.767

0.329

VDss. For the prediction of human PK parameter VDss, the four modeling methods shared the similar results shown in Table 3. Examination of the results in prediction of the independent test set of 248 compounds demonstrated that the SVM model was the best most predictive of the four models with R2train of 0.954 and RMSEtrain of 0.140 for the training set, and R2test of 0.870 and RMSEtest of 0.208 for the test set. RF produced the inferior prediction accuracy, yielding R2test of 0.856 and RMSEtest of 0.219. Meanwhile, the performance of RF in the training set was slightly worse than


Page 13 of 35 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


the GBM and XGB models, achieving lower R2train value and higher RMSEtrain value (for RF, R2train = 0.916 and RMSEtrain = 0.188; for GBM, R2train = 0.922 and RMSEtrain = 0.182; for XGB, R2train = 0.934 and RMSEtrain = 0.167). Figure 3 displayed the relationships between experimental logVDss values and predicted ones for four different algorithms. Almost all the data points located close to the regression line, which indicated that both the training set and test set got the similar prediction results. With nearly none data points owning large difference between predicted and experimental logVDss values, these four models achieved large R2 and low RMSE values and could predict human VDss effectively.

Figure 3. Plots of experimental human logVDss versus in silico prediction by four machine learning methods. The regression line is colored grey. The training set and test set are represented by red triangles and blue dots, respectively. CL. From the statistical results of the human CL models in training set fit, it was obvious that the RF and XGB models would predict CL in humans more effectively than GBM and SVM on account of their larger R2 and smaller RMSE values. For the training set, the RF model obtained the highest R2train value of 0.882 and the lowest RMSEtrain value of 0.239, and the XGB model ranked the second when comparing model performance (R2train = 0.836 and RMSEtrain = 0.280). The GBM and SVM models performed relatively worse but acceptable, with the R2train values larger than 0.64 and RMSEtrain values lower than 0.45. Analysis of the individual test set



prediction results would suggest that the RF algorithm performed the best among the four regression models, achieving the R2test of 0.875 and RMSEtest of 0.103. Contrary to the VDss modeling, SVM got the worst predictive ability in CL estimation and its R2 was only 0.541 for test set, 0.648 for the training set. As clearly demonstrated in the plot of predicted logCL versus experimental ones for the training set and test set (Figure 4), most of the data points in the RF subplot were concentrated to the regression line, while there existed points far away from the fitted line for the SVM subplot, corresponding to the statistics results for human CL predictive models.

Figure 4. Plots of experimental human logCL versus in silico prediction by four machine learning methods. The regression line is colored grey. The training set and test set are represented by red triangles and blue dots, respectively. t1/2. The results in Table 3 highlighted general consistencies in the statistics for the four human t1/2 models. Of all the compounds in the training set, their R2train values were higher than 0.85 and RMSEtrain lower than 0.26. When it was applied to external test set, RF performed the best with the R2test value of 0.832, the only one higher than 0.80. Compared with the statistic results of RF modeling, GBM and XGB performed relatively worse in test set prediction. Both the R2test values of the GBM (0.754) and XGB (0.758) were lower than 0.80, indicating that they shared the similar human t1/2 prediction. Among the four QSAR modes for fu prediction, the SVM model obtained the worst performance with R2test value around 0.70. Considering the performance in training set fit and test set prediction of each model, we concluded that the RF


Page 14 of 35

Page 15 of 35 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


algorithm produced superior predictive ability for human t1/2 prediction. As the most predictive one among the four QSAR models for t1/2 prediction, the scatter plot of RF model (Figure 5) was the most concentrated, which further confirmed the conclusion.

Figure 5. Plots of experimental human logt1/2 versus in silico prediction by four machine learning methods. The regression line is colored grey. The training set and test set are represented by red triangles and blue dots, respectively. fu. The human logfu modeling results displayed in Table 3 may again suggest that the RF model would produce superior predictive ability than the other three methods for fu prediction. When the RF model was applied to fit the training set, the R2train was 0.927 and RMSEtrain was 0.214, and for the prediction of test set, the R2test was 0.818 and RMSEtest was 0.291. All the statistical results were higher than that of the other models, indicating that the selected 121 features could effectively predict logfu values. The other three models owned similar performance in the test set validation, with the R2test values ranging from 0.76 to 0.78, RMSEtest from 0.32 to 0.33. By comparison of the training set statistical values, we concluded that the XGB model could predict human logfu more effectively than GBM and SVM, yielding higher R2train value (0.898) and lower RMSEtrain value (0.253). Figure 6 recorded the scatter of the predicted logfu values and experimental ones by four machine learning methods. For the RF model, the data was the most concentrated to the regression line, while the SVM subplot displayed the most dispersed data distribution.



Page 16 of 35

Figure 6. Plots of experimental human logfu versus in silico prediction by four machine learning methods. The regression line is colored grey. The training set and test set are represented by red triangles and blue dots, respectively. Model Validation. According to the statistics of test set, the most predictive model was selected separately for different PK parameters. The statistical results of four models in predicting test sets were summarized in Table 4. The SVM model for VDss prediction was selected as the best one, yielding an AFE value of 1.018 and a percentage of compounds predicted within 2-, 3-, and 4-fold of the actual value of 87%, 96%, and 99%, respectively. As for CL prediction, the RF model showed the best prediction with an overall AFE of 1.050 and the accuracy of prediction within 2fold of 100%. The t1/2 prediction showed the similar performance as CL prediction, achieving an AFE value of 1.053 and the accuracy of prediction within 2-fold of 99%. Besides, the good predictive results of fu modeling further explored the value of our methods. Table 4. Statistical Parameters on the Selected Best Models for PK Prediction of Human VDss, CL, t1/2, and fu Best Model

Q210-fold

VDss-SVM

0.770

Q2LOO 0.763

AFE

< 2-fold

< 3-fold

< 4-fold

1.018

87%

96%

99%


Page 17 of 35 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


CL-RF

0.764

0.778

1.050

100%

100%

100%

t1/2-RF

0.730

0.729

1.053

99%

100%

100%

fu-RF

0.845

0.865

0.956

81%

93%

96%

As shown in Table 4, these models were internally validated by both 10-fold cross validation and leave-one-out cross validation, gaining Q210-fold values of 0.770, 0.764, 0.730 and 0.845, respectively. For the selected models for PK prediction, Q2LOO values ranged from 0.729 to 0.865. The significant cross-validated coefficient values larger than 0.5 ensures the high predictive ability of these models.44 Y-Randomization was further carried out to validate the robustness of the selected models for pharmacokinetic properties prediction. In this validation, the performance of 10-fold cross-validation was compared with that of models built for shuffled response variables. Figure 7 showed the distributions of Y-randomized Q2 of the predictive models and the fitted curves for four different parameters. Taking the RF model for CL prediction as example, 1014 response variables were randomly permuted, keeping 134 independent variables intact. As shown in Figure 7B, the Q2 values of randomly shuffled models varied -0.101 from to -0.021, which was significantly lower than that of the real model (Q210-fold = 0.764). It was obvious that for prediction from the other three models of VDss, t1/2, and fu, the fitted curve of randomized Q2 values in each subplot was far away from that of original model, represented by the vertical line. The low Q2 values indicated that the good results in our original models were not due to a chance correlation or structural dependency on the training set.



Figure 7. Distributions of Q2 of randomized models constructed in Y-randomization test. On the left side of each subplot presents the distribution of Q2 values of randomized models while the red vertical line on the right is the Q2 value of the original model (A: VDss, B: CL, C: t1/2, and D: fu). The domain of applicability is an important concept for QSAR models which provides the estimation of uncertainty in predicting a query molecule based on the similarity with the compounds used in the modeling. In this study, the Williams plot was employed to investigate whether the molecules fall in the defined applicability domain or become outliers. Figure 8 presented the Williams plots of the four models for different PK parameters, VDss, CL, t1/2, and fu. As illustrated by the case of CL prediction, Figure 8B clearly showed that almost all the molecules in the training sets fell within the defined applicability domain (h* = 0.396). Besides, some compounds in the test set (36 compounds) fell outside the applicability domain. Compared with the total number of test set (254 compounds), prediction for 84% of the compounds were acceptable, similar with the other three models (Table S6).


Page 18 of 35

Page 19 of 35 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 8. Williams plots of leverage values versus prediction errors for training sets (marked in blue dots) and test sets (represented by red triangles) of four pharmacokinetic attributes (A: VDss, B: CL, C: t1/2, and D: fu). Of each subplot, the grey horizontal line represents the leverage value (h*) and two grey vertical lines mean the prediction errors, three times of variance (±3σ). Model Interpretation. Due to the large number of features used and the complex modeling procedure, it is a hard task to thoroughly understand these models. The analysis of important descriptors and fingerprints for modeling would provide more information about the models. Evaluated by the “Normalized permutation importance” criterion in Boruta algorithm,45 the important descriptors and fingerprints used to build models (Table 5-9) were given importance scores. Among these features, the top 10 features obtaining the highest importance scores were marked with names and taken into further analyses (Figure 9). It is noteworthy that the database utilized for constructing our models cover several drugs containing transition metals (e.g. molecule 187, and 812 etc.) and lanthanides (e.g. molecule 518, and 519 etc.), which may lead to bias to the feature interpretations. Fortunately, the removal of those compounds can also result in comparable models with only slight difference in R2 and RMSE with the models built on the whole dataset, further proving the structural inclusiveness and extensiveness of our models. The correlation between each PK parameter and their corresponding features with great significance was shown in Figure S1-S4, respectively.



Figure 9. Important variables distribution of four pharmacokinetic attributes modeling (A: VDss, B: CL, C: t1/2, and D: fu). Each column represents one variable and its height depends on the importance score calculated by Boruta algorithm. For the top 10 features with the highest scores, their names are marked next to the columns. VDss. Estimation of the apparent volume of distribution of a drug in the human body at steady state (VDss) has become a key pharmacokinetic problem, owing to its ability to predict the plasma concentration for a given amount of drug in the body and determine the half-life together with clearance. As shown in Table 5, 10 features with the largest importance scores among 209 features were further analyzed for VDss modeling. SubFP297, the only one fingerprint in the VDss prediction, is deemed as an important variable with the highest absolute value of correlation coefficient (-0.507), which emphasized the potentiality of fingerprints in representing compounds. Ranked by the importance score, the 3D descriptor, vsurf_Wp2, is the sixth valuable feature containing information about polar volume. The remaining eight descriptors are defined as 2D molecular descriptors which


Page 20 of 35

Page 21 of 35 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


further classified as different types. FCharge, logP(o/w) and logS make up the first type, physical properties. FCharge refers to the total charge of the molecule. Owing to it second highest importance score (12.333) and third largest correlation coefficient (0.451), it is selected as one important descriptor to accurately predict the logVDss values. logP(o/w) represents the log value of the octanol/water partition coefficient (including implicit hydrogens) and its correlation coefficient is 0.326, positively correlated with logVDss values, meaning that the greater the lipophicity of a molecule is, the easier it is to make the VDss prediction with larger values. logS is defined as the log value of the aqueous solubility (mol/L) of molecules. To the contrary of logP(o/w), the correlation coefficient between logS and logVDss is negative. The second type of 2D descriptors are Hueckel Theory descriptors, consisted of h_emd and h_ema. These descriptors are based on a modified Hueckel Theory calculation and do not depend on the particular input resonance form and consistently treat electron withdrawal effects. h_emd is defined as the sum of hydrogen bond acceptor strengths, while h_ema contains information about hydrogen bond acceptor strengths. The last three descriptors represent three types of different 2D descriptors, for GCUT_SLOGP_0 refers to the logP value calculated by the GCUT descriptors, a new class of molecular descriptors that encode not only topological information but also atom-based information related to ligand-receptor interaction strength.46 As a Pharmacophore Atom Type descriptors, a_acid refers to the number of acidic atoms of a molecule, considering only about the heavy atoms. PEOE_VSA-6 refers to the sum of van der Waals (vdw) surface area of atoms when |𝑞𝑖| < ―0.3. All of these three descriptors achieved high importance scores and absolute values of correlation coefficient, identifying their significance in the prediction of logVDss values. Table 5. Detailed Information on the Most Essential Variables for VDss Modeling Importance

Correlation

score

coefficient

LogP GCUT (0/3)

12.943

-0.332

2D

Sum of formal charges

12.333

0.451

a_acid

2D

Number of acidic atioms

11.732

-0.453

vsurf_Wp2

3D

Polar volume at -0.5

10.405

0.079

Rank

Name

Class

Description

1

GCUT_SLOGP_0

2D

2

FCharge

3 4



5

PEOE_VSA-6

2D

Total negative 6 vdw surface area

Page 22 of 35

10.333

-0.401

10.318

-0.507

9.995

0.129

9.943

0.326

9.898

-0.420

9.843

-0.119

Presence of SMARTS 6

SubFP297

FP

patterns for functional hroup classification by Christian Laggner

7

h_emd

2D

8

logP(o/w)

2D

9

h_ema

2D

10

logS

2D

Sum of hydrogen bond donor strengths Log octano/water partition coefficient Sum of hydrogen bond acceptor strengths Log solubility in water

CL. Drug clearance can be divided into three categories: metabolic transformation, renal excretion, and hepatobiliary excretion. A multitude of biochemical and physiological mechanisms may influence the clearance rate of each type.47 Physicochemical properties are usually the representation of the rate-limiting clearance mechanism, with lipophilic molecules tending to metabolism and hydrophilic, polar molecules tending toward passive or active excretion.48 As listed in Table 6, ASA-, PEOE_VSA-6, and other eight descriptors are of great importance exhibited the information about polar surface area (PSA), water solubility and hydrogen bond of the molecules. Among the top 10 descriptors, ASA- of a molecule refers to the water accessible surface area of all atoms with partial negative charge, and accessible surface area refers to the water accessible surface (in Å2) area using a probe radius of 1.4 Å. ASA_H is defined as the total hydrophobic surface area (|qi| < 0.2, qi is the partial charge of atom i). ASA_P denotes the surface area sum over all polar atoms (|qi| ≥ 0.2), and DASA is the absolute value of the difference between ASA+ and ASA-, where ASA+ refers to the positive accessible surface area. These descriptors are classified as 3D descriptors closely related to the partial charge of each molecule and its conformation, and contain information about water contact polar surface area. Same as VDss modeling, GCUT_SLOGP_0, h_ema and PEOE_VSA-6 are important descriptors for predicting logCL. GCUT_SLOGP_0 represents the size of the compound logP (lipid-water partition coefficient). h_logS refers to the water


Page 23 of 35 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


solubility of the compound, and h_ema contains information about the ability to form hydrogen bonds. In some studies, hydrogen bonding is considered to be an important factor affecting drug permeability.49,50 These three descriptors have influence on the clearance rate of the compound by affecting its permeability. PEOE_VSA-6 is one partial charge descriptor, depending on the partial charge of each atom of a chemical structure require calculation of those partial charges and SlogP_VSA0 means the sum of vdw surface area when Li ≤ ―0.4, where Li represents the contribution of logP(o/w) at the time of calculating the descriptor SlogP of atom i.51 In MOE, SlogP shares the same definition as the logP(o/w), which are both the representation of lipid-water partition coefficient of the compound. SlogP_VSA2 is defined as the sum of vdw surface area when Li is between -0.2 and 0. The above 2D descriptors reveal the information on the molecular surface properties, especially the polar surface area, which were found to be the calculation filters for membrane permeability in the early stages of drug discovery.52,53 Drugs with high-permeability will reabsorb in the renal tubules and cross the cell membrane at high passive transport rates, in which case kidney clearance is reduced and these compounds will be primarily cleared by metabolism. Table 6. Detailed Information on the Most Essential Variables for CL Modeling Rank

Name

Class

1

ASA-

3D

2

PEOE_VSA-6

2D

3

SlogP_VSA0

2D

4

h_ema

2D

5

h_logS

2D

6

SlogP_VSA2

2D

7

ASA_H

3D

8

DASA

3D

9

GCUT_SLOGP_0

2D

Importance

Correlation

score

coefficient

11.310

-0.219

10.573

-0.207

10.459

-0.113

10.374

-0.245

Log solubility in water

10.342

0.112

Bin 2 SlogP (-0.20, 0.00]

9.934

-0.164

9.624

-0.035

9.574

0.144

9.260

-0.211

Description Negative accessible surface area Total negative 6 vdw surface area Bin 0 SlogP (-10, -0.40] Sum of hydrogen bond acceptor strengths

Total hydrophobic surface area Absolute difference in surface area LogP GCUT (0/3)



10

ASA_P

3D

Total polar surface area

Page 24 of 35

9.122

-0.133

t1/2. 10 descriptors with importance scores higher than 8.10 were summarized in Table 7. Compared with other four models, t1/2 model is unique owing to the existence of only 2D descriptors, including 4 autocorrelation descriptors (ATSC4s, AATSC0i, ATSC4dv, AATS0v), 2 physical properties (logS and logP(o/w)), 1 partial charge descriptor (PEOE_VSA-6), 1 BCUT descriptor (BCUTv-1l), 1 adjacency matrix descriptor (VE1_A), and 1 Estate descriptor (SdssC). Autocorrelation descriptors are a class of molecular descriptors based on the statistical concept of spatial autocorrelation applied to the molecular structure.54 ATSC4s and ATSC4dv refer to centered moreau-broto autocorrelation of lag 4 with different weight. Both AATSC0i and AATS0v represent averaged moreau-broto autocorrelation of lag 0, weighting by different index, ionization potential and vdw volume, respectively. All of them are positive correlated with logt1/2, suggesting that chemicals with shorter t1/2 tend to have smaller autocorrelation values. logS and logP(o/w) are both physical properties existed in the fu and t1/2 models, proving that the basic properties ought to be used for molecule representation for modeling. PEOE_VSA-6 is an important parameter to characterize the negative surface area of chemicals and shows weakly negative correlation with logt1/2 (-0.166). BCUTv-1l, a member of BCUT metrics found to provide unique information regarding molecular structure.55 VE1_A is classified as one adjacency matrix descriptor which may positively influence the logt1/2 values. SdssC, the only Estate descriptor, has a positive contribution in t1/2 prediction. Table 7. Detailed Information on the Most Essential Variables for t1/2 Modeling Importance

Correlation

score

coefficient


10.621

-0.264

VE1 of adjacency matrix

9.889

0.242

9.634

0.282

8.888

0.182

8.675

-0.166

Rank

Name

Class

Description

1

logS

2D

2

VE1_A

2D

3

logP(o/w)

2D

Log octano/water partition coefficient Centered moreau-broto

4

ATSC4s

2D

autocorrelation of lag 4 weighted by intrinsic state

5

PEOE_VSA-6

2D

Total negative 6 vdw surface


Page 25 of 35 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


area First lowest eigenvalue of 6

BCUTv-1l

2D

Burden matrix weighted by

8.565

-0.092

8.506

0.067

8.248

0.148

8.193

0.124

8.129

0.040

vdw volume Averaged and centered moreau-broto 7

AATSC0i

2D

autocorrelation of lag 0 weighted by ionization potential Centered moreau-broto

8

ATSC4dv

2D

autocorrelation of lag 4 weighted by valence electrons

9

SdssC

2D

Sum of dssC Averaged moreau-broto

10

AATS0v

2D

autocorrelation of lag 0 weighted by vdw volume

fu. As clearly illustrated in Table 8, 2D structural descriptors almost covered the important-variable lists, with only one 3D structural descriptors, ASA-, included. By counting its frequency in the lists of four models, ASA- is deemed as an essential descriptor for CL and fu modeling, further illustrating the essentiality of the conformation dependent charge descriptors. The rest descriptors reveal information about hydrogen bond, partial charge of each atom and octanol/water partition coefficient of the molecule. h_logD represents the octanol/water distribution coefficient at pH 7, classified into the group of Hueckel Theory descriptors report quantities based on averages over all protonation states at pH 7. Obtaining both the highest absolute value of correlation coefficient (-0.617) and the highest importance score (18.967), h_logD is considered as the most important descriptor for fu prediction. h_logS, h_emd, and h_log_pbo were divided into another group of Hueckel Theory descriptors, carrying information about aqueous solubility and hydrogen bond of the molecule. As mentioned above, both logS and logP(o/w) are representation of physical properties. Calculated by different methods, logS shares the same description with h_logS and both of them rank top in Table 8, emphasizing the value of water solubility in fu modeling. Besides, a_aro refers to the count of aromatic atoms and has a strong positive relationship with logfu values owing to its



Page 26 of 35

correlation coefficient higher than 0.50. PEOE_VSA_NEG and PEOE_VSA-1 are representation of polar surface area, both of them obtain negative correlation coefficients with logfu values, meaning that when PSA increase in amount, the human logfu values will decrease. Table 8. Detailed Information on the Most Essential Variables for fu Modeling Rank

Name

Class

Description

Importance

Correlation

score

coefficient

18.967

-0.617

Octanol/water 1

h_logD

2D

distribution coefficient (pH = 7)

2

logS

2D

Log Solubility in Water

15.163

0.615

3

h_logS

2D


15.053

0.533

4

logP(o/w)

2D

14.659

-0.554

5

ASA-

3D

11.643

-0.431

6

h_log_pbo

2D

10.447

-0.455

7

a_aro

2D

10.440

-0.504

8

PEOE_VSA_NEG

2D

10.327

-0.469

9

PEOE_VSA-1

2D

9.438

-0.457

10

h_emd

2D

8.911

0.183

Log octanol/water partition coefficient Negative accessible surface area Sum of log (1+p-bond orders) Number of aromatic atoms Total negative vdw surface area Total negative 1 vdw surface area Sum of hydrogen bond donor strengths

Comparative Analysis with Published Models. For further estimating the predictability of the built models, comparative analysis was conducted with previously published models for pharmacokinetic parameters. Over the years, substantial progress has been made in the computational modeling for the human PK parameters, VDss, CL, t1/2, and fu, even in the absence of large datasets. As shown in Table 9, most of the previous models employed the MLR algorithm for modeling, while in our research, RF and SVM models were respectively selected as the best ones due to their robustness and prediction ability. Besides, the dataset quality and quantity may have a great effect on the model performance and this hypothesis was


Page 27 of 35 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


corroborated by the modeling results of VDss, CL and t1/2. As for VDss modeling, the dataset capacity of published models ranged between 121 and 642 and their R2test values varied from 0.60 to 0.82. In this study, 1236 records on human VDss were employed to build regression models and SVM achieved the maximum R2test value of 0.870. With regard to CL and t1/2 models, both of our models made use of larger dataset to construct models, 1268 and 1253, respectively, and achieved higher R2test or RMSEtest values than existing ones. Overall, compared with published models for PK parameters prediction, our models applied larger dataset for modeling and performed better in the prediction. In addition, relatively few computational studies have been undertaken on the prediction of human fu. Thus, it is of great value to predict these four important pharmacokinetic parameters including VDss, CL, t1/2, and fu. Table 9. Summary of Computational Models for Human PK Parameters Prediction PK endpoint

Year

Authors

Dataset capacity

Modeling

Model

Model

method

performance

evaluation

Descriptors

R2test = 0.612, 2006

Paul et al.7

199

Physicochemical

CARTBNN-PLS

RMSEtest =

Y-rand

0.479 2009

2009

2011 VDss 2013

2013

2015

Sui et al.8 Giuliano et al.9 Zhivkova et al.10 Gombar et al.11

Del et al. 13

Zhivkova et al.12

121

Physicochemical

MLR

R2test = 0.817

669

MOE, VolSurf+

MLR

Q2 = 0.54

MLR

R2test = 0.687

MLR

R2test = 0.782

PLS

Q2 = 0.70

MLR

R2test = 0.602

126

569

642

216

ACD/LogD, MDL QSAR

E-state

ACDlabs, VolSurf+, MOE ACD/LogD, MDL QSAR MOE,

2019

Wang et al.

1236

Mordred,

SVM

PaDEL


R2test = 0.870,

LCO Y-rand, LCO LOO, LMO

LOO

Y-rand, AD LOO, LCO

Y-rand, AD, LOO

RMSEtest =


Page 28 of 35

0.208, Q2 = 0.770

1999

2000

2001

Schneider et al.14

Ekins and Obach16 Zuegge et al.19

R2test = 0.86, 22

-

ANN

LOO Q2 = 0.79 R2test = 0.88,

18

Cerius2

MLR

LOO Q2 = 0.79

22

-

BPN

R2test = 0.84

LOO

Q2 = 0.682, 2002

Wajima et al.17

68

Physicochemical

MLR

RMSEtest

LOO

=0.350 CL

2009

2010

2013

Li et al. 18

Paulo et al.15

Gombar et al.11

49

89

525

TASR

ALOGP, EDragon

E-state

MLR

R2test = 0.73

LOO

R2test = 0.646, ANN

RMSEtest

-

=0.544

MLR

R2test = 0.70

LOO

R2test = 0.875, MOE, 2019

Wang et al.

1268

Mordred,

RF

PaDEL

RMSEtest =

Y-rand,

0.103,

AD, LOO

Q2 = 0.764

2013

Arnot et al.20

R2test = 0.73, 1105

Physicochemical

MLR

LMO RMSEtest = 0.75 R2test = 0.820,

2016

Lu et al.21

1105

Codessa, Pipeline Pilot

GBM

RMSEtest =

AD

0.555

t1/2

R2test = 0.832, MOE, 2019

Wang et al.

1253

Mordred,

RF

PaDEL

RMSEtest = 0.154, Q2 = 0.730


Y-rand, AD, LOO

Page 29 of 35 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Del et al. 13

2013

541

ACDlabs, Volsurf+, MOE

PLS

Q2 = 0.54

Y-rand, AD

R2test = 0.818, MOE,

fu 2019

Wang et al.

878

Mordred,

RF

PaDEL

RMSEtest =

Y-rand,

0.291,

AD, LOO

Q2 = 0.845

CART, classification and regression trees; BNN, Bayesian neural networks; PLS, partial least squares; MLR, multiple linear regression; SVM, support vector machine; BPN, back propagation neural; ANN, artificial neural networks; RF, random forest; GBM, gradient boosting machines; R2, squared Pearson’s correlation coefficient; RMSE, root-mean- square error; Q2: squared crossvalidated correlation coefficient; Y-rand, Y-randomization; LCO, leave-class-out; LOO, leaveone-out; LMO, leave-many-out; AD, applicability domain.

CONCLUSION The prediction of human pharmacokinetic parameters for new chemical entities is widely used in predicting dosing regimen appropriately. However, traditional experimental approaches for the assessment of important PK parameters in humans such as extrapolation of multiple preclinical PK measurements to human beings, are limited by time, cost and resource. In this work, a dataset of human pharmacokinetic parameters values for 1352 drug compounds recently published was explored and four machine learning algorithms were employed for computational prediction of human PK parameters. As for VDss, the SVM algorithm performed best, yielding R2test of 0.870 and RMSEtest of 0.208 for the test set. For the estimation of other parameters (CL, t1/2, and fu), RF showed better performance than other models. The robustness of these models was further examined using 10-fold cross-validation, leave-one-out cross-validation, Y-randomization test and defined the applicability domains. Lastly, the most important 10 features affecting each PK parameter most were successively selected for analysis to interpret the models clearly. By comparison with previously published models for human pharmacokinetic parameters, our models applied to larger dataset for modeling and performed better in the prediction. It can be concluded that our research is an important application of in silico models for the entirely prediction of four human pharmacokinetic parameters and provides useful guidance for the early-stage drug discovery.



SUPPORTING INFORMATION Further details on the models including the initial dataset, molecular descriptors, and correlation maps regarding top 10 variables (DOCX)

CORRESPONDING AUTHORS For Yanmin Zhang: Tel.: +86-25-86185163. Fax: +86-25-86185182. E-mail: [email protected] For Yadong Chen: Tel.: +86-25-86185163. Fax: +86-25-86185182. E-mail: [email protected] ORCID Yanmin Zhang: 0000-0002-4075-7556 Yadong Chen: 0000-0002-2189-9356 Author Contributions Yuchen Wang and Haichun Liu contributed equally to this work and should be considered as co-first authors.

CONFLICT OF INTEREST The authors declare no competing financial interest. ACKNOWLEDGMENTS This work was financially supported by National Natural Science Foundation of China (No. 81803370), Natural Science Foundation of Jiangsu Province (No. BK20180559), State Key Laboratory Innovation Research and Cultivation Fund (No. SKLNMZZCX201812), and “Double World-classes” Construction Program of China Pharmaceutical University (No. CPU2018GF02).

ABBREVIATIONS USED VDss, volume of distribution at steady state; CL, clearance; t1/2, terminal half-life; fu, fraction unbound in plasma; SVM, support vector machine; RF, random forest; GBM,


Page 30 of 35

Page 31 of 35 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


gradient boosting machine; XGB, XGBoost; PK, pharmacokinetics; ADME, absorption, distribution, metabolism, excretion; QSAR, quantitative structure-activity relationship; BNN, Bayesian neural networks; CART, classification and egression trees; PLS, partial least squares; MLR, multiple linear regression; ANN, artificial neural networks; BPN, back propagation neural; SMILES, Simplified Molecular Input Line Entry System; PCA, principal component analysis; MOE, Molecular Operating Environment; 2D, two-dimensional; 3D, three-dimensional; R2, squared Pearson’s correlation coefficient; RMSE, root-mean-square error; FE, fold error; AFE, average fold error; Q2, squared cross-validated correlation coefficient; LOO, leave-one-out; OECD, the Organization for Economic Cooperation and Development; AD, applicability domain; vdw, van der Waals; Y-rand, Y-randomization; LCO, leaveclass-out; LMO, leave-many-out.

REFERENCES (1) Norinder, U.; Bergström, C. A. Prediction of ADMET Properties. ChemMedChem. 2010, 37, 920937. (2) Lombardo, F.; Berellini, G.; Obach, R. S. Trend Analysis of a Database of Intravenous Pharmacokinetic Parameters in Humans for 1352 Drug Compounds. Drug metab. dispos. 2018, 46, 1466-1477. (3) Di, L.; Feng, B.; Goosen, T. C.; Lai, Y.; Steyn, S. J.; Varma, M. V.; Obach, R. S. A Perspective on the Prediction of Drug Pharmacokinetics and Disposition in Drug Research and Development. Drug metab. dispos. 2013, 41, 1975-1993. (4) Ward, K. W.; Smith, B. R. A Comprehensive Quantitative and Qualitative Evaluation of Extrapolation of Intravenous Pharmacokinetic Parameters from Rat, Dog, and Monkey to Humans. I. Clearance. Drug Metab. Disposition. 2004, 32, 603-611. (5) Toutain, P. L.; Bousquet-Mélou, A. Plasma Terminal Half-Life. J. Vet. Pharmacol. Ther. 2010, 27, 427-439. (6) Franco, L.; R Scott, O.; Shalaeva, M. Y.; Feng, G. Prediction of Volume of Distribution Values in Humans for Neutral and Basic Drugs Using Physicochemical Measurements and Plasma Protein Binding Data. J. Med. Chem. 2002, 45, 2867-2876. (7) M Paul, G.; Waters, N. J.; Paine, S. W.; Davis, A. M. In Silico Human and Rat Vss Quantitative Structure-Activity Relationship Models. J. Med. Chem. 2006, 49, 1953. (8) Sui, X.; Sun, J.; Li, H.; Wang, Y.; Liu, X.; Zhang, W.; Chen, L.; He, Z. Prediction of Volume of Distribution Values in Human Using Immobilized Artificial Membrane Partitioning Coefficients, the Fraction of Compound Ionized and Plasma Protein Binding Data. Eur. J. Med. Chem. 2009, 44, 44554460. (9) Giuliano, B.; Clayton, S.; Waters, N. J.; Franco, L. In Silico Prediction of Volume of Distribution in Human Using Linear and Nonlinear Models on a 669 Compound Data Set. J. Med. Chem. 2009, 52, 4488-4495. (10) Zhivkova, Z.; Doytchinova, I. Prediction of Steady-State Volume of Distribution of Acidic Drugs by Quantitative Structure–Pharmacokinetics Relationships. J. Pharm. Sci. 2011, 101, 1253-1266. (11) Gombar, V. K.; Hall, S. D. Quantitative Structure-Activity Relationship Models of Clinical



Pharmacokinetics: Clearance and Volume of Distribution J. Chem. Inf. Model. 2013, 53, 948-957. (12) Zhivkova, Z.; Tsvetelina, M.; Doytchinova, I. Quantitative Structure-Pharmacokinetics Relationships Analysis of Basic Drugs Volume of Distribution. J. Pharm. Pharm. Sci. 2015, 18, 515 527. (13) Del, A. E. M.; Ghemtio, L.; Xhaard, H.; Yliperttula, M.; Urtti, A.; Kidron, H. Applying Linear and Non-linear Methods for Parallel Prediction of Volume of Distribution and Fraction of Unbound Drug. PloS One. 2013, 8, e74758. (14) Schneider, G.; Coassolo, P., T Combining in Vitro and in Vivo Pharmacokinetic Data for Prediction of Hepatic Drug Clearance in Humans by Artificial Neural Networks and Multivariate Statistical Techniques. J. Med. Chem. 1999, 42, 5072. (15) Paulo, P.; Gouveia, L. F.; Morais, J. A. G. Prediction of the in Vitro Intrinsic Clearance Determined in Suspensions of Human Hepatocytes by Using Artificial Neural Networks. Eur. J. Pharm. Sci. 2010, 39, 310-321. (16) Ekins, S.; Obach, R. S. Three-Dimensional Quantitative Structure Activity Relationship Computational Approaches for Prediction of Human in Vitro Intrinsic Clearance. J. Pharmacol. Exp. Ther. 2000, 295, 463-473. (17) Wajima, T.; Fukumura, K.; Yano, Y.; Oguma, T. Prediction of Human Clearance from Animal Data and Molecular Structural Parameters Using Multivariate Regression Analysis. J. Pharm. Sci. 2002, 91, 2489-2499. (18) Li, H.; Sun, J.; Sui, X.; Liu, J.; Yan, Z.; Liu, X.; Sun, Y.; He, Z., He First-Principle, StructureBased Prediction of Hepatic Metabolic Clearance Values in Human. Eur. J. Med. Chem. 2009, 44, 1600-1606. (19) Zuegge, J.; Schneider, G.; Coassolo, P.; Lavé, T. Prediction of Hepatic Metabolic Clearance: Comparison and Assessment of Prediction Models. Clin. Pharmacokinet. 2001, 40, 553. (20) Arnot, J. A.; Brown, T. N.; Frank, W. Estimating Screening-Level Organic Chemical Half-Lives in Humans. Environ. Sci. Technol. 2013, 48, 723-30. (21) Lu, J.; Lu, D.; Zhang, X.; Bi, Y.; Cheng, K.; Zheng, M.; Luo, X. Estimation of Elimination HalfLives of Organic Chemicals in Humans using Gradient Boosting Machine. Biochim. Biophys. Acta, Gen. Subj. 2016, 1860, 2664-2671. (22) Macielag, M. J., Chemical Properties of Antimicrobials and Their Uniqueness. Antibiotic Discovery and Development: Springer, US, 2012. (23) Veber, D. F.; Johnson, S. R.; Hung-Yuan, C.; Smith, B. R.; Ward, K. W.; Kopple, K. D. Molecular Properties that Influence the Oral Bioavailability of Drug Candidates. J. Med. Chem. 2002, 45, 2615. (24) Mercader, A. G.; Duchowicz, P. R.; Fernández, F. M.; Castro, E. A. Modified and Enhanced Replacement Method for the Selection of Molecular Descriptors in QSAR and QSPR Theories. Chemometrics Intellig. Lab. Syst. 2008, 92, 138-144. (25) Vilar, S.; Cozza, G.; Moro, S. Medicinal Chemistry and the Molecular Operating Environment (MOE): Application of QSAR and Molecular Docking to Drug Discovery. Curr. Top. Med. Chem. 2008, 8, 1555-1572. (26) Su, B. H.; Shen, M.; Esposito, E. X.; Hopfinger, A. J.; Tseng, Y. J. In Silico Binary Classification QSAR Models Based on 4D-Fingerprints and MOE Descriptors for Prediction of hERG Blockage. J. Chem. Inf. Model. 2010, 50, 1304-1318. (27) Moriwaki, H.; Tian, Y. S.; Kawashita, N.; Takagi, T. Mordred: A Molecular Descriptor Calculator. J. Cheminform. 2018, 10, 4. (28) Yap, C. W. PaDEL-Descriptor: An Open Source Software to Calculate Molecular Descriptors and Fingerprints. J. Comput. Chem. 2011, 32, 1466-1474. (29) Kursa, M. B.; Jankowski, A.; Rudnicki, W. R. Boruta-A System for Feature Selection. Fund. Inform. 2010, 101, 271-285. (30) Cherkassky, V.; Ma, Y. Practical Selection of SVM Parameters and Noise Estimation for SVM Regression. Neural Netw. 2004, 17, 113-126. (31) Cutler, A.; Cutler, D. R.; Stevens, J. R. Random Forests. Mach. Learn. 2004, 45, 157-176. (32) Friedman, J. H. Greedy Function Approximation: A Gradient Boosting Machine. Ann. Stat. 2001,


Page 32 of 35

Page 33 of 35 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


29, 1189-1232. (33) Chen, T.; Guestrin, C., XGBoost: A Scalable Tree Boosting System. In Proceedings of the 22Nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, ACM: San Francisco, 2016. (34) Lombardo, F.; Jing, Y. In Silico Prediction of Volume of Distribution in Human. Extensive Data Set and the Exploration of Linear and Non-linear Methods Coupled with Molecular Interaction Fields Descriptors. J. Chem. Inf. Model. 2016, 56, 2042-2052. (35) Wen, Z.; Zhang, R.; Ramamohanarao, K.; Yang, L. Scalable and Fast SVM Regression Using Modern Hardware. World Wide Web. 2017, 21, 261-287. (36) Bredensteiner, E. J.; Bennett, K. P. Multicategory Classification by Support Vector Machines. Computational Optimization and Applications. 1999, 12, 53-79. (37) Natekin, A.; Knoll, A. Gradient Boosting Machines, a Tutorial. Front. Neurorobot. 2013, 7, 21. (38) Tedesco, M.; Pulliainen, J.; Takala, M.; Hallikainen, M.; Pampaloni, P. Artificial Neural Networkbased Techniques for the Retrieval of SWE and Snow Depth from SSM/I Data. Remote Sens. Environ. 2004, 90, 76-85. (39) Houston, J. B.; Carlile, D. J. Prediction of Hepatic Clearance from Microsomes, Hepatocytes, and Liver Slices. Drug Metab. Rev. 1997, 29, 891-922. (40) Klopman, G.; Kalos, A. N. Causality in Structure-Activity Studies. J. Comput. Chem. 2010, 6, 492–506. (41) Sahigara, F.; Mansouri, K.; Ballabio, D.; Mauri, A.; Consonni, V.; Todeschini, R. Comparison of Different Approaches to Define the Applicability Domain of QSAR Models. Molecules. 2012, 17, 4791-4810. (42) Netzeva, T. I.; Worth, A. P.; Aldenberg, T.; Benigni, R.; Cronin, M. T. D. Current Status of Methods for Defining the Applicability Domain of (Quantitative) Structure-Activity Relationships : The Report and Recommendations of ECVAM Workshop 52. Altern. Lab. Anim. 2005, 33, 155-173. (43) Kamari, A.; Bahadori, A.; Mohammadi, A. H.; Zendehboudi, S. New Tools Predict Monoethylene Glycol Injection Rate for Natural Gas Hydrate Inhibition. J. Loss Prev. Process Indust. 2015, 33, 222231. (44) Tropsha, A. Best Practices for QSAR Model Development, Validation, and Exploitation. Mol. Inform. 2010, 29, 476-488. (45) Zhu, L.; Zhao, J.; Zhang, Y.; Zhou, W.; Yin, L.; Wang, Y.; Fan, Y.; Chen, Y.; Liu, H. ADME Properties Evaluation in Drug Discovery: In Silico Prediction of Blood-Brain Partitioning. Mol. Diversity. 2018, 22, 979-990. (46) Pearlman, R. S.; Smith, K. M. Novel Software Tools for Chemical Diversity. Perspect. Drug Discov. Des. 1998, 9-11, 339-353. (47) Giuliano, B.; Waters, N. J.; Franco, L. In Silico Prediction of Total Human Plasma Clearance. J. Chem. Inf. Model. 2012, 52, 2069. (48) Lu, Z.; Wang, J.; Wientjes, M. G.; Au, J. L. Intraperitoneal Therapy for Peritoneal Cancer. Future Oncol. 2010, 6, 1625-1641. (49) Winiwarter, S.; Ax, F.; Lennernäs, H.; Hallberg, A.; Pettersson, C.; Karlén, A. Hydrogen Bonding Descriptors in the Prediction of Human in Vivo Intestinal Permeability. J. Mol. Graph. Model. 2003, 21, 273-287. (50) Abraham, M. H.; Chadha, H. S.; Martins, F.; Mitchell, R. C.; Bradbury, M. W.; Gratton, J. A. Hydrogen Bonding Part 46: A Review of the Correlation and Prediction of Transport Properties by an lfer Method: Physicochemical Properties, Brain Penetration and Skin Permeability. Pest Manage. Sci. 2015, 55, 78-88. (51) Wildman, S. A.; Crippen, G. M. Prediction of Physicochemical Parameters by Atomic Contributions. J. Chem. Inf. Comput. Sci. 1999, 39, 868-873. (52) Pickett, S. D.; Mclay, I. M.; Clark, D. E. Enhancing the Hit-to-Lead Properties of Lead Optimization Libraries. J. Chem. Inf. Comput. Sci. 2000, 40, 263-272. (53) Lu, L.; Zhan, T.; Ma, M.; Xu, C.; Wang, J.; Zhang, C.; Liu, W.; Zhuang, S. Thyroid Disruption by Bisphenol S Analogues via Thyroid Hormone Receptor β: in Vitro, in Vivo and Molecular Dynamics



Simulation Study. Environ. Sci. Technol. 2018, 52, 6617-6625. (54) Lv, X.; Pan, L.; Wang, J.; Lu, L.; Yan, W.; Zhu, Y.; Xu, Y.; Guo, M.; Zhuang, S. Effects of Triazole Fungicides on Androgenic Disruption and CYP3A4 Enzyme Activity. Environ. Pollut. 2017, 222, 504-512. (55) Stanton, D. T. Evaluation and Use of BCUT Descriptors in QSAR and QSPR Studies. J. Chem. Inf. Model. 2010, 39, 11-20.


Page 34 of 35

Page 35 of 35 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Table of Contents In Silico Prediction of Human Intravenous Pharmacokinetic Parameters with Improved Accuracy Yuchen Wang a,1, Haichun Liua,1, Yuanrong Fana, Xingye Chena, Yan Yanga, Lu Zhua, Junnan Zhaoa, Yadong Chena,* and Yanmin Zhanga,* a School

of Science, China Pharmaceutical University, 639 Longmian Avenue, Nanjing, Jiangsu 211198, China. 1 The authors contributed equally to this work and should be considered as co-first authors


In Silico Prediction of Human Intravenous Pharmacokinetic

Recommend Documents