QSAR models for prediction study of HIV protease inhibitors using support vector machines, neural networks and multiple linear regression

2.2.1 Multiple linear regression

MLR is a statistical tool that regresses independent variables against a dependent variable. The objective of MLR is to find a linear model of the property of interest, which takes the form below: $$y=a_0+{\displaystyle \sum _{i=1}^n}{a}_i{x}_i$$ where y is the property which is the dependent variable, x_i represents molecular descriptors, a_i represents the coefficients of those descriptors and α₀ is the intercept of the equation.

2.2.2 Artificial neural networks

ANNs are artificial systems simulating the function of the human brain. Three components constitute a neural network: the processing elements or nodes, the topology of the connections between the nodes, and the learning rule by which new information is encoded in the network. While there are a number of different ANN models, the most frequently used type of ANN in QSAR is the three-layered feed-forward network (Zupan and Gasteiger, 1993). In this type of networks, the neurons are arranged in layers (an input layer, one hidden layer and an output layer). Each neuron in any layer is fully connected with the neurons of a succeeding layer and no connections are between neurons belonging to the same layer.

According to the supervised learning adopted, the networks are taught by giving them examples of input patterns and the corresponding target outputs. Through an iterative process, the connection weights are modified until the network gives the desired results for the training set of data. A back-propagation algorithm is used to minimize the error function. This algorithm has been described previously with a simple example of application (Cherqaoui and Villemin, 1994) and a detail of this algorithm is given elsewhere (Freeman and Skapura, 1991).

2.2.3 Support vector machine

SVM is gaining popularity due to many attractive features and promising empirical performance. It originated from early concepts developed by Cortes and Vapnik (1995). This method has proven to be very effective for addressing general purpose classification and regression problems. The main advantage of SVM is that it adopts the structure risk minimization (SRM) principle, which has been shown to be superior to the traditional empirical risk minimization (ERM) principle (Burges, 1998), employed by conventional neural networks. SRM minimizes an upper bound of the generalization error on Vapnik–Chernoverkis dimension (“generalization error”), as opposed to ERM that minimizes the training error. So SVM is usually less vulnerable to the overfitting problem. Since various introductions into SVM were already stated before (Cristianini and Shawe-Taylor, 2000), only the main ideas about SVM are given in this paper.

SVM can be applied to regression problems by the introduction of an alternative loss function that is modified to include a distance measure. Considering the problem of approximating the set of data $G={\{(x_i\text{,}{d}_i)\}}_{i=1}^n$ (x_i is the input vector, d_i is the desired value, and n is the total number of data patterns). In SVM method, the regression function is approximated, in a feature space F, by the following function:

In Eq. (2), the first term $\frac{1}{2}{w}^{T}w=\frac{1}{2}\Vert w{\Vert }^2$ is called regularized term. Minimizing $\frac{1}{2}\Vert w{\Vert }^2$ will make a function as flat as possible, thus playing role of controlling the function capacity. The second term is the empirical error measured by the ε-insensitive loss function, which is defined by Eq. (3). This defines a ε tube so that if predicted value is within the tube, the loss is zero, while if predicted point is outside the tube, the loss is the magnitude of the difference between the predicted value and the radius ε of the tube. C is penalty parameter, which is a regularized constant to determine the trade-off between training error and model flatness. To get the estimations of w and b, Eq. (2) is transformed into the primal objective Eq. (4) by introducing ξ_i and ${\xi }_i^{\ast }$ (slack variables representing upper and lower constraints on the outputs of the system).

In Eq. (6), α_i and ${\alpha }_i^{\ast }$ are the introduced Lagrange multipliers. They satisfy the equality ${\alpha }_i{\alpha }_i^{\ast }=0\text{,}$ α_i ⩾ 0, ${\alpha }_i^{\ast }⩾0$ (i = 1, …, n) and are obtained by maximizing the dual form of Eq. (4) which has the following form:

