QSAR studies of novel 1-(4-methoxyphenethyl)-1H-benzimidazole-5-carboxylic acid derivatives and their precursors as antileukaemic agents

Peer review under responsibility of Taibah University

Abstract

In an attempt to find a potent inhibitor of antileukaemic activity, we performed two-dimensional quantitative structure–activity relationship studies. The QSAR models of the 1-(4-methoxyphenethyl)-1H-benzimidazole derivatives were developed by means of utilizing simulated annealing for the purpose of identifying an effective inhibitor. In the present study, the best 2D QSAR model developed had a cross-validated value of (q²) = 0.7584 and a coefficient of determination of (r²) = 0.8477 with partial least squares analysis. The models also suggest that the inclusion of fluorine, chlorine, methyl, and nitro substituents would enhance the antileukaemic activity. The resulting descriptors produced by the QSAR models were used to identify the physico-chemical features relevant to antileukaemic agents.

Keywords:

• Benzimidazole-5-carboxylic acid
• QSAR
• K562 cells
• CEM cells
• Simulated annealing
• Partial least squares
• Antileukaemic agents

1 Introduction

Recently, multi-target drugs, which are designed as single molecules to modulate multiple physiological targets simultaneously, have increasingly attracted the attention of medicinal chemists [1], and they represent a very promising way to enhance the efficacy and decrease the adverse effects of drugs, especially in the treatment of complex diseases such as cancers, cardiovascular diseases, and neurodegenerative diseases [2,3]. The leukaemias of childhood are common cancers of the hematopoietic system, primarily involving the malignant transformation of lymphoid progenitor cells and, less commonly, the transformation of myeloid progenitor cells [4]. Cancer is now believed to result from the unlimited growth of a given cell [5]. Cancer has been recognized as a disease of aberrant cellular proliferation, with traditional cancer therapies aiming to exploit the proliferation machinery. Cells die in two ways: necrosis and apoptosis. The regulation of apoptosis is crucial for development and sustained health [6]. The dysregulation of apoptosis will result in a variety of clinical disorders, including cancer. The benzimidazole ring displays an important heterocyclic pharmacophore in drug discovery. Compounds carrying different substituents in the benzimidazole structure are associated with a wide range of biological activities, including anticancer [7], antiviral [8], antibacterial [9], antifungal [10], anti-inflammatory [11], proton pump inhibitor [12], antioxidant [13], antihypertensive [14] and anticoagulant [15] properties.

Quantitative structure–activity relationship (QSAR) studies are mathematical equations that quantitatively correlate chemical structures with biological activity. These relationship models have been helpful in understanding the influence of molecular properties on the biological activity of different compounds, ultimately providing rational clues for the development of new compounds with desirable biological properties. Because they provide valuable information for molecular design and medicinal chemistry, QSAR studies have been widely used in drug design and discovery [16].

To gain insight into the structural and molecular requirement influencing antileukaemic activity, we herein describe the QSAR analysis of the novel compound 1-(4-methoxyphenethyl)-1H-benzimidazole, focusing on the development of a QSAR models for antileukaemic activity. These results should serve as a guideline in designing compounds with more potent antileukaemic activity.

2 Materials and methods

2.1 Computational details

The computational studies were performed on a computer using HP Windows 7 Home Basic with an Intel^® core processor. The molecular structures of the compounds in the data set were sketched using V-life MDS (Molecular Design Suite)™ 3.5 software, supplied by V-life Sciences Technologies Pvt. Ltd., Pune, India [17].

2.2 Obtaining biological data and training and test sets

The QSAR studies were performed using a series of substituted 1-(4-methoxyphenethyl)-1H-benzimidazole analogues reported in the literature [18]. Out of 22 molecules, four molecules were discarded for which the precise data were not available. The biological activities were represented by the IC₅₀ divided by the mean IC₅₀ in the two tumour cell lines K562 and CEM. The IC₅₀ values were converted into the negative logarithm of pIC₅₀ values and used as dependent variables in the QSAR analysis. The test compounds were selected considering the structural diversity and wide range of activity within data set. The chemical structures and their biological activities are represented in Table 1. The sphere exclusion method [19] was adopted for division of the training and test data sets comprising 13 and five molecules, respectively, with a dissimilarity value of 1.0, where the dissimilarity value gives the sphere exclusion radius. The compounds in the test set allowed us to use one test compound over four training compounds, thus resulting in a more rigorous validation of the training model. Five compounds, namely, 3, 7, 11, 16 and 20, were used as the test set, while the remaining molecules were used as the training set.

Table 1 The structures of 1-(4-methoxyphenethyl)-1H benzimidazole-5-carboxylic acid derivatives with their K562 and CEM cells. .

