Comparative QSAR Analysis of 3,5-bis (Arylidene)-4-Piperidone Derivatives: the Development of Predictive Cytotoxicity Models

Introduction

The principal aim of this laboratory is the synthesis of conjugated unsaturated ketones as candidate antineoplastic agents. These compounds interact with thiols but in general, they have little or no affinity for amino and hydroxyl groups which are found in nucleic acids (1-3). Hence thiol alkylators may not have the genotoxic properties associated with a number of contemporary anticancer drugs (4).

Initially various series of compounds possessing one conjugated enone group were developed. However after an initial chemical insult, certain neoplasms are more vulnerable to a subsequent cytotoxic effect than various non-malignant cells (5, 6). Hence by mounting the 1,5-diaryl-3-oxo-1,4-pentadienyl pharmacophore [AR-C=C-C(O)-C=C-AR] on heterocyclic and cycloaliphatic scaffolds, two sequential interactions with cellular thiols can take place which may be more detrimental to tumours than normal tissues. Such considerations led to the development of 3,5-bis(benzylidene)-4-piperidones 1a-d which demonstrated potent cytotoxic properties with the IC₅₀ values in the low micromolar range against human Molt 4/C8 and CEM T-lymphocytes as well as murine L1210 lymphocytic leukemia cells (7, 8). The hypothesis of sequential cytotoxicity was advanced that the 1,5-diaryl-3-oxo-1,4-pentadienyl group interacts at a primary binding site and a side chain on the piperidine nitrogen may align at an auxiliary binding site which could enhance cytotoxic potencies. In order to evaluate this hypothesis, a novel series of N-aroyl-3,5-bis(benzylidene)-4-piperidone derivatives 2-7 were synthesized (8). In these compounds, the side chain contains atoms and groups capable of forming hydrogen and van der Waals bonds and also an ionic bond in series 3-6 which could bind at an auxiliary binding site. When the potencies of the compounds in series 3-6 were compared with that of the analogs 1a-d which have the same aryl substituent, the IC₅₀ values of the amides 3-6 were lower in 48% of the comparisons made while equipotency was noted in 35% of the cases (8). A number of these amides have submicromolar IC₅₀ values.

The proposed primary and auxiliary binding sites of cytotoxic N-aroyl-3,5-bis(benzylidene)-4-piperidone analogs

Figure 1

The objectives of the present investigation are two folds. First, to gain some insight into the physicochemical properties of the 3,5-bis(benzylidene)-4-piperidone derivatives 1-7 that influence cytotoxic potencies, and second, to derive predictive 2D-QSAR models that will be helpful to design new analogs with improved bioactivity.

Experimental

Data set

The in-vitro cytotoxic data used in this study were taken from the literature (8). The reported IC₅₀ values were converted into logarithmic scale (pIC₅₀) and then were used for QSAR analysis as dependent variables. The data set (n = 25) was divided randomly into two groups: calibration set (n = 19) and prediction set (n = 6). (Table 1) (Asterisk sign in the table demonstrate the prediction set).

Descriptor generation

The chemical structure of molecules was constructed using a Hyperchem package (Version 7; Hypercube Inc., Gainesville, FL, USA) and primarily were energetically minimized (100 steepest descent steps using MM⁺ force field with a gradient convergence value of 0.1 kcal/mol) using HyperChem software. The Z-matrices of the structures were provided by the software and were then transferred to the GAUSSIAN 98 program. (9). Complete geometry optimization was performed based on most extended conformations as the starting geometries. GAUSSIAN 98 was used for semi-empirical molecular orbital calculations (AM1) of the structures.

A large number of molecular descriptors were calculated using HyperChem, GAUSSIAN 98 and Dragon Packages (10). Different quantum chemical descriptors were calculated using GAUSSIAN 98. Different functional descriptors such as topological, geometrical, and constitutional descriptors were calculated using Dragon software. Moreover, some chemical parameters including molecular volume (MV), molecular surface area (SA), hydrophobicity (logP), heat of formation and hydration energy (HE) were calculated using Hyperchem software.

Data processing and modeling

The calculated descriptors were collected in a data matrix and generated descriptors were scaled using auto scaling procedure. Firstly, the descriptors were checked for constant or near-constant values and those detected were removed from the original data matrix. Then the correlation of descriptors with each other’s and with the potency data was determined. Among the detected collinear descriptors (r² > 0.8), one of them demonstrating the highest correlation with potency was retained, and the rest were omitted. In the model development step, leave-one-out cross-validation was used to optimize the models (e.g., selection of variables or no. of principal components). The final models were validated using a set of 6 molecules in the test set, which did not contribute to the model development.

