Molecular docking, QSAR and ADMET studies of withanolide analogs against breast cancer

Abstract

Withanolides are a group of pharmacologically active compounds present in most prodigal amounts in roots and leaves of Withania somnifera (Indian ginseng), one of the most important medicinal plants of Indian traditional practice of medicine. Withanolides are steroidal lactones (highly oxygenated C-28 phytochemicals) and have been reported to exhibit immunomodulatory, anticancer and other activities. In the present study, a quantitative structure activity relationship (QSAR) model was developed by a forward stepwise multiple linear regression method to predict the activity of withanolide analogs against human breast cancer. The most effective QSAR model for anticancer activity against the SK-Br-3 cell showed the best correlation with activity (r²=0.93 and rCV² =0.90). Similarly, cross-validation regression coefficient (rCV²=0.85) of the best QSAR model against the MCF7/BUS cells showed a high correlation (r²=0.91). In particular, compounds CID_73621, CID_435144, CID_301751 and CID_3372729 have a marked antiproliferative activity against the MCF7/BUS cells, while 2,3-dihydrowithaferin A-3-beta-O-sulfate, withanolide 5, withanolide A, withaferin A, CID_10413139, CID_11294368, CID_53477765, CID_135887, CID_301751 and CID_3372729 have a high activity against the Sk-Br-3 cells compared to standard drugs 5-fluorouracil (5-FU) and camptothecin. Molecular docking was performed to study the binding conformations and different bonding behaviors, in order to reveal the plausible mechanism of action behind higher accumulation of active withanolide analogs with β-tubulin. The results of the present study may help in the designing of lead compound with improved activity.

Keywords:

• ADMET
• breast cancer
• QSAR
• molecular docking
• withanolides

Introduction

Breast cancer, a heterogeneous group of tumors, is the most commonly diagnosed cancer in women globally. In spite of improvements in early diagnosis and development of several targeted therapeutic methods, breast cancer-related morbidity still remains high. The existing therapeutic approaches are associated with high toxicity, low efficacy, therapeutic resistance and therapy-related morbidity. Recent studies have reported that the natural products owing to their cancer preventive potential will pay a way for more effective, non-endocrine, nontoxic therapeutic approaches for anticancer therapies. Historically, phytomolecules have played a key role in the discovery and development of novel anticancer agents.^1,2 Withania somnifera (commonly known as ashwagandha or Indian winter cherry) is one such medicinal plant, whose all parts are used as ayurvedic remedies for healing various diseases, including inflammation, arthritis, asthma and hypertension.^3–5 The root extract of W. somnifera is known as withanolide, which is composed of 14 bioactive compounds.^6,7 Preclinical experimental data indicates that W. somnifera leaf and root extracts have anticancer potential.^8,9 For example, chemically induced tumorigenesis in the stomach and skin of mice was inhibited significantly following administration of W. somnifera root.⁹ Although no controlled clinical trials of ashwagandha have been reported for any indication, it appears to have a relatively low toxicity profile based on a single human study¹⁰ and toxicological studies in mice.¹¹ In ayurvedic medicine, ashwagandha has also been claimed as an effective agent against arthritis, anxiety, insomnia and stress.¹² W. somnifera is currently regulated in the US and Europe as a dietary supplement. The first withanolide to be characterized was withaferin A (WFA), isolated in 1965,¹³ with its cytotoxic activities being the focus of interest.^14,15 The noncytotoxic, anti-inflammatory¹⁶ and immunomodulatory mechanisms¹⁷ of WFA have however thus far remained rather than poorly considered.

WFA, an important lactone of the withanolide class of phytomolecules, is a highly oxygenated steroidal prototype found in W. somnifera and its related Solanacea species, and exhibits potential anticancer activities.¹⁸ It has been reported to have shown decreases in mammary tumors and pulmonary metastasis in an MMTV-neu transgenic model and is associated with increased apoptosis.¹⁹ WFA-induced apoptosis involves the production of reactive oxygen species²⁰ and induction of FOXO3a and Bim,²¹ while WFA effectively inhibits oncogenic transcription factors such as STAT3,²² resulting in growth inhibition. It has been reported that WFA promotes Notch signaling, plays an oncogenic role and is also often hyperactive in breast cancer cells.²³ Herein, we report the development of a quantitative structure activity relationship (QSAR) model for the antiproliferative activity of withanolide analogs against human breast (SK-Br-3 and MCF7/BUS) cancer cell lines. Moreover, druggability of the compounds was evaluated using Lipinski’s “rule of five” and a series of absorption, distribution, metabolism, excretion and toxicity (ADMET) properties. Molecular docking simulations of these compounds against β-tubulin were also conducted to characterize their binding affinity and interactions.

Materials and methods

Molecular modeling

Molecular modeling studies of withanolide analogs were carried out using Sybyl-X 2.1 (Tripos International, St Louis, MO, USA). The biological activity of all the withanolide analogs measured as inhibitory concentration (GI₅₀) was converted to negative logarithmic concentration in moles (pGI₅₀). Drawing of structures and geometry optimization was performed using ChemBioOffice Suite Ultra v12.0 (2012; Cambridge Soft Corp., Cambridge, UK).^24,25 The Surflex-Dock module in Sybyl was used to construct the binding orientation of withanolide analogs.^26,27

QSAR model development

