A quantitative structure-activity relationship (QSAR) study of some diaryl urea derivatives of B-RAF inhibitors.

In the current study, both ligand-based molecular docking and receptor-based quantitative structure activity relationships (QSAR) modeling were performed on 35 diaryl urea derivative inhibitors of V600EB-RAF. In this QSAR study, a linear (multiple linear regressions) and a nonlinear (partial least squares least squares support vector machine (PLS-LS-SVM)) were used and compared. The predictive quality of the QSAR models was tested for an external set of 31 compounds, randomly chosen out of 35 compounds. The results revealed the more predictive ability of PLS-LS-SVM in analysis of compounds with urea structure. The selected descriptors indicated that size, degree of branching, aromaticity, and polarizability affected the inhibition activity of these inhibitors. Furthermore, molecular docking was carried out to study the binding mode of the compounds. Docking analysis indicated some essential H-bonding and orientations of the molecules in the active site.

INTRODUCTION

RAF is among tyrosine kinase type receptors with serine/threonine kinase activity (1). Its contribution is in mitogen activated protein kinase (MAPK) signaling pathway, which conducts signals from membrane-based receptors to the nucleus to mediate cell proliferation, differentiation, and survival (2). Numerous cancers are related to the constitutive activation of the above signaling pathway (3). B-RAF is one of the isoforms of the RAF kinase family that can regulate multiple downstream molecules and is also regulated by a variety of signaling molecules (45). In about 7% of human cancers, the mutation of B-RAF has been detected (678910).

Some small molecule RAF kinase inhibitors by diverse scaffolds such as ureas, urea bioisosteres, imidazoles, benzamides, oxindoles, and aza-stilbenes have emerged in the recent past decades (1112). But diaryl urea have been most extensively investigated because of sorafenib success in clinical for renal and hepatocellular carcinoma (131415).

It is of great importance to introduce computer-aided drug design (CADD) approach to accelerate the time-consuming process of conventional drug discovery (16). Quantitative structure activity relationships (QSAR) and molecular docking are two of the helpful methods of CADD for drug design and prediction of drug activity (1718). In QSAR large number of compounds are usually evaluated resulting in models that can predict the potency or activity of new or even non-synthesized compounds (19).

When the three-dimensional structure of the target protein is available or can be modeled, molecular docking is often used for screening of compound libraries. Molecular docking predicts the conformation of a protein-ligand complex and calculates the binding affinity and investigates protein–ligand interactions (2021).

In this study aimed to develop a robust and accurate model for the inhibitory activity of inhibitors in order to design potential B-RAF kinase inhibitor. We used different method to connect between structural parameters and B-RAF kinase inhibitory. These methods included multiple linear regressions (MLR) as linear method and partial least squares least squares support vector machine (PLS-LS-SVM) as a nonlinear approach. The latter method was used to carry out non-linear mappings on the physicochemical and biological descriptors of the molecules. In Support vector machines, nonlinear kernel based functions were used to solve both regression and classification problems. An advantage of this method is its reproducibility in data mapping (22). Our aim in this study was to develop more examples of modeling based on this approach. Finally docking study was performed to suggest a binding mode for the inhibitors on B-RAF target.

MATERIALS AND METHODS

All calculations were made using an Intel Core i5 (CPU 2.6 GHZ) laptop running on windows 7 operating system. The m-files for MATLAB calculations were developed in our group.

Dataset and descriptor generation

The dataset used in this study was taken from the work of Menard, et al. (23). Chemical structure of 35 studied compounds is provided in Table 1. This set contains diarylurea derivatives with inhibition potency against B-RAF kinase. The chemical structures of molecules were drawn and optimized by HyperChem 7.0 software (HyperCube Inc. USA). Energy minimizations for all compounds were performed by AM1 semi-empirical method with Polark-Ribiere algorithm until the root-mean-square gradient of 0.01 Kcal/mol was reached. The resulted geometries were transferred into Dragon program (developed by Milano Chemometrics and QSAR Group) to calculate descriptors. The physicochemical parameters were calculated utilizing HyperChem and Dragon softwares. Molecular descriptors computed using the Dragon software were constitutional, functional, topological, and geometrical groups. Hyperchem was used to obtain chemical descriptors such as Log P, hydration energy, polarizability, molar refractivity, molecular volume, and molecular surface area. Gaussian 98 W package was employed to use HF method at 6-31G* basis set for optimization and calculation of different quantum chemical descriptors including dipole moment, local charge on atoms, high-occupied molecular orbital (HOMO), and low-unoccupied molecular orbital (LUMO) energies. Indices of electronegativity, electrophylicity, hardness, and softness were calculated from the energies of HOMO and LUMO.

Table 1

Chemical structure of B-RAF kinase inhabitor in this study