Two different methods were used for the development of QSAR equations (i) stepwise multiple linear regression (MLR), (iii) genetic algorithm–partial least squares (GA-PLS). These methods are well documented in the literature (11).

In stepwise regression, a multiple-term linear equation was built step by step. The basic procedures involve (i) identifying an initial model, (ii) iteratively 'stepping', that is, repeatedly altering the model at the previous step by adding or removing a predictor variable in accordance with the 'stepping criteria' (in our case, probability of F = 0.05 for inclusion; probability of F = 0.1 for exclusion for the forward selection method), and (iii) terminating the search when stepping is no longer possible given the stepping criteria, or when a specified maximum number of steps have been reached. Specifically, at each step, all variables are reviewed and evaluated to determine which one will contribute most to the equation. This variable will then be included in the model, and the process starts again. A limitation of the stepwise regression search approach is that it presumes there is a single 'best' subset of X variables and seeks to identify it. There is often no unique 'best' subset, and all possible regression models with a similar number of X variables as in the stepwise regression solution should be fitted subsequently to study whether some other subsets of X variables might be better. In the present study, MLR with stepwise selection and elimination of variables was applied for developing QSAR models using SPSS 11.5 software. The resultant models were validated by a leave-one-out cross-validation procedure (using MATLAB software) to check their predictive power and robustness.

A genetic algorithm uses genetic rules such as reproduction, cross-over, and mutation to build pseudo organisms that are then selected, on the basis of a fitness criterion to survive and pass information on to the next generation. To select the most relevant set of descriptors, the evolution of population was simulated (12-15). PLS is a generalization of regression, which can handle data with strongly correlated and ⁄ or numerous X variables (15). It gives a reduced solution which is statistically more robust than MLR. The linear PLS model finds 'new variables' (latent variables or X scores) that are a linear combination of the original variables. To avoid over fitting, a strict test for the significance of each consecutive PLS component is necessary and then stopping when the components are non-significant. Cross-validation is a practical and reliable method of testing this significance (14). Application of PLS thus allows the construction of larger QSAR equations while avoiding over fitting and eliminating most variables. PLS is normally used in combination with cross-validation to obtain the optimum number of components (15)19]. In the GA-PLS procedure, the optimum number of latent variable must be determined. Here, for each subset of descriptors (i.e., for each chromosome of the GA), a PLS model was developed separately and consequently the number of latent variables was optimized. The PLS regression method used was the NIPALS-based algorithm found in the chemometrics toolbox of MATLAB software (version 6.5; Math work Inc., Natick, MA, USA). Based on the Haaland and Thomas F-ratio criterion we performed the leave-one-out cross-validation procedure as it was applied in our previous article (16). The MATLAB PLS toolbox developed by eigenvector company was used for PLS and GA modeling.The predicted power of obtained equations was validated using cross-validated squared correlation coefficient (Q2). In order to calculate Q² following equation was used.

Q² = 1 − Σ(Y_pred − Y_act)2/ Σ (Y_act − Y_mean)

Where Y_pred, Yact, and Ymean are predicted, actual, and mean values of the pIC₅₀, respectively. Σ(Ypred - Yact)² is the predictive residual error sum of squares (PRESS) (17). Moreover, the predictive squared correlation coefficients “R² _pred” was used for further evaluation of the predictive potential of obtained equations. To calculate the R² _pred following equation was used.

R² _pred = 1 − Σ(Y_pred(Test) − Y_Test)²/Σ(Y_Test − Y_Training)

where _YPred(Test) and Y_Test are predicted and observed activity values, respectively, of test set compounds, and Y_Training is the mean activity value of training set (17).

Variable importance in the projection

In order to evaluate the relative importance of the variable appeared in the final GA-PLS, variable importance in projection (VIP) was employed (18). VIP values reflect the importance of terms in the PLS model. VIP is the sum over all model dimensions of the contributions variable influence (VIN). For a given PLS dimension, A (VIN)ak² is equal to the squared PLS weight (w_ak)2 of that term, multiplied by the explained SS of that PLS dimension (19).

Y-randomization test

