An In Silico Method for Screening Nicotine Derivatives as Cytochrome P450 2A6 Selective Inhibitors Based on Kernel Partial Least Squares

Nicotine and a variety of other drugs and toxins are metabolized by cytochrome P450 (CYP) 2A6. The aim of the present study was to build a quantitative structure-activity relationship (QSAR) model to predict the activities of nicotine analogues on CYP2A6. Kernel partial least squares (K-PLS) regression was employed with the electro-topological descriptors to build the computational models. Both the internal and external predictabilities of the models were evaluated with test sets to ensure their validity and reliability. As a comparison to K-PLS, a standard PLS algorithm was also applied on the same training and test sets. Our results show that the K-PLS produced reasonable results that outperformed the PLS model on the datasets. The obtained K-PLS model will be helpful for the design of novel nicotine-like selective CYP2A6 inhibitors.


Introduction
Cytochrome P450 2A6 (CYP2A6), the major coumarin 7-hydroxylase present in human liver (Cashman, etc., 1992;Pearce, etc., 1992;Shimada, etc., 1996), is known to metabolize a variety of compounds including quinoline (Reigh, etc., 1996), nicotine (Nakajima, etc., 1996), cotinine (Nakajima, etc., 1996), and various N-nitroso compounds present in cigarette smoke (Guengerich, etc., 1994). Hepatic CYP2A6 catalyses the major route of nicotine metabolism via the intermediacy of the aldehyde oxidase-catalyzed iminium ion that is converted to the metabolite, cotinine. (Cashman, etc., 1992;Tricker, 2003;Hukkanen, etc., 2005). The efficiency of CYP2A6-mediated metabolism of nicotine is closely related to the specific concentration of nicotine in blood for keeping addiction liability. Potent and specific inhibitors of the CYP2A6 enzyme might improve nicotine bioavailability and thus make oral nicotine administration feasible in smoking cessation therapy. The inhibition of CYP2A6 may decrease the number of cigarettes a person needs to smoke to obtain their desired blood nicotine concentration. Nowadays, a number of compounds tested as CYP2A6 inhibitors possess strong inhibitory effects (Draper, etc., 1997;Maenpaa, etc., 1993;Fujita, etc., 2003). However, to our knowledge, no compounds have been characterized as both potent and selective CYP2A6 inhibitors. In the present study QSAR models were established based on a series of nicotine derivatives, with the ultimate aim of aiding the prediction and development of a potent and specific CYP2A6 inhibitor. The in silico methods were built employing electrotopological state descriptors by using kernel partial least squares (K-PLS), a relatively novel method in chemometrics compared to the partial least squares (PLS) method.
The partial least squares method (Wold, 1975;Wold, etc., 1984) has been a popular modeling, regression, discrimination and classification technique in its domain of origin chemometrics. In its general form PLS creates orthogonal score vectors by using the existing correlations between different sets of variables while also keeping most of the variance of all sets. It is a statistical tool specifically designed to deal with multiple regression problems, where the number of observations is limited, the missing data are numerous and the correlations between the predictor variables are high.
PLS has proven to be useful in situations where the number of observed variables is much greater than the number of observations and high multicollinearity among the variables exists. This situation is quite common in the case of kernel-based learning where the original data are mapped to a highdimensional feature space corresponding to a reproducing kernel Hilbert space. Too high dimensions also cause problems like overfitting, thus leading to the decrease of the prediction accuracy of the external data. As an alternative to PLS, a nonlinear PLS has been newly developed based on kernel methods, i.e., kernel partial least squares. In the next section, a detailed description of K-PLS was offered.
The outline of the paper is as follows. The kernel partial least squares analysis was introduced based on an optimization-derived method. QSAR models were built for nicotine analogues employing K-PLS for a library of 58 nicotine analogues as CYP2A6 selective inhibitors (Denton, etc., 2005). Finally, PLS and K-PLS were compared to determine which exhibits superior performance.