S. No	R	K562 IC50	K562(pIC50)	CEM IC50	CEM (pIC50)
1		3	5.5228	4	5.3979
2		78	4.1079	175	3.7569
3^a		120	3.9208	122	3.9136
4		180	3.7447	175	3.7569
5		136	3.8664	138	3.8601
6		70	4.1549	68	4.1674
7^a		84	4.0757	85	4.0705
8		>250	–	>250	–
9		56	4.2518	55	4.2596
10		>250	–	>250	–
11^a		64	4.1938	65	4.1870
12		135	3.8696	132	3.8794
13		180	3.7447	182	3.7399
14		250	3.6020	245	3.6108
15		>250	–	>250	–
16^a		200	3.6989	195	3.7099
17		40	4.3979	42	4.3767
18		80	4.0969	85	4.0705
19		110	3.9586	112	3.9507
20^a		240	3.6197	244	3.6126
21		124	3.9065	120	3.9208
22		>250	–	>250	–

a Indicates the compounds considered in the test set in QSAR.

2.3 Structure generation

The molecular structures of all of the molecules were sketched using the 2D builder module within VLife MDS 3.5 software, and the structures were then converted into the 3D space for further analysis. All of the compounds were batch optimized for the minimization of energies and optimization of geometry using Merck Molecular Force Field, followed by adjusting the parameters to include a distance-dependent dielectric constant of 1.0, convergence criterion or root mean square (RMS) gradient at 0.01 kcal/mol Å and an iteration limit of 10,000 [20].

2.4 Descriptor calculation

A large number of theoretical 2D individual descriptors such as Chi, Path-cluster, Element Count, Estate number, Polar surface area, and alignment-independent topological descriptors have been computed from the chemical structures with an aim of identifying the structure/activity relationship of 1-(4-methoxyphenethyl)-1H-benzimidazole-5-carboxylic acid compounds which would, in turn, predict their biological activity.

The preprocessing of the independent variables (i.e., descriptors) was performed by removing invariable descriptors (constant column), which resulted in a total of 230 descriptors to be used for QSAR analysis. The descriptors having the same or almost the same value, or those that were highly correlated with other descriptors, were removed initially. The various alignment-independent (AI) descriptors [21] were also calculated.

2.5 Cross validation

Models generated by 2D-QSAR studies were cross validated using a standard LOO procedure. The cross validated r² (q²) value was calculated using Eq. (1)(1) $q^2=1-\frac{\sum (y_i-{\hat{y}}_i)}{\sum {(y_i-y_{\text{mean}})}^2}$(1) , where y_i and ŷ_i are the actual and predicted activities of the ith molecule, respectively, and y_mean is the average activity of all of the molecules in the training set [22]. Because the calculation of the pairwise molecular similarities, and hence the predictions, were based upon the current trial solution, the q² obtained is indicative of the predictive power of the current model. The q² was calculated using Eq. (1)(1) $q^2=1-\frac{\sum (y_i-{\hat{y}}_i)}{\sum {(y_i-y_{\text{mean}})}^2}$(1) .(1) $q^2=1-\frac{\sum (y_i-{\hat{y}}_i)}{\sum {(y_i-y_{\text{mean}})}^2}$(1)

The external predictive power of the model is assessed by predicting the pIC₅₀ value of five test set molecules, which are not included in the QSAR model development. The predictive ability of the selected model is also confirmed by pred_r².(2) $\mathrm{Pr}ed\_r^2=1-\frac{\sum (y_i-{\hat{y}}_i)}{\sum {(y_i-y_{\text{mean}})}^2}$(2)

External validation of the generated models was carried out by predicting the activity of the test set of compounds. The predicted r² (pred_r²) value was calculated using Eq. (2)(2) $\mathrm{Pr}ed\_r^2=1-\frac{\sum (y_i-{\hat{y}}_i)}{\sum {(y_i-y_{\text{mean}})}^2}$(2) , where y_i and ŷ_i are the actual and predicted activities of the ith molecule in the test set, respectively, and y_mean is the average activity of all of the molecules in the training set.

3 Results and discussion

For the development of QSAR models for 1-(4-methoxyphenethyl)-1H-benzimidazole-5-carboxylic acid derivatives analogues, PLS methodology was used based on the simulated annealing (K562 and CEM cell) feature selection method. The best 2D-QSAR models were chosen for discussion.

pIC₅₀ (K562-cell line) = +0.9643 (±0.247) StNE-index +0.2487 (±0.0085) SddsN (nitro) count −0.3521 (±0.0483) T_N_O_1 + 0.7453 (±0.1574) SsNH₂-index