A QSAR model was developed to screen potential leads against human breast cancer cell lines SK-Br-3 and MCF7/BUS within a training set (Tables S1 and S2).²⁸ A total of 52 physicochemical descriptors were calculated for each compound using Scigress Explorer v.7.7.0.47 (Fujitsu Ltd., Poland; Tables S3). The dataset division into training and test was based on structural/pharmacophore or chemical class similarity in order to include only diverse compounds. Similarly, to select the best subset of physicochemical properties, highly correlated chemical descriptors were excluded through covariance analysis using a correlation matrix (Tables S4 and S5). Finally, a robust QSAR model equation was derived by using multiple linear regression; irrelevant chemical descriptors were removed through a forward stepwise approach leading to a selection of 4 out of 52 descriptors in the final QSAR regression equation. The resulting QSAR model equation showed high regression coefficients (r²) of 0.93 and 0.91 with activity data from the SK-Br-3 and MCF7/BUS cell lines, respectively. The developed QSAR model was validated by predicting the experimental activity of compounds in a test set. Furthermore, the robustness of the developed models was accessed by cross-validation regression coefficient (rCV²) using leave-one-out, test set prediction and randomization methods.^29–31

Molecular docking

Molecular docking simulations were performed to predict the bioactive conformations and explore the binding interactions of withanolide analogs using the Surflex-Dock module of Sybyl X 2.1. The crystal structure of β-tubulin (Protein Data Bank [PDB] code 4IHJ)³² was used as the potential anticancer drug target, with the dataset of compounds as ligands. The Surflex-Dock scoring function was used to assess the strength of the ligand–protein interactions. This scoring function includes energy terms to account for hydrophobicity, polarity, repulsiveness, entropy and solvation. Docking simulations were conducted allowing ligand flexibility but maintaining the protein structure rigid, while all docking parameters were set to default values.

QSAR model validation

A reliable and predictive QSAR model should 1) be statistically significant and robust; 2) provide accurate predictions for external datasets not used during the model development and 3) have its application boundaries defined.

Cross-validation techniques

To explore the reliability of QSAR models, a leave-one-out procedure was performed. The outcome from the cross- validation procedure was cross-validation regression coefficient (rCV²), which was used as a criterion of both robustness and predictive ability of the model. It was calculated using the following equation:

${\text{rCV}}^2=1-\frac{\sum {\left({\mathrm{Y}}_{\text{obs}}-{\mathrm{Y}}_{\text{pred}}\right)}^2}{\sum {\left({\mathrm{Y}}_{\text{obs}}-\overline{\mathrm{Y}}\right)}^2}$

In the above equation, $\overline{\mathrm{Y}}$ means the average activity value of the training dataset, while Y_obs and Y_pred represent the observed and predicted activity values. A high rCV² (>0.5) suggests a reasonably robust model.

Y-randomization test

In this technique, the dependent variable (biological activity) is randomly shuffled, and a new QSAR model is developed using the original independent variable. The new QSAR models (after several trials) are expected to have low r² and rCV² values; if the reverse happens, then an acceptable QSAR model cannot be obtained for the specific modeling method and data. These techniques ensure the robustness of a QSAR model.

Estimation of the predictive ability of a QSAR model

To ensure the high predictive power of a QSAR model, it should be estimated by an external test set of compounds that are not used in building of the QSAR model. The external validation or predictive capacity of the derived model was judged by predictive $r^2(r_{\text{pred}}^2)$ as shown in the following equation:

$r_{\text{pred}}^2=1-\frac{\sum {\left({\mathrm{Y}}_{\text{pred}(\text{test})}-{\mathrm{Y}}_{(\text{test})}\right)}^2}{\sum {\left({\mathrm{Y}}_{(\text{test})}-{\mathrm{Y}}_{(\text{training})}\right)}^2}$

where Y_pred (test) and Y_(test) indicate the predicted and observed activity values, respectively, for test set compounds and ${\overline{\mathrm{Y}}}_{(\mathrm{\text{training}})}$ indicates the average bioactivity of compound in the training set. An acceptable predictive power of a QSAR model $(r_{\text{pred}}^2)$ should be >0.6 for the test set molecules.^33–35

Defining the model – applicability domain (AD)

AD was calculated by the using 23 and 24 inhibitors of MCF7BUS and Sk-Br-3 cell lines used in bioassay model building. It was designed in the form of range values of Lipinski’s parameters, connectivity index (order 0, standard), dipole vector X (debye) and molar refractivity shape index (basic kappa, order 2) in Figure 1 and the connectivity index (order 0, standard), dipole vector X (debye) and shape index (basic kappa, order 2) inhibition in Figure 2. The input parameters for queried compounds were screened through AD decided by the training set.

Figure 1 Plot experimental vs predicted log GI₅₀ with the SK-Br-3 cell line for training set (black color dots) and test set (red color dots) compounds.

Figure 2 Plot of experimental vs predicted log GI₅₀ with the MCF7/BUS cell line for training set (black color dots) and test set (red color dots) compounds.

Prediction of pharmacokinetic (PK) parameters

It would be advantageous in the process of drug discovery, if the ADMET properties of drug molecules were predicted earlier. This information helps the chemist to assess the PK profile of molecules. The PK properties depend on chemical descriptors of drugs, which determine their ADMET. The PK parameters were calculated by QikProp v3.2³⁶ module of Schrödinger Suite 2011. A number of mathematical predictive ADMET models for different PK parameters were available, namely, aqueous solubility, apparent Caco-2 and MDCK permeability, log K_p for skin permeability, blood–brain barrier (BBB; log BB), the volume of distribution and plasma protein binding (shown as log K_hsa for serum protein binding), which were used to quantitatively predict the ADMET properties of the withanolide analogs. Similarly, the drug-likeness (DL) parameters, such as molecular weight (MW) of ≤500 Da, a logP value of ≤5, hydrogen bond donor ≤5 and hydrogen bond acceptor site (N and O atoms) ≤10, and Veber’s parameters, namely, topological polar surface area (TPSA) (≤140 Ų), which assess the passive molecular transport through membranes and number of rotatable bonds (≤10), used for flexibility measurement were calculated.^35,37–39