Kernel partial least squares
As a generic kernel regression method, kernel partial least squares has been proven to be more competitive, and even more stable than other kernel regression algorithms such as support vector machines (SVM) and kernel ridge regression, and this method is also much more easily implemented (John and Nello, 2004).
The idea of the kernel PLS is developed based on the mapping of the original Ξ-space data into a high-dimensional feature space. A kernel is a continuous function κ: Ξ × Ξ → Ρ for which there exists an Φ inner product space as a representation space and a map φ : Ξ → Φ such that for all x, y∈ Ξ This definition allows us to perform calculations in the Φ space in an implicit way, by substituting the scalar product operation with its corresponding kernel version.
In the following part, a derivation of Direct Kernel Partial Least Squares (DK-PLS) based on the optimization algorithm (Bennett and Embrechts, 2003) for nonlinear regression is introduced. The DK-PLS is developed on the basis of a direct factorization of the kernel matrix. DK-PLS has the advantage that the kernel does not need to be square, which factorizes the kernel matrix directly and then the final regression function is computed based on this factorization. We provide here the simplified algorithm for one response variable, which is more popular in QSAR modeling.
Lets consider the data sample (X, Y) where ; X and Y represent the variable matrix and the response matrix (normally a one-dimensional vector), respectively. First to define a Gram matrix in feature space: 0 ( ) ( ') Let K c be the centered form of K 0 , the Y' = y has been normalized to have mean 0 and standard deviation 1. Let M be the desired number of latent variables. It should be noted that the test data should be centralized before, according to the following formula: where 1 is the vector of element 1, I is the unit matrix. As we can see that this algorithm is easy to be complemented using C-or other languages. This derivation should make the PLS algorithm more accessible to machine learning researchers and popularly used for chemometrics applications.
Meanwhile, in order to compare the performances of K-PLS and PLS methods on the data set, the Partial Least Squares regression using the SIMPLS algorithm is also proposed (Jong, 1993). The same training and test sets are applied for both K-PLS and PLS models.

Data set
In the present study, we used a data set of 55 nicotine analogues whose selective inhibition on CYP2A6 was reported in the literature (Denton, etc., 2005). All these compounds were shown in Tables 1 and 2. The relative potency of the analogues, expressed by K i values, on the functional activity of cDNA-expressed human CYP2A6 were determined by examining coumarin 7-hydroxylation (Denton, etc., 2005). Several molecules (Tables 1, 2) with un-deterministic chemical structure such as molecule 38b in the original paper (Denton, etc., 2005) were omitted in this work. In order to guarantee the linear distribution of the biological data, the K i values were transformed into -LogK i .   The molecular descriptors in Table 3 were computed using MolconnZ program. The model has been trained by a training set (3/4 of the whole data) and validated by an independent test set (1/4) ( Table 1).

Molecular descriptors
In the past decade, electrotopological state (E-state) indices have been used for correlating a variety of physicochemical and biological properties of chemical compounds. The E-state indices are computed for each atom in a molecule and encode information about both the topological environment of that atom and the electronic interactions due to all other atoms in the molecule (Kier and Hall, 1990). E-state indices have been found to be very useful in building QSAR models (Wang, etc., 2004;Wang, etc., 2005a;Wang, etc., 2005b). In this work, the E-state descriptors with detailed definitions are indicated in Table 3. For the present data, the sum of intrinsic state (sumI) and the sum of delta-I values (sum-delI), the extreme atom level E-State values (Hmax, Gmax and Hmin), as well as the number of hydrogen bond (H-bond) donor and acceptor are found to be useful in construction of a reliable K-PLS model. The intrinsic state encodes the valence state electronegativity of the atom as well as its local topology, which are particularly useful for describing the chemical features of a series of compounds. It has been found that the H-bonding and aromatic-aromatic interactions are very important for the binding of nicotine analogues with CYP2A6 enzyme (Yano, etc, 2006). Our present model, similarly, demonstrates that two descriptors of nHBa and nwHBa describing the hydrogen bonding interaction, as well as another descriptor of SwHBa describing aromatic carbons, are crucial for the functioning of nicotine inhibitors. All these descriptors possibly revealed that the nicotine derivatives play a main role of hydrogen-bond donor when interacting with the P450 enzyme. The importance of these descriptors also proved the previous results obtained by Yano, etc. (Yano, etc, 2006).

K-PLS parameters
K-PLS performs as well as or better than support vector regression for moderately-sized problems with the advantages of simple implementation, less training cost, and easier tuning of parameters. The most critical and demanding phase of any K-PLS model is the definition of kernels and the determination of parameters.
From the functions available, three types of kernels are popularly used in both SVM and K-PLS, i.e., linear, polynomial (a quadratic kernel function is normally applied) and radial basis function (Gaussian kernel), or to obtain complex kernels by combining simpler ones. The Gaussian kernel is possibly the simplest and effective kernel functions used in many cases. Therefore, in this work in the case of the kernel transformations we used a Gaussian RBF kernel function, which has the form: Before generating the kernel, all the data have been firstly Mahalanobis scaled to have mean 0 and standard deviation 1. The value of w (width) for the Gaussian kernel should be tuned before the calculations proceed. In this work, the w values varying from 1 to 8 are assigned for the Gaussian function. The number of components was randomly assigned as 3, as this value did not influence the optimal choice of w values. Correlation coefficients (R) of predicted versus measured -logK i s, as well as the mean squared errors (MSE) were determined for each method to reflect their bias and precision, respectively. Fig. 1 illustrates that the MSE and the R vary with the w value for the training and test data. As can be seen from Fig. 1, with the increase of the w value, the regression errors and coefficients of both data sets come approaching to each other with small fluctuations. Although their MSEs are not identical there is no real difference in their performance. These experiments illustrated that K-PLS was less sensitive to the tuning procedure. From this figure, we can find that the K-PLS model performs best for the present case when w =4.5, with the coefficients of 0.95 and 0.70, and errors of 0.07 and 0.63 for the training and test sets, respectively (Table 4). A second aspect for application of K-PLS regression analysis is the optimal choice of the number of latent components (N). The optimal parameters could result in a better K-PLS performance. Fig. 2 illustrates what happens if a different choice of the number of the latent variables of K-PLS is made. When the number of the latent variables ranges from 1 to 6, reliable number is detected, i.e., N=3, which is reasonable for both training and test sets. From Fig. 2, one can find that the correlation coefficients for the training sets increase with the increase of the number of latent variables, which result in the decrease of regression errors. However, for test sets, R keeps almost constant, whereas resulting in a continuous increase of MSE.

Interpretation of the K-PLS model
The structure of the optimum K-PLS achieving the highest R coefficient was determined. Meanwhile, a leave-one out cross-validated Q 2 (0.41) was also obtained for the model. Fig. 3 shows the performance of this model. As can be seen from this figure, all compounds of training and test sets are equally distributed around the diagonal line y = x. The results indicate that the proposed K-PLS based model can be used in virtual screening or optimization of nicotine-like lead compounds for the inhibition of CYP2A6. From this figure, we can find that the most potent compounds like S29, S30 and S37 in the training set, or like S10 and S44 in the test set are correctly modeled. However, we also find that the prediction errors of the model for compounds S50 and S51 are big. One major reason is that the two compounds are the ones with the weakest inhibitory effects on CYP2A6. Thus, the chemical space of the model might not be big enough to cover these two compounds, although in the training sets several compounds with the same biggest K i values (S2 and S13) were deliberately included. However, even for a series of synthesized compounds, it is possible that they are sparsely distributed through the chemical space, thus making the model resulted from the study of these compounds inapplicable to other molecules (Sun, 2006). Being renovated by addition of new data in the future, the model may expand its coverage to a new applicability. Another possible reason is that those compounds with the same biggest K i values are structurally different. It is just those molecules with different structures but same activities in one data set that might cause difficulties for the derived-model to correctly predict the activity using structure-based method. In addition, the two compounds possess negative charges at physiologic pH, which may also cause the prediction in-capablility of the model, since the descriptors applied in the present model do not work with negatively charged compounds.
Based on the obtained model, we have attempted the prediction of lots of new virtual compounds for their binding abilities. Two compounds (P1, P2) with their structures shown in Table 2 were obtained with relatively potent binding affinities with CYP2A6, and their predicted pK i values are -1.35 and -0.80 respectively. The prediction attempt might be useful for advancing our work for synthetic studies of this series of compounds.

Comparisons between K-PLS and PLS
A kernel version of PLS has some important advantages, such as the ability to find non-linear, global solutions and to work with high dimensional input vectors. Different from the PLS involving two orders of correlation for the latent components, K-PLS has three or more orders of correlation for the nonlinear components. As a relatively new method K-PLS has not gained the popularity as PLS in the field of chemometrics and other relevant fields. For a comparison of performance of both PLS and K-PLS, PLS approach was also applied to build QSAR models using the same training and test tests in the present work. The number of latent components was assigned 4 based on the optimum R and MSE obtained for both training and test sets (data not shown). Finally, the structure of the optimum PLS achieving the highest R coefficient was determined. Upon inspecting the results the first thing one notices is that the nonlinear K-PLS outperforms its linear conversion. Fig. 4 depicts the optimum PLS modeling results and all of the statistical results were shown in Table 4. PLS has been widely used in the modeling of biochemical databases, but the technique is often unsuitable for predicting very complex phenomena such as the ADME/T properties of drugs. Basically, partial least squares regression is an extension of the multiple linear regression method. However, the present case is quite complex, where many compounds are structurally different, but with identical activities, such as the K i values for S2, S13 and S51 are all 67, for S5 and S40 are both 1.4, for S10, S24 and S48 are all 0.25, and for S42 and S47 are both 0.17. This fact indicates that the relationship between the structure and activity of nicotine analogues may be nonlinear. And a linear technique is usually inapplicable for the study of data sets with nonlinear relationships. This might be the reason why the K-PLS model is successful but PLS fails to produce reasonable results on the data sets. For these data sets, and even for the training set, PLS model performs badly. Based on the results depicted in Figures 3 and 4, one might conclude that K-PLS is a preferable method to PLS on these datasets, where K-PLS exhibits obvious advantages over PLS.

Conclusion
The main goal of this paper was to build a QSAR model for nicotine derivatives as selective CYP2A6 inhibitors. Another goal was also to compare the performances of kernel partial least squares and partial least squares analysis methods when being applied to QSAR modeling. Due to the nonlinearity of the data, K-PLS outperforms PLS in the present work. The above successful application of K-PLS method on nicotine derivatives will be helpful for quantitative design of nicotine analogues as selective CYP2A6 inhibitors.
This work also proposes a derivation of K-PLS based on optimization algorithms, which makes the K-PLS approach more easily applied for chemometrics field, and also more accessible to machine learning researchers. All these will promote the kernel partial least squares algorithm, a relatively novel method, to gain popularities in chemometrics applications and other fields.