Quantitative structure-activity relationship study of P2X7 receptor inhibitors using combination of principal component analysis and artificial intelligence methods.

P2X7 antagonist activity for a set of 49 molecules of the P2X7 receptor antagonists, derivatives of purine, was modeled with the aid of chemometric and artificial intelligence techniques. The activity of these compounds was estimated by means of combination of principal component analysis (PCA), as a well-known data reduction method, genetic algorithm (GA), as a variable selection technique, and artificial neural network (ANN), as a non-linear modeling method. First, a linear regression, combined with PCA, (principal component regression) was operated to model the structure-activity relationships, and afterwards a combination of PCA and ANN algorithm was employed to accurately predict the biological activity of the P2X7 antagonist. PCA preserves as much of the information as possible contained in the original data set. Seven most important PC's to the studied activity were selected as the inputs of ANN box by an efficient variable selection method, GA. The best computational neural network model was a fully-connected, feed-forward model with 7-7-1 architecture. The developed ANN model was fully evaluated by different validation techniques, including internal and external validation, and chemical applicability domain. All validations showed that the constructed quantitative structure-activity relationship model suggested is robust and satisfactory.

INTRODUCTION

The P2X7 receptor is extremely expressed by cells of hemopoietic origin. This receptor is an ATP-gated ion channel belonging to the family of ATP-sensitive ionotropic P2X receptors (12). The P2X7 receptor has a 595-amino acid sequence within average of 40% homology to other members of purinergic P2X receptors, while contains a carboxyl-terminal which is considerably longer than that of other P2X receptors (1).

The P2X7 protein possesses features such as low affinity for adenosine 5’-triphosphate (ATP) differentiating it from other P2X receptors (1). These features possibly reflect the preference of this protein for ATP as its endogenous ligand. Neumerous possible biological activities including maturation and release of interleukin-1 beta (IL-1β) and ATP-induced apoptosis for P2X7 receptor have been proposed (13). In brief, studies have demonstrated the potential role of P2X7 receptor in modulating IL-1β and possibly glutamate, to decrease nociception in neuropathic pain models (45).

In modern medicinal chemistry, it becomes to a greater extent essential to handle huge sets of structural data (6). The analysis of large data sets can be interesting as may consist of over a thousand descriptors calculated from different molecules each sampled from thousands molecular geometries. Even when the raw descriptors are converted into a data matrix, it could be dealt with more than a hundred descriptor vectors for each molecule. Such a large number of variables lead to multicollinearity, and to redundancies among the descriptors (6). As a result, it becomes more complicated to disclose patterns present in the original data. Advanced data mining methods dealing with such difficulties become tremendously more important (67). In this study, techniques are used that permit us to better understand the structure of large sets of structural data.

Data mining can be defined as the procedure of extracting helpful information from large data sets (8). Until now, a number of data mining approaches have been developed, but often a single data mining method is insufficient and, instead, more than a few methods must be used to support a single application (8). However, using different procedures to large databases causes a computational problem. A simple solution would be to reduce the amount of data by taking a subset of representative molecules from a given data set (8). Alternatively, a data compression method such as principal component analysis (PCA) can be used.