Results and discussion

QSAR model for antiproliferative activity against the SK-Br-3 cell line

All the withanolide analogs were evaluated for their predictivity through the developed QSAR models. A forward stepwise variable selection strategies on 52 physicochemical properties (descriptors) using 23 known drugs/compounds in the training set (Tables S1) identified four chemical descriptors, connectivity index (order 0, standard), dipole vector X (debye), molar refractivity and shape index (basic kappa, order 2), responsible for the anticancer activity against the SK-Br-3 cell line.

QSAR model equation 1$\begin{array}{l}\text{Predicted}\text{log}{\text{GI}}_{50}(\mu \mathrm{M})=\\ -0.848\times \text{Connectivity}\text{Index}\\ (\mathrm{\text{order}}0,\text{standard})-0.106\times \text{Dipode}\text{Vector}\mathrm{X}\\ (\text{Debye})+0.247\times \text{Molar}\text{Refractivity}-1.117\\ \times \text{Shape}\text{Index}(\text{basic}\text{kappa},\text{order}2)-2.226\\ r^2=0.934;\mathrm{S}=0.159;\mathrm{F}=31.758;P<0.0001;\\ {\text{rCV}}^2=0.904;\text{Spress}=0.218;\text{SDEP}=0.181;\\ r_{\text{pred}}^2=0.868;r_{(\text{random})}^2=0.314;{\text{rCV}}_{(\text{random})}^2=0.008\end{array}$

$\begin{array}{l}\text{Predicted}\text{log}{\text{GI}}_{50}(\mu \mathrm{M})=\\ -0.848\times \text{Connectivity}\text{Index}\\ (\mathrm{\text{order}}0,\text{standard})-0.106\times \text{Dipode}\text{Vector}\mathrm{X}\\ (\text{Debye})+0.247\times \text{Molar}\text{Refractivity}-1.117\\ \times \text{Shape}\text{Index}(\text{basic}\text{kappa},\text{order}2)-2.226\\ r^2=0.934;\mathrm{S}=0.159;\mathrm{F}=31.758;P<0.0001;\\ {\text{rCV}}^2=0.904;\text{Spress}=0.218;\text{SDEP}=0.181;\\ r_{\text{pred}}^2=0.868;r_{(\text{random})}^2=0.314;{\text{rCV}}_{(\text{random})}^2=0.008\end{array}$

Spress is the standard predictive residual sum of squares; SDEP is the standard deviation of error of prediction.

The derived QSAR equation 1$\begin{array}{l}\text{Predicted}\text{log}{\text{GI}}_{50}(\mu \mathrm{M})=\\ -0.848\times \text{Connectivity}\text{Index}\\ (\mathrm{\text{order}}0,\text{standard})-0.106\times \text{Dipode}\text{Vector}\mathrm{X}\\ (\text{Debye})+0.247\times \text{Molar}\text{Refractivity}-1.117\\ \times \text{Shape}\text{Index}(\text{basic}\text{kappa},\text{order}2)-2.226\\ r^2=0.934;\mathrm{S}=0.159;\mathrm{F}=31.758;P<0.0001;\\ {\text{rCV}}^2=0.904;\text{Spress}=0.218;\text{SDEP}=0.181;\\ r_{\text{pred}}^2=0.868;r_{(\text{random})}^2=0.314;{\text{rCV}}_{(\text{random})}^2=0.008\end{array}$ showed a significant relationship between log GI₅₀ (dependent variable) and the chemical descriptors (independent variables). The value of the regression coefficient (r²=0.934) indicates the existence of ~93% correlation between the activity and the chemical descriptors in the training dataset, while the value of the cross-validation regression coefficient (rCV²=0.903) suggests ~90% prediction accuracy of this QSAR model (Figure 1). A Y-randomization of 100 trials with original descriptor producing an average value of 0.314 and 0.008 for r² and rCV², respectively, suggests the robustness and significant accuracy of the developed QSAR model (Table S6). Furthermore, test set predictivity of 0.868 supported the robustness and accuracy of the model. It is evident among the molecular descriptors that molar refractivity is positively correlated with activity, that is, if molar refractivity increases, the biological activity also increases. However, the connectivity index (order 0, standard), dipole vector X (Debye) and shape index (basic kappa, order 2) are inversely correlated showing an increase in the biological activity with the decrease in the magnitude of these descriptors and vice versa. Based on contributing descriptors, the predictivity of 29 withanolide analogs through a QSAR model indicated that 2,3-dihydrowithaferin A-3-beta-O-sulfate, withanolide 5, withanolide A, withaferin A, CID_10413139, CID_11294368, CID_53477765, CID_135887, CID_301751 and CID_3372729 showed significant anticancer in vitro activity compared to the reference drug camptothecin (CPT) against the SK-Br-3 cell line (Table S7 and Figure S1).

QSAR model for antiproliferative activity against MCF7/BUS

The QSAR model for antiproliferative activity against the MCF7/BUS cell line was developed with a training set of 24 drugs/compounds. A total of 52 chemical descriptors were considered during model development (Table S2), and forward stepwise variable search techniques resulted in four parametric models that included all atom count, dielectric energy (kcal/mole), total energy (hartree) and heat of formation (kcal/mole) correlating significantly against the MCF7/BUS cell line.