Y-randomization test is the technique ensures the robustness of a QSAR model. In this procedure, for each created models, the dependent variable vector pIC50 is randomly shuffled, and a new QSAR model is developed using the original independent variable matrix. For acceptable specific QSAR model, the new obtained QSAR equation (after several repetitions) is expected to have low R² and q² values (16).

Results and discussion

The structural features and the experimental cytotoxic potencies (represented as pIC₅₀) of the 3,5-bis(benzylidene)-4-piperidone derivatives 1-7 used in this study are shown in Table 1 (8). A SAR study of these compounds demonstrated that the attachment of the N-aroyl group to the 3,5-bis(arylidene)-4-piperidones resulted in enhanced cytotoxic potencies toward murine L1210 lymphocyte leukemia cells and two human T-lymphocytes Molt4/C8 and CEM. It was suggested that in the case of various members of series 2-7, alignment of the N-aroyl groups with auxiliary binding sites takes place reinforcing the interaction of the 1,5-diaryl-3-oxo-1,4-pentadienyl pharmacophore with a primary binding site.

In order to find important structural features of the N-aroyl-3,5-bis(benzylidene)-4-piperidones which contribute to the cytotoxic potencies, a QSAR study was performed using different molecular descriptors and chemometrics tools. In this study, the chemometrics methods such as stepwise multi-linear regression (MLR) and genetic algorithm-partial least squares (GA-PLS) were used for modeling the relationship between the biological activity and molecular descriptors and the results were compared to find the best predictive cytotoxicity models.

MLR modeling

Separate stepwise selection-based MLR analyses were performed using different types of descriptors, and then a MLR equation was obtained using the pool of all calculated descriptors. As there are 19 molecules in the training set and according to the rule of thumb (the ratio of 5:1 for molecule ⁄ variable; Toplis ratio) (16), MLR models with a maximum number of 3 or 4 variables were selected. The results are summarized in Table 2. A small difference between the conventional and cross-validated correlation coefficients of the different MLR equations (Table 2) reveals that none of the models are over fitted, which can be partially attributed to the absence of collinearity between the variables in one hand and using no extra variables on the other hand. The correlation coefficient (r²) matrix for the descriptors used in different MLR equations shows that no significant correlation exists between pairs of descriptors (Table 3).