N_training = 13, N_test = 5, r² = 0.8477, q² = 0.7584, F test = 47.945, r²_se = 0.3965, q²_se = 0.3485, pred_r² = 0.7855

The model-1 simulated annealing-PLS (SA-PLS) method can explain the coefficient of determination (r²), 84.7%, of the variance in the observed activity values. The cross-validated squared correlation coefficient (q² = 0.7584) was 75%, and the predictivity for the external test set (pred_r² = 0.7855) was approximately 78%. The model is validated by the values α_ran_r² = 0.00369, α_ran_q² = 0.001, α_ran_pred_r² = 0.0001, best_ran_r² = 0.015, best_ran_q² = 0.02471, Z_score_ran_r² = 0.0176 and Z_score_ran_q² = 0.5234. The correlation matrix is shown in Table 2. Model-1 showed that evaluation of the descriptor SddsN (nitro) count suggests that the total number of nitro groups connected with one single and two double bonds at the R position will lead to improved activity. The presence of electron-withdrawing groups in compounds 1, 3, 12 and 14 resulted in highly active antileukaemic activity. The descriptors StNE-index and SsNH₂Count (the total number of –NH₂ groups connected with one single bond) suggest that the presence of cyano group and amino atoms (∼30) result in an increase in the activity at the R position in compounds 2, 4, 13 and 15, respectively. The contributing descriptor T_N_O_1 indicates that the presence of methoxy (compounds 5 and 16) at the R position of the ring would also affect the activity. The contribution chart for the QSAR model is shown in Fig. 1a and plots of the observed vs. predicted values of pIC₅₀ are shown in Fig. 1b.

Table 2 Correlation matrix between descriptors present in the best QSAR model 1.

Parameter	pIC50	StNE-index	SddsN (nitro) count	T_N_O_1	SsNH2-index
pIC50	1.0000
StNE-index	0.5421	1.0000
SddsN (nitro) count	0.3852	0.4731	1.0000
T_N_O_1	0.4732	0.5476	0.7485	1.0000
SsNH2-index	0.3986	0.4866	0.6197	0.8055	1.0000

Fig. 1 (a) Contribution charts of the descriptors for the 2D model-1. (b) Graphs of observed vs. predicted activity of 2D model-1. (c) Contribution charts of the descriptors for the 2D model-2. (d) Graphs of observed vs. predicted activity of 2D model-2. (e) Contribution charts of the descriptors for the 2D model-3. (f) Graphs of observed vs. predicted activity of 2D model-3. (g) Contribution charts of the descriptors for the 2D model-3. (h) Graphs of observed vs. predicted activity of 2D model-4.

pIC₅₀ (K562 cell line) = 0.6588 (±0.0531) SsOHcount +0.4985 (±0.0587) T_N_F_1 − 0.8753 (±0.1498) HOMO energy

N_training = 13, N_test = 5, r² = 0.7953, q² = 0.7336, F test = 38.095, r²_se = 0.3487, q²_se = 0.3127, pred_r² = 0.7069

Equation 2 explains 79% of the variance in activity with low standard error estimation. It also has an internal (q²) and external (pred_r²) predictive ability of ∼73 and ∼70%, respectively. The contribution chart of the selected descriptors is represented in Fig. 1c and plots of the predicted vs. observed values of pIC₅₀ are shown in Fig. 1d. The QSAR model showed a positive contribution of the SsOHcount indices for the total number of -OH groups connected with one single bond and indicated that the activity was increased with the presence of hydroxy groups in 2nd to 4th position in the benzimidazole moiety. T_N_F_1 (∼27%) indicated that the presence of a fluorine atom (i.e., -F) at position site R on the nucleus (as in compounds 1, 2, 9, 12, 13, 18 and 20) would increase the antileukaemic activity. The descriptor of HOMO energy contributes negatively to the activity and substituent hydroxy and methoxy groups on the ring.

pIC₅₀ (CEM cell line) = 0.3436 (±0.0195) smr +0.5985 (±0.0968) SsCH₃Count −0.4874 (±0.0854) Mol. Wt

N_training = 13, N_test = 5, r² = 0.7487, q² = 0.6853, F test = 28.5973, r²_se = 0.3653, q²_se = 0.3218, pred_r² = 0.6542

The same data set subjected to the SA-PLS method resulted in a coefficient of correlation of 0.7487 and an internal predictive power of 68%, with an external predictivity of 65%. The descriptor smr (∼20%) suggests that an increase in the R position may lead to an increase in the activity. The descriptor SsCH₃Count implies that a bulky hydrophobic group at the R position is favourable for the antileukaemic activity of the compound. The descriptor of mol. wt. indicates that the presence of a bulky group at the R position is detrimental for the antileukaemic activity. The contribution chart of the selected descriptors is represented in Fig. 1e and plots of the predicted vs. observed values of pIC₅₀ are shown in Fig. 1f.