QSAR model equation 2$\begin{array}{l}\text{Predicted}\text{log}{\text{GI}}_{50}(\mu \mathrm{M})=\\ -0.235\times \text{Atom Count}(\text{all atoms})+1.843\\ \times \text{Dielectric Energy}(\text{kcal}/\text{mole})+0.080\\ \times \text{Total Energy}(\text{Hartree})-0.005\\ \times \text{Heat of Formation}(\text{kcal}/\text{mole})+5.588\\ r^2=0.915;\mathrm{S}=0.176;\mathrm{F}=34.757;P<0.0001;\\ {\text{rCV}}^2=0.856;\text{Spress}=0.245;\\ \text{SDEP}=0.214;r_{\text{pred}}^2=0.632;\\ r_{(\text{random})}^2=0.236;{\text{rCV}}_{(\text{random})}^2=0.004\end{array}$

$\begin{array}{l}\text{Predicted}\text{log}{\text{GI}}_{50}(\mu \mathrm{M})=\\ -0.235\times \text{Atom Count}(\text{all atoms})+1.843\\ \times \text{Dielectric Energy}(\text{kcal}/\text{mole})+0.080\\ \times \text{Total Energy}(\text{Hartree})-0.005\\ \times \text{Heat of Formation}(\text{kcal}/\text{mole})+5.588\\ r^2=0.915;\mathrm{S}=0.176;\mathrm{F}=34.757;P<0.0001;\\ {\text{rCV}}^2=0.856;\text{Spress}=0.245;\\ \text{SDEP}=0.214;r_{\text{pred}}^2=0.632;\\ r_{(\text{random})}^2=0.236;{\text{rCV}}_{(\text{random})}^2=0.004\end{array}$

This QSAR equation shows the relationship between in vitro experimental activities (log GI₅₀) as the dependent variable and the four chemical descriptors mentioned earlier as independent variables. The value of the regression coefficient (r²=0.91) shows the existence of ~91% correlation between the activities and the chemical descriptors of the training dataset compounds, while the value of cross-validation regression coefficient (rCV²=0.85) points to the 85% prediction accuracy of the QSAR model (Figure 2). A Y-randomization of 100 trials with the original descriptor producing an average value of 0.236 and 0.004 for r² and rCV², respectively, suggest the robustness of the developed QSAR model (Table S6). Furthermore, test set predictivity of 0.632 supported the robustness of the model. In the QSAR model, the molecular descriptors, namely, atom count (all atoms), dielectric energy (kcal/mole) and total energy (hartree), correlate positively with activity, that is, if any of these variables increase, the biological activity against the breast cancer also increases. Moreover, the heat of formation (kcal/mol) is negatively correlated with activity, showing that the biological activity decreases if the value of this descriptor increases. Analyzing the most active analogs, CID_73621 and CID_435144 possess higher anti-proliferative activity than the reference drugs 5-fluorouracil (5-FU) and CPT (Table S8 and Figure S1).

Binding affinity and interactions with β-tubulin

Molecular docking simulations were carried out to elucidate the likely binding affinity and binding interactions of withanolide analogs with anti-cancer target β-tubulin receptor protein. The predictions of the docking simulations are summarized in Table 1. Withanolide analogs were found to bind to the same active site as found for X-ray crystal structure. The reference compounds 5-FU and CPT were predicted to have significant binding affinity to β-tubulin with total docking scores of 2.5304 and 4.1837, respectively, while withanolide active analogs CID_11294368 and CID_53477765 showed higher binding affinities with total docking scores of 6.4885 and 6.1862, respectively. The binding interactions of these withanolide analogs within the active site of β-tubulin receptor protein are shown in Figure 3. By contrast, 2,3-dihydrowithaferin A-3-beta-O-sulfate, withanolide 5, withanolide A, withaferin A, CID_73621, CID_435144, CID_10413139, CID_135887 and CID_301751 were predicted to have moderate binding affinities to β-tubulin with total docking scores of 5.6272, 5.958, 5.6691, 5.4153, 5.3516, 5.9275, 5.0615, 5.2200 and 5.2417, respectively. The withanolide analog CID_3372729 was predicted to have low binding affinity with β-tubulin showing a total docking score of 4.6965. The chemical nature of binding amino acid residues in β-tubulin that were predicted to interact with the various withanolide analogs was hydrophobic (eg, Ile-165, Pro-162, Met-166 and Val-257), polar acidic (eg, Asp-163 and Asp-199), polar basic (eg, Arg-158, Arg-164 and Arg-253), polar uncharged (eg, Thr-198) and polar amide (eg, Asn-197). Consequently, the interactions of all analogs with β-tubulin involved predominantly strong hydrophobic interactions, leading to high affinities and likely high anticancer activity. The active analogs also formed one to three H-bonds with β-tubulin, contributing to the strength of their interactions, which may lead to a high inhibitory activity of withanolide analogs with β-tubulin.

Table 1 Comparison of binding affinities and interactions of standard anticancer drugs and active withanolide analogs for breast cancer receptor (PDB code 4IHJ)