Through selecting the appropriate kernel function, the non-linear relation between the building cooling load and its correlative influence parameters based on SVM is established.

Any function satisfying Mercer’s (1909) condition can be used as the kernel function, and the typical kernel functions include linear, polynomial, Gaussian and sigmoid functions. Among these functions, the Gaussian function can map the sample set from the input space into a high dimensional feature space effectively, which is good for representing the complex non-linear relationship between the output and input samples. Moreover, there is only one variable (the width parameter) in it needed to be determined, which ensures the high calculation efficiency. Because of the above advantages, the Gaussian function is used widely. In this paper, Gaussian function is also selected as the kernel function, whose expression is shown as follows:

All SVM models, in our present study, were implemented using the software LIBSVM for classification and regression developed by Chang and Lin. All calculation programs implementing ANN were written in an M-file based on the MATLAB script, developed in our laboratory.

2.3.1 Root Mean Square Error (RMSE)

The residual between observed and estimated data is evaluated. This index assumes that larger estimated errors are of greater importance than smaller ones; hence they are given a more than proportionate penalty. The RMSE is known to be descriptive when the prediction capability among predictors is compared, it is defined by: $$\text{RMSE}=\sqrt{\frac{1}{N}{\displaystyle \sum _{m=1}^N}{(y_m-{\widehat{y}}_m)}^2}$$

2.3.2 Scatter Index (SI)

SI is a standard metric for wave model intercomparison. Essentially, it is a normalized measure of error that takes into account the observed data. It is defined as: $$\text{SI}=\frac{\text{RMSE}}{{\overline{y}}_m}$$

Lower values of the SI are an indication of a better prediction.

2.3.3 Correlation coefficient

The predictive power of the QSAR models developed on the selected training sets is estimated on the predictions of an external test set and also examined by cross validation test to the calculating statistical parameters of Q²: $$Q^2=1-\frac{{\sum }_{m=1}^N{({\widehat{y}}_m-y_m)}^2}{{\sum }_{m=1}^N{(y_m-\overline{y})}^2}\text{.}$$

2.3.4 Average Absolute Relative Error (AARE)

AARE indicates the relative absolute deviation in percent from the experimental values. It is defined as: $$\text{AARE}(\%)=\frac{1}{N}{\displaystyle \sum _{m=1}^{N_s}}\left|\frac{y_m-{\widehat{y}}_m}{y_m}\right|\times 100$$

A lower value implies a better correlation.

In these equations, y_m is the desired output, ${\widehat{y}}_m$ is the predicted value by model, $\overline{y}$ is the mean of dependent variable, and N is the number of the molecules in data set.

3.1.1 Support vector machine

The training of the SVM model included the selection of capacity parameter C, ε of ε-insensitive loss function and the corresponding parameters of the kernel function. Firstly, the kernel function should be decided, which determines the sample distribution in the mapping space. Generally, using RBF kernel function will yield better prediction performance (Nianyi et al., 2004), and it was used as the SVM model's kernel in this study accordingly. The radial basis function used is: $$\exp (-\gamma \Vert \mu -\nu {\Vert }^2)$$ where γ is the parameter of the kernel, μ and ν are two independent variables.

Secondly, corresponding parameters, i.e., γ of the kernel function greatly affects the number of support vectors, which has a close relation with the performance of the SVM and training time. Over many support vectors could produce overfitting and increase the training time. In addition, γ controls the amplitude of the RBF function, and therefore, controls the generalization ability of SVM.

Parameter ε-insensitive prevents the entire training set meeting boundary conditions and so allows for the possibility of sparsity in the dual formulation's solution. The optimal value for ε depends on the type of noise present in the data, which is usually unknown.

Lastly, the effect of capacity parameter C was tested. It controls the trade-off between maximizing the margin and minimizing the training error. If C is too small then insufficient stress will be placed on fitting the training data. If C is too large then the algorithm will overfit the training data. However, Wang et al. (2005) indicated that prediction error was scarcely influenced by C. To make the learning process stable, a large value should be set up for C.