L1210 cell line. Stepwise selection and elimination of variables produced a three-parametric QSAR equation 1, having moderate statistical quality (R² = 0.70, SE = 0.38, Q² = 0.50, RMS_CV = 0.45 and R²_pred = 0.69) (Table 2 and 4). The selected variables demonstrate that chemical (heat of formation), topological and spatial parameters (molecular density and Narumi harmonic topological index (HNar) mainly affect the cytotoxic potential of the 3,5-bis(arylidene)-4-piperidone derivatives against the L1210 cell line.

CEM cell line

The MLR equation 2 obtained from the pool of the calculated descriptors for the CEM cell line possesses good statistical quality (R² = 0.85, SE = 0.25, Q² = 0.77, RMS_CV = 0.29 and R²_pred = 0.81) (Table 2 and 4) and demonstrates that the anticancer property of the compounds is mainly dependent on chemical (heat of formation), topological and spatial (molecular density and average connectivity index chi-4:X4A) properties of the molecule. These results show that the parameters that influence anticancer properties of the compounds against the CEM line are similar to those obtained for the L1210 cell line.

Molt/C8 cell line

The three-parametric QSAR equation 3 (Table 2) demonstrates the quantitative effects of the structural parameters on the cytotoxic potential of 3,5-bis(arylidene)-4-piperidone derivatives on the Molt/C8 cell line. This model demonstrated high statistical quality (R² = 0.81, SE = 0.31, Q² = 0.72, RMS_CV = 0.34, and R²_pred = 0.64) (Table 2 and 4) which reveals the significant effects of quantum (HOMO, dipole moment X (DMx)) and spatial parameter (molecular density) on the cytotoxic activity of the compounds. It is revealed that the molecular density plays a determinant effect on the cytotoxic potential of the compounds on this cell line. Therefore, the molecular density and topological parameter are important descriptors that should be taken into consideration in the designing of potent cytotoxic 3,5-bis (arylidene)-4-piperidones analogs.

GA-PLS modeling

L1210 cell line

In PLS analyses, the descriptors data matrix is decomposed into orthogonal matrices, the scores of which are constrained to have inner relationships with the dependent variables. Therefore the multi-collinearity problem in the descriptors is omitted by PLS analysis. To find a more convenient set of descriptors in PLS modeling, a genetic algorithm was used. To do so, many different GA-PLS runs were conducted using different initial sets of populations. Given 19 calibration samples, the leave-one-out cross-validation procedure was used to find the optimum number of latent variables for each PLS model. For the L1210 cell line, the GA-PLS model that resulted in the best fit contained 7 indices (two of these indices such as molecular density and heat of formation were also obtained by MLR in this cell line). Moreover, similar to the results of MLR analysis, GA-PLS analysis also showed that topological parameters affect the cytotoxic potencies of the compounds in the L1210 cell line. As per GA-PLS modeling, a combination of chemical (heat of formation, melting point), quantum chemical (dipole moment X), conformational (torsion energy) and spatial (molecular density) descriptors account for the cytotoxic potential of the compounds toward L1210 cells. The resulted GA-PLS model possess very good statistical quality (i.e., R² = 0.86, Q²= 0.66). The predictive ability of the model was measured by application to 6 external test set molecules (R²_pred = 0.81 and SE_pred =0.37). The calculated values of pIC₅₀ obtained by the PLS model in the L1210 assay (refined from cross-validation of external prediction set) are shown in Table 1.

To measure the significance of the 7 selected PLS descriptors, the variable importance in projection (VIP) was calculated for each descriptor. According to the method described by Eriksson et al. (20), X-variable (predictor variable) could be classified according to their relevance in explaining Y-variable (predicted variable) i.e. VIP > 1.0 (highly influential), VIP < 0.8 (less influential), and 0.8 < VIP < 1.0 (moderately influential). The VIP analysis of the descriptors of the input variables used in the PLS equation (Figure 2A) shows that the parameters such as topological index (Chi3), molecular density and heat of formation are the most important indices in the QSAR equation derived by PLS analysis for L1210 cell line. The regression coefficients of the selected variables obtained by GA-PLS model of QSAR analysis for three different cell lines are illustrated in Figure 3.

Variable importance for projection (VIP) of the selected variables obtained by GA-PLS model for the QSAR analysis of three different cell lines: (A) L1210, (B) CEM, and (C) Molt/C8

Figure 2

Regression coefficients of the selected variables derived by the GA-PLS model for QSAR analysis of three different cell lines: (A) L1210, (B) CEM, and (C) Molt/C8

Figure 3

CEM cell line

The best GA-PLS model to predict cytotoxic activities of the compounds in series 1-7 toward the CEM cell line includes 10 indices (three of the indices are similar to that obtained by the MLR analysis). A combination of topological (X₁A), Chi-1(average connectivity index), PW3 (Randic shape index), X₃Av (average valance connectivity), SIC₂ (structural information content: symmetry 2-order)), geometrical (SPH (spherocity of molecule)), chemical (LogP), spatial (molecular density), electrotopological state (sum of the E-state value for the CH group in the aromatic ring (S_aaCH) and sum of the E-state value for the methyl group (S_sCH₃)) indices account for the cytotoxic properties of the 3,5-bis(arylidene)-4-piperidones toward the CEM cell line. The resulted GA-PLS model showed high statistical quality (i.e., R² = 0.87 and Q² = 0.71). The PLS estimate of the coefficients for these descriptors are given in Figure 3B. The predictive ability of the model was measured by application to 6 external test set molecules (R²_pred = 0.83 and SE_pred =0.35).

The pIC₅₀ values used in the PLS model (refined from cross-validation of external prediction set) are shown in Table 1. The VIP analysis of PLS equation presented in Figure 2B suggests that among the selected descriptors, topological and electrotopological parameters such as X₁A, PW₃, X₃Av and estate keys (S_sCH₃ and S_aaCH) are important influential descriptors. Similar to the L1210 cell line, the molecular density of the compounds is also a determinant parameter that influence cytotoxicity of 3,5-bis(arylidene)-4-piperidones against the CEM cell line.

Molt 4/C8 cell line