Model building

The correlation of the calculated descriptors with each other was calculated and collinear descriptors (0.85) were specified. Those with higher correlation towards activity vector were retained and the others were eliminated.

Splitting of the matrix into calibration (train) and external (test) set was done using kenard-stone algorithm. Subsequently, stepwise multiple linear regressions and partial least square analysis were performed on the training set for MLR and support vector machine (SVM) methods. In case of PLS-LS-SVM, Gaussian RBF Kernel with two tuning parameters, γ (gama) and σ2 (sigma2) were used. Latent variables of partial least squares (PLS) were entered to the models using stepwise method. The predicted residual error sum of squares (PRESS) was used to select the best models for validity evaluations.

Validation of QSAR models

To evaluate robustness of the models in the calibration subset, leave-one-outcross validation method was used. PRESS and Q2LOO values were thereafter calculated according to the equations 1 and 2.

where, SSY is the sum of squares of the response values. The models obtained from calibration set were then used to predict the response values of the test set. PRESS and R2 of test set were calculated as described above.

Finally, chance correlation test was used by scrambling the response vectors for several times and recalculating error metrics (PRESS and R2). Applicability domain for the best model was also done using Williams plot. For this purpose standard residuals were plotted against leverage for the data set (24). Other parameters of validation including K and K’ were also provided for both modeling approaches.

Molecular docking

Docking was conducted using AutoDock 4.2 software package. For this purpose, the crystal structure of B-RAF was retrieved from Protein Data Bank (PDB code: 1UWH). Grid and docking parameter files was generated by AutoDockTools version 6.5. For the preparation of protein, the original ligand and water molecules were removed from the coordinates and Kollman charges together with polar hydrogens were added. For ligands, Gasteiger charges were assigned and non-polar hydrogens were merged. The grid map of 60 × 60 × 60 points in x, y, and z directions was centered on the binding site (Cys 531) with a spacing of 0.375 Å. Lamarckian genetic algorithm (LGA) was used for docking with the following settings: a maximum number of 25,000,000 energy evaluations, an initial population of 300, a maximum number of 27,000 generations, and 100 independent docking runs. Other parameters of. dpf file were remained as default. All parameters of docking were validated based on re-docking of the cognate ligand (BAY 43-9006) and measuring the root-mean-square deviation (RMSD) value between the best pose of the docked structure with its conformation at crystal state. The RMSD value of less than 2 Ε can validate the protocol of docking. Results were ranked according to the binding free energy and mode of interactions. Ligplot program was used to visualize hydrophobic and hydrogen-bonding interactions between the ligand and receptor (25).

RESULTS

Sixty two remained descriptors constructed an input file for stepwise selection based on MLR analysis. Comparison of statistical parameters was employed to result the best model as seen in the models bellow (Table 2). The numbers of training set compounds in the first and second model were 31 and 28, respectively.

Table 2

Comparison of statistical parameters of the models

pIC50 = -34.171 (± 5.597) - 0.006 (± 0.002). αzz - 0.091 (± 0.009). G (N...N) + 3.260 (± 0.379). TI2 - 0.110 (± 0.016). DISPm + 22.523 (± 10.470). PW3 + 7.682 (± 1.542). BLI + 50.35 (± 17.976). PW4 + 2.055 (± 0.750). (PJI3). (3)

where, pIC50 is -log (IC50), N = 31, F = 26.102, R2 = 0.905, Se = 0.3131.

pIC50 = -9.576 (± 8.537) + 0.062 (± 0.12). G (N...N) + 1.981 (± 0.891). DISPm + 66.775 (± 12.261). BLI - 34.512 (± 10.247). (FDI). (4)

where, N = 28, F = 13.441, R2 = 0.837, Se = 0.5144.

In this equation, the values in the parentheses show the standard deviation of the coefficients. N, R, Se, and F are numbers of components, coefficient, standard error of regression, and Fisher's F-ratio, respectively.

As a result for PLS-LS-SVM, the first three latent variables resulted in a model with reasonable R2 values (0.93 and 0.84).

Internal validation of the models was revealed by means of leave-one-out cross-validation. The results for cross-validation studies of both modeling approaches are listed in Table 2.

The more robustness of LS-SVM has been therefore shown based on both PRESS and R2 values. Furthermore, LS-SVM was more resistant to the number of compounds in training and test set with respect to MLR. Chance correlation test by means of Y-permutation ensured that the presented models were not achieved randomly.

As depicted in Fig. 1a, the molecules were not more than the cut-off h value of leverage. Applicability domain studies for the model based on PLS-LS-SVM have revealed that the model can be applied for structurally relevant compounds.

Fig. 1

(a) applicability domain based on williams plot, (b) interaction model of docked inhibitor 32 on B-RAF kinase.