The grid optimization of LIBSVM was used to find the optimal values of the C, γ and ε parameters when using the radial basis function kernel in the SVM mode. In the grid for data set, a series of C values ranging from 50 to 600 with incremental steps of 50, γ in the range from 0.005 to 0.05 with incremental steps of 0.005 and ε from 0.05 to 0.5 with incremental steps of 0.05 have been exploited. The optimal values of C, γ and ε are identified to be 150, 0.02 and 0.5, respectively. The values of R² and RMSE are 0.94 and 0.275, respectively. Because the grid search was performed over three parameters, the RMSE could not be shown in one plot to be visualized easily. Figs. 2a–2c show the influence of each parameter with the other two fixed to the optimal values on the model performance.

Close

Figure 2a

RMSE versus C (γ = 0.01, ε = 0.25).

Export to PPT

Close

Figure 2b

RMSE versus γ (C = 100, ε = 0.25).

Export to PPT

Close

Figure 2c

RMSE versus ε (C = 100, γ = 0.01).

Export to PPT

3.1.2 Artificial neural networks

All the feed-forward ANNs used in this paper are three-layer networks with 10 units (10 molecular descriptors) in the input layer, a variable number of hidden neurons, and one unit (log(1/K_i)) in the output layer. A bias term was added to the input and hidden layers. Each neuron in any layer is fully connected with the neurons of a succeeding layer. There are neither connections between the neurons within a layer nor any direct connection between those of the input and the output layers. Input and output data are normalized between 0.1 and 0.9. The sigmoid function was used as the transformation function (Freeman and Skapura, 1991). The weights of the connections between the neurons were initially assigned with random values uniformly distributed between −0.5 and +0.5 and no momentum was added. The learning rate was initially set to 1 and was gradually decreased during training. The back-propagation algorithm was used to adjust those weights.

One major problem in neural network is how to determine the number of nodes in the hidden layer. Though there is no rigorous rule to rely on (Golbraikh and Tropsha, 2002), a practical way is to use a ratio, ρ, to determine the number of hidden units (Andrea and Kalayeh, 1991). ρ is defined as the following: $$\rho =\frac{\text{number of compounds presented to the network}}{\text{number of connections in the network}}\text{.}$$ The reasonable value of ρ should be between 1.0 and 3.0. If ρ < 1, the network simply memorizes the data, whereas if ρ > 3, the network is not able to generalize (Golbraikh and Tropsha, 2002). In our study, the four selected variables were used as input and the analgesic activity was used as output, so a 4-x-1 (x represents the number of hidden neurons) network was constructed and the suitable value of x could be defined from 2 to 5 (ρ was located between 1.0 and 2.5). Each architecture was trained with 10 different initial random sets of weights and with the number of cycles limited to 1000. In all cases 100 cycles were enough to obtain stable results. The results are reported in Table 2. Among all architectures of ANN, the best one is 4-4-1 (R² = 0.92 and RMSE = 0.319).

3.1.3 Multiple linear regression

The linear function constructed from the four molecular descriptors and the training set has the following form: $$\log \left(\frac{1}{K_i}\right)=5.82+0.98\times \Delta -0.08\times \mathit{MgVol}+2.79\times I_a-0.96\times I_o$$ $$N=30\text{,}{R}^2=0.81\text{,}\text{RMSE}=0.535\text{.}$$

3.2.1 Internal validation

The internal validation technique used is cross-validation (CV). CV is a popular technique used to explore the reliability of statistical models. Based on this technique, a number of modified data sets are created by deleting in each case one or a small group (leave-some-out) of objects. For each data set, an input–output model is developed, based on the utilized modeling technique. The model is evaluated by measuring its accuracy in predicting the responses of the remaining data (the ones that have not been utilized in the development of the model). The leave-one-out (LOO) procedure was utilized, in this study, which produces a number of models by deleting one from the whole data set (38 cyclic-urea derivatives).