A combination of 8 selected variables such as chemical (hydration energy), topological (X₂A), geometrical (PW₃, J₃D (3D-balaban index), SPH), spatial (molecular density) and electrotopological state indices (sum of the E-state value for the CH group in the aromatic ring (SaaCH) and sum of the E-state value for the methyl group (S_sCH₃)) influence cytotoxic potencies of the compounds in series 1-7 against the Molt 4/C8 cell line as indicated by the best GA-PLS model. Most of these variables contribute significantly to the GA-PLS model as discussed previously for the CEM cell line. The resultant GA-PLS model demonstrated very high statistical quality [R² = 0.91 and Q² = 0.77] [Figure 3C]. The predictive ability of the model was measured by applying to 6 external test set molecules and the squared correlation coefficient for prediction (R²_pred) was found to 0.82.

The VIP analysis of the GA-PLS model presented in Figure 2C showed that the molecular density, topological (X₂A), geometrical (SPH, J₃D) and electrotopological state indices influence the cytotoxicity of the compounds in series 1-7. This result is similar to that obtained for the L1210 cell line.

The plots of the predicted pIC₅₀ values (cross-validated) against the experimental values for the L1210, CEM and Molt4/C8 cell lines by GA-PLS analysis showed R² values 0.86, 0.87 and 0.91, respectively (Figure 4) which indicates that the QSAR models are reliable to predict cytotoxicity of the 3,5-bis(benzylidene)-4-piperidones 1-7.

Plots of the cross-validated predicted cytotoxic potency values of activity by GA-PLS against the experimental values for three different cell lines: (A) L1210, (B) CEM, and (C) Molt/C8 cell line

Figure 4

All of the generated GA-PLS models were further validated by applying the Y-randomization test. Several random shuffles of the Y vector were performed and the results are shown in Table 5. The low R² and Q² values indicate that the good results in our original model are not because of a chance correlation or structural dependency of the training set.

The statistical parameters of QSAR analysis by MLR and GA-PLS are demonstrated in Table 4. The results indicate that GA-PLS analysis is more accurate for predicting the cytotoxic potential of 3,5-bis(arylidene)-4-piperidones in all three cell lines. The cross-validation statistics reported in Table 4 suggest the higher prediction ability of the GA-PLS model. This can be a result of the more number of descriptors used by GA-PLS with respect to MLR analysis. As far as the less parametric model is obtained by MLR analysis, the results of MLR analysis are more descriptive for interpretation of the structure-cytotoxic relationship. However for the prediction of the cytotoxic potential of novel compounds of these series, the GA-PLS model is more useful.

Conclusion

It is clearly understood that some of the spatial, (electro) topological and chemical descriptors are important structural parameters that mainly affects the cytotoxic potencies of these compounds in all three cell lines. The “molecular density” of compounds is one of the most important spatial descriptors that should be considered in the design of novel compounds of these series. As the results of the SAR study of various members of these series suggest, alignment of the N-aroyl groups with auxiliary binding sites takes place reinforcing the interaction of the 1,5-diaryl-3-oxo-1,4-pentadienyl pharmacophore with a primary binding site (8). This result can be interpreted by the important structural descriptor “molecular density”. It is obvious that the orientation of N-aroyl group for better interaction with a corresponding binding site is affected by molecular density. Other topological descriptors also affect this desirable orientation of the ligand and therefore play a determinant role in the binding of the ligands with its target receptors.

The other important electro-topological parameter that should be considered is the electro-topological state indices which encode the electronic state of each atom type in the molecule: “atom-type E-state indices”. The results demonstrate that the sum of the electro-topological state of “aromatic CH” (S_aaCH) and the sum of the electro-topological state indices of methyl groups (E-state keys: S_sCH₃) are among the important influential descriptors in the cytotoxic potential of 3,5-bis(arylidene)-4-piperidones especially toward CEM and Molt/C8 cell lines. The sum of the electro-topological state indices of the aromatic CH and methyl groups especially should be considered in the employment of different N-acyl substituents.

The results demonstrate that similar structural descriptors affect the cytotoxic potential of N-aroyl-3,5-bis(arylidene)-4-piperidones in three different cell lines. This result indicates that similar molecular mechanisms are involved in the cytotoxic activity of the compounds in these cell lines. In addition, the previous quantitative-cytotoxic activity study of asymmetrical 2,6-bis(benzylidene)cyclohexanones against MCF-7, SK-N-MC and MDA-MB-231 cell lines (20), reported the spatial and topological parameters as determinant structural descriptors in cytotoxic potential of these compounds against the cell lines. This finding also confirms that the same molecular mechanisms are involved in the cytotoxic activity of different analogues of this scaffold.