Docking analysis revealed important H-bonding interactions between urea moiety with Glu500 and Asp593 residues in allosteric pocket. H-bonding between tow nitrogen atoms of imidazopyridine moiety with Cys531 residue of hinge binding region was also significant.

DISCUSSION

R2, Q2, and F ratio of two models were compared, model 1 revealed the best statistical parameters R2 = 0.905, Q2 = 0.8145, F = 26.102 (Table 2). The more power of LS-SVM compared to MLR has been verified based on external validation including K and K’ as displayed in Table 3. SVM model has passed the validation criteria and showed to be a more powerful modeling approach.

Table 3

The results for QSAR modeling using MLR and PLS-LS-SVM

The best MLR model contained 8 descriptors including 4 topological parameters (TI2, BLT, PW3, PW4), 3 geometrical parameters (G(N..N), PJI3, DISPm), and one chemical parameter (αzz). Based on descriptor coefficients, increasing αzz, G (N...N), and DISPm can decrease pIC50 and increasing TI2, PJI3, BLI, PW3, and PW4 can increase the extent of pIC50 of the structures.

TI2 is a topological descriptor belonging to the Mohar indices which can affect solubility. TIs can be sensitive to size, shape, symmetry, branching, and cyclicity.

G (N...N) is a 3D geometrical descriptor relating to the geometrical distance between nitrogen and nitrogen in molecules with more than one N-atom.

DISPm is related to the molecular geometry as well as molecular size and describes the conformational features of the molecules. Negative contribution of DISPm revealed that rigidity of molecules can increase inhibitory activity due to the absence of long and flexible substituents and the presence of unsaturated bonds. Kier benzene-likeliness index (BLI) is an aromaticity index calculated from molecular topology. Aromaticity might influence the transport or distribution of the chemical across cell membrane.

PW3 and PW4 as the path/walk count ratio are independent of molecular size. These descriptors can be considered as shape descriptors whose value increases with increased branching in the vertices.

PJI3 (Petitjean shape indices), a first Petitjean shape index, is a topological anisometry descriptor. αzz is the polarizability along zz axis of the molecule which represents how readily the molecular charge distribution is distorted by external electric fields. Almost 7 descriptors appeared in the SMLR equation belongs to the geometrical and topological groups of descriptors.

In Support vector machines, nonlinear kernel based functions are used to solve both regression and classification problems. To reduce the complexity of optimization, least squares support vector machines has emerged in the recent years. In the following regression equation, ω and b are the slope and the intercept of the equation, respectively.

F(x) = ωTφ (xi) + bi = 1, 2,… (5)

The function φ(xi) is used to map the input space into a higher dimensional space.

To investigate the role of descriptors on the first three variables, correlation of all descriptors with the loading values of the three variables were calculated. It was observed that the factors 1, 2, and 3 have higher loadings for TE, Ms, and FDI respectively. This result was found in both modes with 28 and 31 compounds. Total energy (TE) of a molecular system is among quantum chemical descriptors and is the sum of the total electronic energy (Eee) and the energy of intern clear repulsion (Enr).

TE = Eee + Enr (6)

The mean electro topological state (Ms) is among constitutional descriptor which is calculated by dividing Ss by the number of non-hydrogen atoms (nSK):

Ms = Ss / nSK (7)

Folding degree index (FDI) is a geometrical descriptor. Negative contribution of TE revealed that decrease in total energy for molecular systems can increase their inhibitory activity. Whereas positive values of Ms and FDI descriptors show that the indicated descriptor contribute positive effect on the values of pIC50. According to experimental QSAR results and docking binding energies performed in this study, the modification of the phenyl middle ring using hydrophobic group such as halogen, thiomethyl substitution (2,6 compounds) or replacement by naphthyl group (24-35, 22, 19, and 5 compounds) can enhance the inhibitory activity due to an increase in the potency of binding and this result was in accord with the previously reported SAR result (23).

The binding energy of all compounds is reported in Table 4 and the compound with naphthyl moiety have the best binding energy. In Fig. 1b, orientation of the interactions for compound 32 with the lowest binding energy is depicted.

Table 4

Binding energy resulted from Autodock software

CONCLUSION

In this work, molecular docking and structure-based QSAR method were performed to explore structural features, orientations, and conformations of some inhibitors interacting with B-RAF kinase. A structure-based QSAR model using MLR and LS-SVM were developed. The result indicated that the performance of LS-SVM model was better than the MLR model.

The selected descriptors indicated that size, degree of branching, aromaticity, and polarizability affected the inhibition activity of these inhibitors. Docking analysis indicated some essential H-bonding and orientations in active site. This information will guide us through design of novel B-RAF kinase inhibitors.