Five ANN architectures of 4-x-1 (x = 1–5) have been tested. The results of QSAR done by these ANN architectures, by MLR analysis and by SVM method are listed in Table 3. The quality of the fitting is estimated by the RMSE and by the statistical parameter Q. As it can be seen in Table 3, high correlation coefficient (Q² = 0.89) and low RMSE = 0.171, SI = 0.022 have been obtained by means of the SVM. According to this table, it is clear that the performance of SVM is better than those obtained by ANN and MLR techniques. Indeed, in every case, the SVM's correlation coefficient is greater and its standard deviation is lower than those of the ANN and MLR.

3.2.2 External validation

In order to estimate the predictive power of SVM, MLR and ANN, we must use a set of compounds which have not been used for training set (used for establishing the QSAR model). The models established in the computation procedure, by using the 30 cyclic-urea derivatives, are used to predict the activity of the remaining eight compounds. The plot of predicted versus experimental values for data set is shown in Fig. 3a (SVM), Fig. 3b (ANN) and Fig. 3c (MLR). Among all these figures, the first one shows that the activity values calculated by the SVM are very close to the experimental ones. The statistical parameters of the three models are shown in Table 4. As can be seen from this table, the statistical parameters of SVM model are better than the other ones.

Close

Figure 3a

log(1/K_i) observed experimentally versus log(1/K_i) predicted by SVM.

Export to PPT

Close

Figure 3b

log(1/K_i) observed experimentally versus log(1/K_i) predicted by ANN.

Export to PPT

Close

Figure 3c

log(1/K_i) observed experimentally versus log(1/K_i) predicted by MLR.

Export to PPT

3.2.3 Y-randomization test

Y-randomization is an attempt to observe the action of chance in fitting given data. In other words it is applied to exclude the possibility of chance correlation. This technique ensures the robustness of a QSAR model (Tropsha et al., 2003; Tropsha and Golbraikh, 2007). The dependent variable vector [y = log(1/K_i)] is randomly shuffled and a new QSAR model is developed using the original molecular descriptors. The new QSAR models (after several repetitions) are expected to have low R² values. If the opposite happens then an acceptable QSAR model cannot be obtained for the specific modeling method and data.

In this work, 10 random shuffles of the y vector were performed for SVM, ANN and MLR. The results are shown in Table 5. For each technique, the mean value of random models is significantly lower than the corresponding value of the non-random model. This suggests that the models are not obtained by chance.

3.2.4 Domain of applicability

The domain of application (Eriksson et al., 2003) of a QSAR model must be defined if the model is to be used for screening new compounds. Predictions for only those compounds that fall into this domain may be considered reliable. Extent of Extrapolation (Gramatica, 2007) is one simple approach to define the applicability of the domain. It is based on the calculation of the leverage h_i for each chemical, for which QSAR model is used to predict its activity: $$h_i=x_i^T{(X^{T}X)}^{-1}{x}_i\text{,}i=1\text{,}\dots \text{,}n$$ where x_i is the descriptor vector of the considered compound and X is the descriptor matrix derived from the training set descriptor values. The superscript T refers to the transpose of the matrix/vector. The warning leverage h^∗ is, generally, fixed at 3(k + 1)/N, where N is the number of training compounds and k is the number of model parameters. A leverage greater than the warning leverage h^∗ means that the predicted response is the result of substantial extrapolation of the model and, therefore, may not be reliable.

The Williams plot for the presented SVM model is shown in Fig. 4. From this plot, the applicability domain is established inside a squared area within (±3s) standard deviations and a leverage threshold h^∗ of 0.5. As shown in the Williams plot (Fig. 4), h_i values of all the compounds in the training and test sets are lower than the warning value (h^∗ = 0.5). None of the compounds are particularly influential in the model space and the training set has great representativeness. For all the compounds in the training and test sets, their standardized residuals are smaller than three standard deviation units (2s). This means that all predicted values are acceptable.

Close

Figure 4

Williams plot of the current QSAR model.

Export to PPT