Compound name	Total score	Amino acids in active pocket within 4 Å	Interacting amino acids	H-bond length (Å)	No of H-bonds
2,3-Dihydrowithaferin A-3-beta-O-sulfate	5.6272	Arg-158, Pro-162, Asp-163, Met-166, Glu-196, Asn-197, Thr-198, Asp-199, Arg-253, Val-257	Arg-158	2.1	1
12-Deoxy-withastramonolide	5.3729	Arg-158, Pro-162, Arg-164, Met-166, Glu-196, Asn-197, Thr-198, Asp-199, Arg-253, Val-257	Arg-158	2.1	1
Withanolide 1	4.0823	Arg-158, Pro-162, Asp-163, Arg-164, Met-166, Asn-197, Thr-198, Asp-199, Arg-253, Val-257	Arg-158 Arg-164	2.1 1.8	2
Withanolide 2	4.6426	Arg-158, Pro-162, Asp-163, Arg-164, Ile-165, Met-166, Asn-197, Thr-198, Asp-199, Arg-253, Val-257	Asn-197	2.0	1
Withanolide 3	4.6662	Arg-158, Pro-162, Met-166, Glu-196, Asn-197, Thr-198, Asp-199, Arg-253, Val-257, His-266	Pro-162 Arg-158 Arg-253	1.8 1.8 1.9	3
Withanolide 4	5.0296	Arg-158, Pro-162, Asp-163, Arg-164, Ile-165, Met-166, Asn-197, Thr-198, Asp-199, Arg-253, Val-257	Arg-158	2.0	1
Withanolide 5	5.958	Arg-158, Pro-162, Met-166, Val-195, Glu-196, Asn-197, Thr-198, Asp-199, Arg-264, His-266	Pro-162	2.0	1
Withanolide A	5.6691	Arg-158, Pro-162, Asp-163, Arg-164, Ile-165, Met-166, Asn-197, Thr-198, Asp-199, Arg-253, Val-257	Asp-199	1.7	1
Withaferin A	5.4153	Arg-158, Pro-162, Asp-163, Arg-164, Ile-165, Met-166, Asn-197, Thr-198, Asp-199, Arg-253, Val-257	–	–	–
Withanolide D	6.5523	Arg-158, Pro-162, Asp-163, Met-166, Asn-197, Thr-198, Asp-199, Arg-253, Val-257	Pro-162	1.8	1
CID_73621	5.3516	Arg-158, Pro-162, Asp-163, Arg-164, Ile-165, Met-166, Asn-197, Thr-198, Asp-199, Arg-253, Val-257	Asp-199 Arg-158 Pro-162	1.7 2.1 1.9	3
CID_135887	6.2234	Arg-158, Pro-162, Asp-163, Arg-164, Ile-165, Met-166, Asn-197, Thr-198, Asp-199, Arg-253, Val-257	Asp-199	2.0 1.7	2
CID_301754	4.4921	Arg-158, Val-195, Glu-196, Asn-197, Thr-198, Asp-199, Arg-253, Val-257, His-266	Val-195 Arg-158	1.8 1.9	2
CID_435144	5.9275	Arg-158, Pro-162, Asp-163, Arg-164, Ile-165, Met-166, Asn-197, Thr-198, Asp-199, Arg-253, Val-257	Asp-199	2.1 1.9	2
CID_3034071	5.8237	Arg-158, Pro-162, Asp-163, Arg-164, Ile-165, Met-166, Asn-197, Thr-198, Asp-199, Arg-253, Val-257	Arg-253 Asp-199	2.1 2.0	2
CID_5315323	5.5121	Arg-158, Pro-162, Arg-164, Met-166, Asn-197, Thr-198, Asp-199, Arg-253, Val-257	–	–	–
CID_10161347	4.8312	Arg-158, Pro-162, Asp-163, Ile-165, Met-166, Asn-197, Thr-198, Asp-199, Arg-253, Val-257	–	–	–
CID_10413139	5.0615	Arg-158, Pro-162, Asp-163, Arg-164, Ile-165, Met-166, Asn-197, Asp-199, Val-257	Asp-199	1.8 2.0	2
CID_11294368	6.4885	Arg-158, Pro-162, Asp-163, Arg-164, Ile-165, Met-166, Asn-197, Thr-198, Asp-199, Arg-253, Val-257	Asp-199 Pro-162	2.1 1.6	2
CID_53477765	6.1862	Arg-158, Pro-162, Asp-163, Ile-165, Met-166, Asn-197, Thr-198, Asp-199, Arg-253, Val-257	Asp-199	1.9	1
CID_161671	4.6248	Arg-158, Pro-162, Asp-163, Arg-164, Ile-165, Met-166, Glu-196, Asn-197, Thr-198, Asp-199, Arg-253, Val-257	Pro-162	1.8	1
CID_135887	5.2200	Arg-158, Pro-162, Asp-163, Arg-164, Ile-165, Met-166, Asn-197, Thr-198, Asp-199, Arg-253, Val-257	Asp-199	2.1	1
CID_161671	6.3312	Arg-158, Pro-162, Asp-163, Arg-164, Ile-165, Met-166, Asn-197, Thr-198, Asp-199, Arg-253, Val-257	Pro-162	2.0	1
CID_301751	5.2417	Arg-158, Ile-165, Met-166, Asn-197, Thr-198, Asp-199, Val-257	Asp-199	1.8	1
CID_3372729	4.6965	Arg-158, Pro-162, Asp-163, Arg-164, Ile-165, Met-166, Glu-196, Asn-197, Thr-198, Asp-199, Arg-253, Val-257	Asp-199	1.9	1
CID_301751	5.6487	Arg-158, Pro-162, Asp-163, Arg-164, Ile-165, Met-166, Asn-197, Thr-198, Asp-199, Arg-253, Val-257	Arg-253	1.9	1
SID_50526634	4.576	Arg-158, Pro-162, Asp-163, Ile-165, Met-166, Asn-197, Thr-198, Asp-199, Arg-253, Val-257	–	–	–
CID_11294368	5.6777	Arg-158, Pro-162, Asp-163, Arg-164, Ile-165, Met-166, Asn-197, Thr-198, Asp-199, Arg-253, Val-257	Asp-199	1.8	–
CID_161671	6.9272	Arg-158, Pro-162, Asp-163, Arg-164, Met-166, Asn-197, Thr-198, Asp-199, Arg-253, Val-257	Pro-162 Arg-158	1.7 2.1	–
5-FU	2.5304	Ile-154, Arg-158, Tyr-161, Pro-162, Arg-164, Met-166, Asn-197	–	–	–
CPT	4.1837	Ile-154, Ile-157, Arg-158, Tyr-161, Pro-162, Ile-165, Met-166, Asn-197, Thr-198, Asp-199, Arg-253	–	–	–