PCA has been extensively used in data mining to study data structure (6). In PCA, new orthogonal variables (latent variables or PC's) are calculated by maximizing variances of the data (6). The number of the latent variables (factors) is much lower than the number of original descriptors, so that the data can be visualized in a low-dimensional PC defined space (6910). While PCA really reduces the dimensionality of the space, it does not reduce the number of the original descriptors (the independent variables in a typical quantitative structure–activity relationship (QSAR) study), as it uses all the original descriptors to produce the new latent variables (principal components) (6910). For interpretation purposes and future investigations or model building, it would often be very useful to reduce the number of variables. PC selection can be attained either by choosing informative PC's (PC's with maximum variance) or using stochastic methods such as genetic algorithm. Several approaches exist and most of them carry out feature reduction using stepwise forward and/or backward techniques (6910). Jolliffe (11) compared a number of methods, mainly working on preserving most of the variation of the data. McCabe (12) developed techniques to remain as much information as possible by optimizing four numerical criteria (6). Rannar and coworkers (13) chose variables that span the original space as well as possible by a combination of PCA and partial least squares. In data mining, it is of importance to select a small subset of variables that can reproduce as closely as possible the structure of the complete data (6). Krzanowski (14) developed such a method based on Procrustes analysis. As the technique seeks variables by a stepwise procedure (backward elimination), there is no assurance to find the best global subset. Moreover, with hundreds or thousands of independent variables, as is often the case in data mining, intensive calculation is needed to perform PCA in each elimination step (6). In this study, a method is presented that uses a genetic algorithm (GA) to search for the best subset instead of a classical variable selection such as backward elimination procedure (6).

QSAR models can be generated employing a number of methods, including a variety of statistical methods (e.g., principal component regression (PCR)). For predicting biological activity, PCR has emerged as the statistical method of choice (1516). Artificial neural network (ANN) as a representative artificial intelligence method stands for a non-linear technique that has emerged as a potential alternative to linear regression techniques such as PCA (6171819).

ANN are not constrained by a known mathematical equation between dependent and independent variables, and have the power to model any arbitrarily complicated nonlinear relationship (16). Developers of ANN QSAR models do not require formal training in statistical methodology, and models can be generated by users with a minimum of theoretical and mathematical knowledge (16). There are a large number of researches suggesting that ANN models may offer significantly better predictive performance than traditional statistical approaches such as multiple linear regression (MLR) for certain problems such as QSAR (9102021). Thus, ANNs may represent an attractive alternative to PCA as a statistical modeling technique under certain circumstances.

In the present study both ANN and PCR techniques were used for modeling of the observed P2X7 antagonist activity of 49 studied molecules. The capability of the developed QSAR model was assessed by means of the prediction of P2X7 receptor antagonist activity of test set for which biological activity data have been reported.

In this study, we applied semi-empirical method to derive structural-chemical descriptors for the QSAR study of the 49 P2X7 antagonist activities of purine analogues. First, a linear regression, combined with PCA was operated to model the structure–activity relationships, and after that a combination of PCA and ANN algorithm was employed to predict the biological activity of the P2X7 antagonist accurately.

In this study, it is shown that ANN was superior to linear PCR in providing a good prediction of P2X7 antagonist activity of purine analogues.

MATERIALS AND METHODS

Data sets

A data set including 49 P2X7 receptor inhibitors was collected from literature (2223) (Table 1). The majority of the tested inhibitors are efficient P2X7 receptor inhibiting agents showing pIC50 values from 5.559 to 7.619. The 2D structure of each molecule was drawn by Chem Draw version 8.0 (2004) and then converted into 3D structure by Chem3D. The resulted structures were optimized using parametric method. The generated 3D structure of each molecule was visually examined to guarantee that the chirality of the chiral agent is correctly prepared and structure of molecules was not duplicated (24). Molecules were further separated into the training (32 compound), validation (7 compound) and test (10 compounds) sets based on Kennard and Stone algorithm (2526).

Table 1

Main structure and details of the compounds used in this study.

The best situation of this stage of model building is dividing data set to guarantee that both training and test sets individually cover the total space occupied by original data set. Then ideal splitting of the data set as each of objects in test set is close to at least one of the objects in the training set. Various methods were used as tools for splitting the whole original data set to the training and test sets.

The Kennard–Stone algorithm (25) selects a set of molecules in studied set of data, which are ‘uniformly’ distributed over the space defined by the candidates.

This is a classic technique to extract a representative set of molecules from a given data set. In this technique the molecules are selected consecutively. The first two objects are chosen by selecting the two farthest apart from each other. The third sample chosen is the one farthest from the first two objects, etc. Supposing that m objects have already been selected (m

Where, n stands for the number of samples in the training set, d, j=1,…,m are the squared euclidean distances from a candidate sample r, not yet included in the representative set, to the m samples already included in the representative set. One more benefit of the Kennard–Stone method is that it may be used to any matrix of predictors; there are no restrictions regarding the matrix multicollinearity. The other advantage is that the test molecules all fall inside the measured region and the training set molecules map the measured region of the input variable space completely with respect to the induced metric.

Overfitting problem or poor generalization capability occurs when a typical ANN model overlearns during the training phase. A too well-trained model may not carry out well property prediction on unseen data set due to its lack of generalization capability. A technique to solve the problem is the early stopping method in which the training process is concluded as soon as the overtraining signal appears. This technique needs the data set to be divided into three subsets including training set, test set, and validation sets. The training and the validation sets are the norm in all model formation procedures. The test set is employed to test the trend of the prediction accuracy of the model trained at some point of the training process. At later training steps, the validation error increases.

This is the moment when the model should stop to be trained to solve the overfitting problem. To achieve this task, the extracted PCs matrix was classified as discussed above. Then, the training and validation sets were employed to optimize the ANN performance.

Molecular descriptors

Structural descriptors have been routinely employed for quantitative description of different features of compounds such as topological and geometrical characteristics (2426). In this study, a total of 285 theoretical structural descriptors were calculated.

This set of descriptors was manually selected from descriptors calculated in Dragon software (version 2.1) by eliminating those descriptors that were obviously redundant or irrelevant to the prediction of pharmaceutical agents (2627).

The name and number of calculated theoretical descriptors included 12 descriptors in the class of functional group, 218 descriptors in the class of topological descriptors, 28 descriptors in the class of geometrical descriptors, and 27 descriptors in the class of constitutional indexes.

These theoretical descriptors were computed from the optimized 3D structure of each molecule, but not all of the descriptors are essential for the QSAR modeling.

In order to decrease interference of multicollinearity before model building, these molecular descriptors were preprocessed. The procedure includes: (i) exclusion of descriptors which have the identical value for more than 90% of the molecules; (ii) exclusion of descriptors with relative standard deviation less than 0.05; (iii) for each pair of descriptors with Pearson correlation coefficient over 0.9, only one descriptor, which has the higher correlation with the biological activity, remained (28).

Principal component analysis

Calculated descriptors were exported to the MATLAB environment for the purpose of PCA. The complete data set, defined by the descriptors in the columns (in this study, 285 descriptors) and the molecules in the rows, was auto scaled through mean centering by column. PCA models the maximum directions of variation in a data set by projecting the molecules as a swarm of points in a space spanned by PC's. Each PC is a linear function of a number of original descriptors, resulting in a reduction of the original number of descriptors. PC's describe, in decreasing order, the most variation among the molecules, and because they are calculated to be orthogonal to one another, each PC can be interpreted independently. This allows an overview of the data structure by revealing relationships between the molecules as well as the detection of deviating molecules (29). To find these sources of variation, the original data matrix of descriptors, defined by X(n,m), is decomposed into the molecule space, the descriptor space, and the error matrix. The latter represents the variation note explained by the extracted PC's and is dependent on the problem definition. The approach describing this decomposition is presented as:

where, X is the independent descriptor matrix, T is the scores matrix, P is the loadings matrix, E is the error matrix, n is the number of molecules, m is the number of descriptors, and k is the number of PC's used (29).

In a PCR analysis, a model formation step was carried out with stepwise selection and elimination of PCs to model the binding pIC50 relationships with different vectors of scores. On the other hand, in PCR procedure, all calculated scores were collected in a single data matrix and the best subset of descriptors was obtained by stepwise regression.

Genetic algorithm

The GA is used to select the descriptors that are most significant for the molecular data set. GA is a stochastic optimization technique that has been inspired by evolutionary principles (3031). One of the most important aspects of GA is that it studies many possible solution simultaneously, each of which explores different regions in space spanned by input variable (32). In this study, GA was tried for selecting significant PC's. In this way of using GA, an individual in the population is represented by the string of bits that encoding the selected descriptor. The first step in a GA is to create a gene pool of n individuals. Each individual (chromosome) contains some PC's that in the first generation these PC's are selected randomly from a table including calculated PC's and in a way such that no two individuals can be found that contain exactly the same set of descriptors (33). This individual was used as an input of ANN. The fitness score of each individual in this generation is determined by mean square error (MSE) of the network.

In the next step regeneration happen, so that the new offspring contains characteristics from both of its parents. Two individuals are selected probabilistically on the basis of their fitness scores and serve as parents. The selection strategy that used in this program was random selection method. Next step is applying a crossover operator that each parent contributes a random selection of half of its descriptors and an offspring is built by combining these two halves of genetic code. At last, this offspring is subjected to a random mutation in one of its gene, i.e. one variable is replaced by another. This selection–crossover–mutation procedure is repeated until all of the n parents in the gene pool are replaced by their offspring. The fitness score (RMSE and RMSECV) of each member of this new generation is again evaluated by using network, and the reproductive cycle is continued until a desired number of generation or target fitness score is reached. In the routine GA implemented in the MATLAB, some modifications were made (33). Here, for the computation of the fitness score of each chromosome a nonlinear model was built using descriptors consist in each chromosome separately by ANN technique and the values of MSE were calculated by means of this model. This process was used for each chromosome separately (3334). The basic design of the algorithm combined genetic algorithms and ANN is summarized in the flow diagram shown in Fig. 1.

Fig. 1

The basic design of the algorithm combined genetic algorithms and artificial neural network used in this study.

Feature selection with GA–PCA strategy

The GA–PCA based coupling technique used to select the suitable features was performed using our own written routines for GA and PCA in MATALB environment. These m-codes are optimization tools based on the GA strategy in MATALB. Steps were applied in this process can be described as follows: (1) definition and encoding of chromosome; (2) population initialization; (3) evaluation of each chromosome; (4) protection of chromosome; (5) selection of best chromosomes; (6) crossover and mutation operations; (7) stoping if a halt condition is satisfied, otherwise going to step 3 (22). More description about the theory of this GA strategy can be found in the literature (223536).

Artificial neural network

One method for providing a more flexible form of regression is to use a feed-forward neural network with error back-propagation learning algorithm. This is a computational system whose design is based on the architecture of biological neural networks which consists of artificial ‘neurons’ joined so that signals from one neuron can be passed to many others (23).

Clarification of the theory of the artificial neural networks in details has been adequately described elsewhere (37) but little relevant remarks is presented. ANN are parallel computational tools consisting of computing units named neurons and connections between neurons named synapses that are arranged in a series of layers (25).

Back propagation artificial neural network includes three layers. The first layer namely input layer has n neurons whose function is the reception of information (i.e. inputs) and their transferring to all neurons in the next layer called the hidden layer whose number is assined by n. The neurons in the hidden layer calculate a weighted sum of the inputs that is subsequently transformed by a linear or non-linear function. The last layer, the calculated response vector, is the output layer whose neurons handle the output from the network (38). The task of synapses is the connection of input layer to hidden layer and hidden layer to output layer. The manner in which each node transforms its input depends on the “weights” and bias of the node, which are modifiable. On the other hand the output value of each node depends on both the weight, and biases values. In addition, depend on the weighted sum of all network inputs which are normally transformed by a nonlinear or linear transform function determine the outputs of the network (39).

The relation between response, Y, of the network and a vector input, X, can be written as following if the number of neurons in the output layer is equal to 1 (same with our condition in here):

where, b is the bias term, W is the weight of the connection between the Ith neuron of the input layer and the Jth neuron of the hidden layer, and f is the transformation function of the hidden layer. In the training process, the weights and bias of the network which are the adjustable parameters of the network are determined from a set of objects, which is known as training set.

Through the training of the network, the connection weights are regulated so that error of calculated responses and observed values were minimized. For this, a nonlinear transfer function makes a connection between the inputs and the outputs. Commonly neural network is adjusted, or trained, so that a particular input leads to a specific target output.

There are numerous algorithms available for training ANN models. We used back propagation algorithm here for training of the network. In this algorithm several steps for minimizing of networks were performed and the update of weight for the ( pattern is given as:

By using following equation, the descent down the error surface is calculated (40):

where, α and μ are momentum and learning rate, respectively.

With respect to above description, some adjustable parameters including number of nodes in input and hidden layers, transfer function of the hidden and transfer function output layers, momentum (optimum value in this study was 0.16), number of iteration for training of the network (17000 epoch), and learning rate (0.84) evaluated by obtaining those which result in minimum prediction are present in the ANN.

As described above, in order to avoid overfitting or underfitting, a validation set was used in the ANN modeling. Evaluation of ANN was performed on an external set (validation set) that consisted of molecules belonging to neither the training set nor the test set (41).

RESULTS

A lot of descriptors were calculated for each studied molecule using Dragon. In order to calculate a relationship with independent variable, logarithms of the inverse of biological activity (Log 1/IC50) data of 49 molecules were used. After dividing the molecules into training, validation, and test sets, building of QSAR models using training and validation sets was carried out (42).

PCA was performed on the training set after deleting constant descriptors. Results of PCA considering 20 first PC's and also their eigenvalues are reported in Table 2. In this Table, the eigenvalues, the percent of variances explained by each eigenvalue and the cumulative percent of variances are reported (34). After acquiring PC's, linear regression with stepwise factor selection was carried out. The cross-validation technique used was eliminating only one molecule at a time and then performing PCR on the remaining of training set (leave one out technique). The activity of the left-out object was predicted by using this developed model. This process was repeated until each molecule in the training set had been gone out once.

Table 2

The result of principal component analysis on the total descriptors.

For evaluation of the predictive power of the generated PCR, the optimized model was applied for prediction of pIC50 values of all molecules in the training and test sets. The calculated pIC50 for each molecule and statistical error of prediction by model are summarized in Table 3. Experimental versus predicted values for pIC50 values of training and test set, obtained by the PCR modeling is shown graphically in Fig. 2A.

Table 3

The experimental pIC50 and the predicted values of the training set and test set.

Fig. 2

pIC50 estimated by modeling versus experimental values for training and test sets A; PCR, B; GA-PC-ANN.

This model was validated by some statistical parameters such as PRESS and RMSE reported in Table 4. It is clear that this model on the basis of statistical parameters is weak. With respect to these results, it is a good idea to try a nonlinear regression method, artificial neural network, to obtain robust and predictive QSAR model able to describe a relationship between the structure and P2X7 antagonist activity of the studied purine analogues.

Table 4

Statistical parameters obtained for the QSAR models.

In the non-linear model used, a network including a fully connected three layer, feed-forward ANN model trained with a back-propagation learning algorithm was used. The input of the network was the GA selected PC's. Seven PC's were selected by GA and were applied as input of the networks. In order to evaluate the ANN, MSE was used. The values resulting from hidden layer are transferred to the last layer, which contains a single neuron representing the predicted activity. For output layer a linear transfer function was chosen. Various ANN architectures were run with the seven selected PC's as input. In each run, the neuron architecture and parameters were optimized to reach the lowest MSE as the performances of the resulted models.

It must be considered that for inhibition of overfitting in the ANN model, the training of the network for the prediction of activity must be stopped when the MSE of the test set commences to increase while MSE of training set continues to decrease. Therefore, training of the network was stopped when overtraining began (43).

As mentioned above, before training the network, the number of nodes in the hidden layer must be optimized. For this purpose, several training of network was performed with different numbers of hidden nodes from 1 to 15. The MSE for training sets was obtained for different numbers of neurons at the hidden layer, and the minimum value of MSE was verified as the optimum value. A typical plot of MSE for training set versus the number of nodes in the hidden layer is shown in Fig. 3. It is clear that 7 nodes in the hidden layer is the optimum value.

Fig. 3

Plot of MSE for training sets versus the number of nodes in hidden layer.

The network was trained using training data and it was evaluated by prediction molecules. The predicted activity of the ANN calculated values of pIC50 are plotted against the experimental values in Fig. 2B and are reported in Table 3 and as expected, the calculated values are in good agreement with experimental values. The statistical parameters for the nonlinear model are represented in Table 4. As it is observed, the model obtained by the PC-GA-ANN has superior qualities relative to those obtained by PCR. This means that there is nonlinear relationship between the calculated PC's and the activity of the antagonists used in this study.

As a result, it was found that correctly opted and trained neural network could practically represent dependence of the activity of P2X7 receptor antagonist to the extracted PC's from various geometrical, topological, and other calculated descriptors. Then, the optimized neural network could simulate the complicated nonlinear relationship between pIC50 value and the PC's.

DISCUSSION

The statistical parameters obtained by developed ANN model for the training and test sets were reported in Table 4. These results demonstrated some significant differences between PC-GA-ANN and PCR models. It can be seen from this Table that statistical results of the ANN model is better than other method. Also these results reveal that GA is superior method for feature-selection in this QSAR study.

All QSAR models have to constantly be verified for their applicability of domain (AD), in order to generate reliable predicted activities for compounds that are not too structurally dissimilar (44). The Williams plot verified the presence of outliers (i.e. molecules with standardized residuals greater than two standard deviation units) and the compounds that were very influential in determining model parameters i.e. molecules with high leverage value (h) (4445) greater than 3 (m + 1)/n, where m is the number of the variables used in model formation, and n the number of the molecules employed to calculate the model. Also the data predicted by the models were verified for reliability by their leverage; with the intention that only predicted activities for molecules belonging to the chemical domain of the training set would be proposed (44). Actually, the leverage can be employed as a quantitative measure of the model applicability domain appropriate for assessing the degree of extrapolation: it represents a sort of compound “distance” from the model experimental space. The prediction for a compound having a high leverage value (h>h*, the warning limit leverage (WLL) h*= (3m/n) must be considered as unreliable (44). On the other hand, when the molecule has a leverage value lower than the WLL, the chance of accordance between the predicted and the experimental values is as high as that for the training set compounds (46). In this work each developed QSAR model was verified for the AD of the studied compounds to validate the prediction reliability.

Fig. 4A and 4B shows the AD of compounds of PCR and GA-PC-ANN employed in this study. Compounds influencing the structural domain of the model can be defined as molecules characterized by unusual structural features and badly represented in the training set. On the other hand, the outliers, whose standardized residual values exceed the cut off value of 2 standard deviation units, could be associated with the experimental error (447).

Fig. 4

Applicability domains of developed A; PCR and B; GA-PC-ANN.

By examining the AD of the developed models (1) from the Williams plot (Fig. 4), it can be seen that two compounds for PCR model (molecule 28 and 47) and three compounds (molecule 34, 42 and 8) for GA-PC-ANN model are identified as a response outlier based on the 2σ rule for the training set. On the basis of Williams plot leverage method neither of the molecules in the studied set is identified as structurally influential chemicals. According to the results presented in the Williams plot, it is evident that neither of studied compounds is outlier.

CONCLUSION

A novel approach is suggested to select features from high-dimensional data in data mining. A GA procedure is used to search for the best subset in a PCA matrix and results is combined with ANN for P2X7 receptor antagonist activity prediction. The performance of the approach was validated on the P2X7 receptor inhibitory data by comparison of PC regression methods. The obtained results represent that the suggested approach (combination of GA and ANN) leads to a better subset of variables than original linear PCR method. The effectiveness of the GA is explained by the selection of the best set of PC's. The main aspect of developed model based on a neural network is its ability to allow for fiexible mapping of the selected PC's by manipulating their functional dependence implicitly, unlike regression analysis. Models based on ANN handle both linear and nonlinear relationships between dependent and independent variables without adding complexity to the model. This capability offset the larger computing time required by a neural network simulation. This study shows that combination of PCA and artificial intelligence methods are potentially helpful for the prediction of P2X7 receptor inhibitory activity from a set of inhibitors. Also, it indicates that a non-linear method such as ANN is superior to a linear method such as PCR in building prediction models. In addition, the analysis of these GA selected PC's in the developed model can provide helpful clues to the structural and physicochemical characteristics of molecules contributing to the P2X7 receptor inhibitory activity; this may help to provide reference information for ligand-based P2X7 receptor antagonist drug design.