pIC_{50 (CEM)} = 0.7453 (±0.0431) SssOE-index −0.5421 (±0.0665) SsBr-index + 0.4313 (±0.0321) T_2_C_2

N_training = 13, N_test = 5, r² = 0.7232, q² = 0.6976, F test = 21.786, r²_se = 0.3198, q²_se = 0.3068, pred_r² = 0.6637

SA-PLS model-4 reveals that the descriptor SsBr-index reveals that the number of bromines connected with one single bond (i.e., -Br) at the R position on the ring is detrimental for the antileukaemic activity (compounds 10 and 21, with substituents at the R positions producing good activity). SssOE-index, describing the number of oxygen atoms connected with two single bonds, indicated that the activity was increased with the presence of electron donating groups, and T_2_C_2 indicates that presence of substituents with direct attachment of carbon on the aromatic ring (i.e., -COOH) is favourable for the antileukaemic activity. The contribution chart of the selected descriptors is represented in Fig. 1g, and plots of the predicted vs. observed values of pIC₅₀ are shown in Fig. 1h. The overall model is validated by predicting the biological activities of the test molecules, as indicated in Table 3.

Table 3 Observed and predicted activities (pIC₅₀) for the training and test set compounds.

Com. No	2D Model-1		2D Model-2		2D Model-3		2D Model-4
	Pred.	Res.	Pred.	Res.	Pred.	Res.	Pred.	Res.
1	5.4747	0.0481	5.4479	0.0749	5.4209	−0.023	5.4463	−0.0484
2	4.1393	−0.0314	4.0430	0.0649	3.7134	0.0435	3.7979	−0.041
3	3.9169	0.0039	3.9131	0.0077	3.8764	0.0372	3.9008	0.0128
4	3.7440	0.0007	3.7412	0.0035	3.7259	0.031	3.7142	0.0427
5	3.8540	0.0124	3.8488	0.0176	3.8556	0.0045	3.8867	−0.0266
6	4.1247	0.0302	4.1374	0.0175	4.0779	0.0895	4.0105	0.1569
7	4.1172	−0.0415	4.0343	0.0414	4.2723	−0.2018	3.8017	0.2688
8	4.3699	−0.1181	4.1914	0.0604	4.2240	0.0356	4.1479	0.1117
9	4.2037	−0.0099	4.2165	−0.0227	4.0568	0.1302	4.2883	−0.1013
10	3.7502	0.1194	3.9491	−0.0795	3.8920	−0.0126	3.8578	0.0216
11	3.7402	0.0045	3.6953	0.0494	3.7283	0.0116	3.7527	−0.0128
12	3.5938	0.0082	3.5019	0.1001	3.6228	−0.012	3.6056	0.0052
13	3.6953	0.0036	3.6937	0.0052	3.7122	−0.0023	3.8830	−0.1731
14	4.4128	−0.0149	4.4057	−0.0078	4.4458	−0.0691	4.4804	−0.1037
15	4.1183	−0.0214	4.1976	−0.1007	4.1501	−0.0796	4.1972	−0.1267
16	3.9375	0.0211	3.9468	0.0118	3.8865	0.0642	3.8716	0.0791
17	3.6218	−0.0021	3.6033	0.0164	3.6274	−0.0148	3.6099	0.0027
18	3.8747	0.0318	3.8879	0.0186	3.9033	0.0175	3.9113	0.0095

Pred.: Predicted; Res.: Residual.

4 Conclusion

The present study describes the successful application of 2D-computational approaches to identify the essential structural requirements in 3D chemical space for the modulation of antileukaemic agents of substituted benzimidazole derivatives. The QSAR models suggest that -OCH₃ and amine groups are beneficial, indicating that the presence of a bulky substituent at these positions should improve biological activity. The results of analyses for the descriptors SsOHcount and SssOE-index suggest that the hydroxy and methoxy properties of the molecules are responsible for the determination of their activities. At the R position, atoms such as H, F, Cl, Br, alkyl groups and aromatic rings were added. The contribution of each descriptor specifies the properties that should be present in the drug lead to enhance its activity. Further QSAR studies, along with these findings, will helpful for the design of 1-(4-methoxyphenethyl)-1H-benzimidazole-5-carboxylic acid derivatives as antileukaemic agents.

Conflict of interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgements

The author wishes to express gratitude to V-life Science Technologies Pvt., Ltd. for providing the trial version software for the study.

Notes

Peer review under responsibility of Taibah University