Abbreviations: PDB, Protein Data Bank; 5-FU, 5-fluorouracil; CPT, camptothecin.

Figure 3 (A) Predicted interactions of CID_11294368 with anticancer target enzyme β-tubulin (PDB code 4IHJ) with a docking total score of 6.4885, revealing two H-bonds of length 2.1 and 1.6 Å, respectively, to the binding site pocket residues ASP-199 and PRO-162. (B) Predicted interactions of ID_53477765 with anticancer target enzyme β-tubulin (PDB code 4IHJ) with a docking total score of 6.1862, revealing a H-bond of length 1.9Å to the binding site pocket residue ASP-199.

Abbreviation: PDB, Protein Data Bank.

The most active compound, CID_301751, was predicted to have a significant binding affinity as evidenced by a total docking score of 5.2417 and form a hydrogen bond of length 1.8 Å to the polar acidic residue Asp-199. The amino acid residues binding within a radius of 3 Å with ligand were polar acidic residue Asp-199; polar uncharged residue Thr-198; nonpolar hydrophobic residues Ile-165, Met-166 and Val-257; polar amide Asn-197 and polar basic Arg-158. The presence of these strong hydrophobic and polar interactions explains the high affinity of this molecule and its anticancer activity (Figure 4A).

Figure 4 (A) Predicted interactions of CID_301751 with anticancer target enzyme β-tubulin (PDB code 4IHJ) with a docking total score of 5.2417, revealing a H-bond of length 1.8 Å to the binding site pocket residue ASP-199. (B) Predicted interactions of CID_3372729 with the anticancer target enzyme β-tubulin (PDB code 4IHJ) with a docking total score of 4.6965, revealing a H-bond of length 1.9 Å to the binding site pocket residue ASP-199.

Abbreviation: PDB, Protein Data Bank.

Docking of compound CID_3372729 against β-tubulin resulted in a predicted total docking score of 4.6965 and form a hydrogen bond of length 1.9 Å to the polar acidic residue Asp-199. The binding site amino acid residues within a radius of 3 Å with ligand were hydrophobic Ile-165, Pro-162, Met-166 and Val-257; polar uncharged Thr-198; polar basic Arg-158, Arg-164 and Arg-253; polar acidic Asp-163, Asp-199 and Glu-196 and polar amide Asn-197. The presence of these strong hydrophobic and polar interactions with β-tubulin explains the high affinity of this molecule and its anticancer activity (Figure 4B).

PK properties, DL and bioavailability

All the withanolide analogs were predicted for their ADME properties. The predicted properties have shown that all withanolide analogs satisfy the Lipinski’s rule of five for oral bioavailability. The logP values were found to be in the range (logP<5) for all the compounds; however, withanolide 3, CID_73621 and CID_435144 showed low membrane permeability due to somewhat high MW (Table S9). The high value of their MWs in fact slightly violates one of Lipinski’s rule of five, since it affects the drug excretion (elimination from the body). Moreover, previous studies have suggested that molecules with intermediate lipophilicity are more likely to arrive to the receptor site.^37,38 Predictions related to MDCK cell, skin permeability (K_p), Caco-2, CNS activity, metabolic reactions, log GI₅₀ for HERG K+ channel blockage and human oral absorption in the gastrointestinal tract showed that these parameters for the active analogs fall within the standard ranges normally observed for drugs (Table 2). On the other hand, BBB and intestinal absorption penetration were predicted by developing an ADME model using 2D descriptors polar surface area and ALogP98 that define 95% and 99% confidence ellipses. These ellipses define the regions where well-absorbed compounds are expected to be found. All withanolide analogs were predicted to possess ~95% confidence levels for human intestinal absorption and BBB penetration, except active compounds CID_73621, CID_135887 and CID_435144, which showed 99% confidence level for intestinal absorption and 95% confidence level for BBB penetration, indicating good intestinal absorption and BBB penetration ability. The values and plot of polar surface area and ALogP for the withanolide analogs are shown in Table S10 and Figure 5, respectively.

Table 2 Predicted PK (ADME) properties of active withanolide analogs

Compound name	log Khsa serum protein binding	No of primary metabolites	Predicted CNS activity	Log IC 50 for HERG K+ channel blockage	Apparent Caco-2 permeability (nm/sec)	Apparent MDCK permeability (nm/sec)	log Kp for skin permeability	Human oral absorption in GI tract (%) (±20%)	Qualitative model for human oral absorption
2,3-Dihydrowithaferin A-3-beta-O-sulfate	0.336	4	−2	−4.466	210.574	91.84	−4.012	86.092	High
12-Deoxy-withastramonolide	0.51	5	−2	−4.56	263.97	117.252	−3.825	90.335	High
Withanolide 1	0.024	5	−2	−4.227	305.469	137.297	−3.864	72.114	High
Withanolide 2	−0.023	6	−2	−4.541	225.973	99.12	−4.116	81.681	High
Withanolide 3	0.305	6	−2	−4.462	141.57	59.793	−4.417	69.065	High
Withanolide 4	0.675	5	−1	−4.37	864.965	422.915	−3.007	100	High
Withanolide 5	0.48	5	−2	−4.345	285.542	127.642	−3.749	90.626	High
Withanolide A	0.541	6	−1	−4.447	610.905	290.411	−3.209	100	High
Withaferin A	0.355	4	−2	−4.491	241.508	106.506	−3.896	87.681	High
Withanolide D	1.314	6	0	−4.743	2,770.276	1,488.213	−1.735	100	Low
CID_73621	0.114	7	−2	−4.33	273.565	121.865	−3.693	70.883	Medium
CID_135887	0.441	10	−2	−4.327	167.469	71.699	−4.263	83.724	Low
CID_301754	0.382	5	−2	−4.417	451.939	209.668	−3.472	92.408	High
CID_435144	0.417	7	−2	−4.453	214.202	93.551	−3.899	73.332	Medium
CID_3034071	0.616	6	−2	−4.803	442.997	205.188	−3.465	95.507	High
CID_5315323	0.566	4	−1	−4.144	907.445	445.409	−2.967	100	High
CID_10161347	0.571	4	−1	−4.282	648.591	309.822	−3.251	100	High
CID_10413139	0.419	7	−2	−4.47	463.395	215.418	−3.331	93.191	High
CID_11294368	0.543	6	−1	−4.299	443.996	205.688	−3.478	94.259	High
CID_53477765	0.415	5	−2	−4.464	418.697	193.049	−3.524	92.181	High
CID_161671	0.4	5	−2	−4.475	423.394	195.392	−3.515	92.07	High
CID_135887	0.44	10	−2	−4.343	168.358	72.111	−4.259	83.774	Low
CID_161671	0.921	5	−1	−4.441	1,081.735	538.556	−2.839	100	Low
CID_301751	0.42	7	−2	−4.467	443.979	205.679	−3.367	92.789	High
CID_3372729	0.386	7	−2	−4.448	502.94	235.355	−3.263	93.487	High
CID_301751	0.637	6	−1	−4.332	682.735	327.488	−3.104	100	High
SID_50526634	0.584	4	−1	−4.167	819.929	399.165	−3.066	100	High
CID_11294368	0.549	6	−1	−4.409	521.407	244.71	−3.343	95.931	High
CID_161671	0.399	5	−1	−4.272	474.508	221.008	−3.44	93.029	High
5-FU
CPT	−0.157	3	−1	−5.144	527.805	247.957	−2.931	85.866	High
Standard range	−1.5/1.5	1.0/8.0	−2 (inactive), +2 (active)	Of concern <−5	<25 poor, >500 great	<25 poor, >500 great	−8.0 to −1.0, Kp in cm/h	<25% is poor	>80% is high

Abbreviations: PK, pharmacokinetic; ADME, absorption, distribution, metabolism; GI, gastrointestinal; 5-FU, 5-fluorouracil; CPT, camptothecin.

Figure 5 Plot of PSA versus ALogP for withanolide analogs showing the 95% and 99% confidence limit ellipses for BBB and intestinal absorption, respectively.

Abbreviations: PSA, polar surface area; BBB, blood–brain barrier; ADMET, absorption, distribution, metabolism, excretion and toxicity.

Toxicity risk assessment

The US Food and Drug Administration standard toxicity risk predictor software Osiris Toxicity Properties Calculator locates fragments within the structure of a molecule that indicates a potential toxicity risk.^33–35 Toxicity risks parameters such as mutagenicity (MUT), tumorigenicity (TUMO), irritation (IRRI) and reproductive or developmental toxicity were computed for all the withanolide analogs. All compounds possess no risk of TUMO and developmental toxicity, with the exception of withanolide 2, which possesses high-risk developmental or reproductive toxicity. All active withanolide analogs were predicted to have high-risk MUT, except for compounds 2,3-dihydrowithaferin A-3-beta-O-sulfate, 12-deoxy-withastramonolide, withanolide 3, withanolide 5, withaferin A and CID_135887, which were predicted to have medium-risk toxicity. All active analogs were predicted to have high-risk skin IRRI, except for compounds 2,3-dihydrowithaferin A-3-beta-O-sulfate, 12-deoxy-withastramonolide, withanolide 1, withanolide 2, withanolide 3, withanolide 5, withaferin A and CID_135887, which showed medium-risk toxicity potential at high doses or long-term therapeutic use in human beings. The overall drug score (DS), which combines DL, hydrophilicity and toxicity risk parameters, was calculated for all compounds and was within acceptable limits (Table 3). The overall DS for all the active compounds predicted to be moderate to good, compared to reference anticancer compounds 5-FU and CPT.

Table 3 Predicted toxicity risk parameters (MUT, TUMO, IRRI and REP) of active withanolide analogs

Compound	Toxicity risk parameters				Drug-likeness parameters (Osiris)
	MUT	TUMO	IRRI	REP	CLP	S	DL	DS
2,3-Dihydrowithaferin A-3-beta-O-sulfate	Medium risk	No risk	Medium risk	No risk	−0.06	−3.25	1.93	0.36
12-Deoxy-withastramonolide	Medium risk	No risk	Medium risk	No risk	2.55	−4.47	1.24	0.35
Withanolide 1	High risk	No risk	Medium risk	No risk	2.27	−4.43	−1.81	0.16
Withanolide 2	High risk	High risk	Medium risk	No risk	0.92	−3.9	−2.85	0.15
Withanolide 3	Medium risk	No risk	Medium risk	No risk	2.49	−4.4	1.54	0.34
Withanolide 4	High risk	No risk	High risk	No risk	3.5	−4.98	−1.04	0.13
Withanolide 5	Medium risk	No risk	Medium risk	No risk	2.55	−4.47	124	0.35
Withanolide A	High risk	No risk	High risk	No risk	2.3	−4.53	−0.63	0.15
Withaferin A	Medium risk	No risk	Medium risk	No risk	2.55	−4.47	1.69	0.37
Withanolide D	High risk	No risk	High risk	No risk	2.13	−4.48	−0.43	0.16
CID_73621	High risk	No risk	High risk	No risk	0.93	−3.63	0.97	0.21
CID_135887	Medium risk	No risk	High risk	No risk	2.18	−3.96	1.31	0.21
CID_301754	High risk	No risk	High risk	No risk	2.3	−4.53	0.14	0.17
CID_435144	High risk	No risk	High risk	No risk	1.6	−3.76	0.83	0.2
CID_3034071	High risk	No risk	High risk	No risk	2.78	−4.58	−0.19	0.15
CID_5315323	High risk	No risk	High risk	No risk	3.5	−4.98	−0.6	0.14
CID_10161347	High risk	No risk	High risk	No risk	3.5	−4.98	−0.6	0.14
CID_10413139	High risk	No risk	High risk	No risk	2.44	−4.48	−0.41	0.16
CID_11294368	High risk	No risk	High risk	No risk	2.3	−4.53	−0.63	0.15
CID_53477765	High risk	No risk	High risk	No risk	2.02	−4.35	0.06	0.18
CID_161671	High risk	No risk	High risk	No risk	2.3	−4.53	0.14	0.17
CID_135887	High risk	No risk	High risk	No risk	2.71	−4.3	0.7	0.19
CID_161671	High risk	No risk	High risk	No risk	3.6	−5.26	1.65	0.17
CID_301751	High risk	No risk	High risk	No risk	2.44	−4.58	−0.49	0.15
CID_3372729	High risk	No risk	High risk	No risk	1.76	−4.03	−0.41	0.16
CID_301751	High risk	No risk	High risk	No risk	2.44	−4.49	0.55	0.15
SID_50526634	High risk	No risk	High risk	No risk	3.5	−4.96	−0.6	0.14
CID_11294368	High risk	No risk	High risk	No risk	2.98	−4.74	−0.44	0.14
CID_161671	High risk	No risk	High risk	No risk	2.3	−4.53	0.14	0.17
CPT	High risk	High risk	No risk	Medium risk	1.48	−2.74	5.35	0.25

Abbreviations: MUT, mutagenicity; TUMO, tumorigenicity; IRRI, irritation; REP, reproduction; CLP, CLogP; S, solubility; DL, drug-likeness; DS, drug score.

Conclusion

There are numerous statistically proven examples, where combined approaches of QSAR and molecular docking-based prediction have been applied successfully in the field of drug design and discovery. The present work of QSAR and molecular docking-based prediction of withanolide analogs showed that 2,3-dihydrowithaferin A-3-beta-O-sulfate, withanolide 5, withanolide A, withaferin A, CID_10413139, CID_11294368, CID_53477765, CID_135887, CID_301751 and CID_3372729 against Sk-Br-3 and CID_73621, CID_435144, CID_301751 and CID_3372729 possess a significant anticancer activity against the MCF7/BUS. The QSAR results for SK-Br-3 suggested that connectivity index (order 0, standard), dipole vector X (Debye), molar refractivity, shape index (basic kappa, order 2), whereas for MCF7/BUS, atom count (all atoms), dielectric energy (kcal/mole), total energy (hartree) and heat of formation (kcal/mole) correlated well with the activity. In docking studies, active withanolide analogs showed high binding affinity against β-tubulin receptor protein. The withanolide analogs CID_301751 and CID_3372729 showed good predicted activity and binding affinity to β-tubulin receptor protein. The docking results showed that the major influencing factors of molecular interactions between withanolide analogs and β-tubulin were H-bonds and hydrophobic and electrostatic interactions. The in silico prediction of oral bioavailability (rule of five) and ADMET risk profiling were within their acceptable limit for active analogs. These compounds have rationalized the structural requirement and need further lead optimization for designing of novel β-tubulin inhibitors.

Acknowledgments

We acknowledge the Council of Scientific and Industrial Research (CSIR), New Delhi, India for financial support via the networking project BSC0121 at the CSIR-Central Institute of Medicinal and Aromatic Plants, Lucknow. We are thankful to the Director of CSIR-CIMAP, Lucknow for rendering essential research facilities and support. Author DKY is thankful to the Indian Council of Medical Research (ICMR), New Delhi for the award of Senior Research Fellowship (SRF) and Research Associate (RA). This study was supported by the Basic Science Research Program of the National Research Foundation of Korea (NRF), which was funded by the Ministry of Education, Science and Technology (No. 2012R1A6A3A04038302 and 2017R1C1B2003380).

Disclosure

The authors report no conflicts of interest in this work.