QSAR and Molecular Docking Studies for Anticancer Drug Targeting BCL-2

Cancer remains one of the leading causes of mortality worldwide, characterized by uncontrolled cell proliferation, genetic instability, and the evasion of programmed cell death. Among the various mechanisms by which cancer cells acquire survival advantages, dysregulation of apoptosis plays a central role. Apoptosis, or programmed cell death, is a tightly regulated process that maintains tissue homeostasis and eliminates damaged or unwanted cells. One of the major regulators of apoptosis is the B-cell lymphoma 2 (BCL-2) protein family, which consists of both pro-apoptotic (e.g., BAX, BAK, BAD) and anti-apoptotic (e.g., BCL-2, BCL-XL, MCL-1) members. The delicate balance between these proteins determines whether a cell undergoes survival or death. The BCL-2 protein, in particular, functions as a potent anti-apoptotic factor by sequestering pro-apoptotic members and preventing the release of cytochrome c from mitochondria, thereby inhibiting the intrinsic apoptotic pathway. Overexpression of BCL-2 has been observed in numerous cancers, including chronic lymphocytic leukemia (CLL), non-Hodgkin’s lymphoma, breast cancer, prostate cancer, and lung cancer, where it contributes to tumor progression, resistance to chemotherapy, and poor prognosis. Given its crucial role in apoptosis evasion, BCL-2 has become an attractive molecular target for anticancer drug discovery. Several small-molecule inhibitors targeting BCL-2 have been developed, among which venetoclax (ABT-199) has achieved clinical success, particularly in hematological malignancies. Venetoclax is a selective BCL-2 inhibitor approved by the FDA, demonstrating that therapeutic inhibition of BCL-2 is feasible and clinically beneficial. However, limitations such as the development of resistance, dose-dependent toxicity, and off-target effects underscore the need to discover novel, more potent, and selective BCL-2 inhibitors. In recent years, computational approaches have gained significant importance in rational drug design. Quantitative Structure–Activity Relationship (QSAR) analysis establishes mathematical models that correlate chemical structures with biological activities, providing insights into the physicochemical and structural requirements for potency. QSAR models allow virtual screening of chemical libraries and prioritization of promising compounds before experimental testing, thereby reducing both cost and time in the drug discovery pipeline. On the other hand, molecular docking is a structure-based computational technique that predicts the preferred binding orientation of small molecules within a target’s active site. Docking studies provide valuable information about binding affinities, key molecular interactions, and structure–activity trends, which are essential for lead optimization. Integrating QSAR and molecular docking offers a powerful dual approach in the search for new BCL-2 inhibitors. While QSAR highlights descriptors and molecular features associated with enhanced activity, docking rationalizes these findings by visualizing ligand–protein interactions within the BH3-binding groove of BCL-2. Together, these techniques provide a complementary framework for the identification and optimization of novel anticancer agents.

dig 1: Role of BCL-2 in Apoptosis Regulation and Cancer Progression

dig 2: The Bcl-2 Family of Apoptotic Regulators

• Bcl-2 protectors protect cells against apoptosis.
•  Bcl-2 killers (eg, Bax and Bak) are proapoptotic proteins that actively kill cells.
• Bcl-2 regulators (widely known as BH3-only proteins) promote cell killing by Either interfering with the protectors or activating the killers.

Acquired resistance to cell death is a common feature of cancer, which involves Abnormal over-expression of pro-survival BCL-2 proteins or abnormal reduction of Pro-apoptotic BCL-2 proteins. These abnormalities lead to the inhibition of apoptosis And are frequently detected in various malignancies. The pro-survival and pro-apoptotic BCL-2 proteins are critical regulators of apoptosis, making them attractive targets for Developing cancer treatment agents. This review discusses the roles of various BCL-2 Family proteins in normal development and organismal function, and how defects in Apoptosis control contribute to the development and therapy resistance of cancer. Finally, the review explores the development of novel BH3-mimetic drugs, inhibitors Of pro-survival BCL-2 proteins, as agents for cancer therapy. Apoptosis is a crucial cellular process that plays a vital role in the survival, development, and functioning of multicellular organisms Deregulation of apoptosis Is linked to various diseases, spanning from cancer to degenerative disorders Two Established pathways to apoptosis are the mitochondrial (intrinsic) pathway, which is Stress-induced and regulated by BCL-2, and the death receptor-induced (extrinsic) Pathway The BCL-2 protein family can be categorized into three sub-groups According to their amino acid sequence similarity and functions: the BH3-only pro- apoptotic proteins (BIM, BID, PUMA, BMF, NOXA, BIK, BAD, HRK), the pro- survival proteins (BCL-2, BCLXL, BCL-W, MCL-1, A1/BFL-1), and the apoptosis Effectors (BAX, BAK, BOK)

Fig 3: The BCL-2 protein family regulates the intrinsic pathway of apoptotic Cell death

The Bcl-2 protein plays a key role in preventing programmed cell death in cancer Cells, making it an important target for anticancer drug development. Researchers have Recently focused on designing BH3 domain mimetics to inhibit Bcl-2, which has led to The approval of Venetoclax (ABT-199) for treating chronic lymphocytic leukaemia. In This study, the researchers extended their previous work on indole-based heterocycles As Bcl-2 inhibitors to investigate quinolin-4-yl based oxadiazole and triazole analogues. The researchers synthesized the target compounds via a common intermediate and Found that some of the quinoline-based oxadiazole analogues showed potent anticancer Activity against Bcl-2-expressing cancer cell lines. Computational molecular modelling Was used to rationalize the Bcl-2 targeted anticancer activity of the most active Analogue (a) and suggested possibilities for designing further potent and selective Bcl- 2 inhibitory heteroaromatics with therapeutic potent

MATERIALS AND METHODS:

Dataset for analysis:

The dataset contain 40 quinazoline based derivative drugs were retrieved from various literature sources with structure elucidated from marine sponge quinazolin derivative drug HEQ-1 was added in the dataset .All ligand chemical structures were designed and converted from 2D structure to 3D structure using Chem Draw software .The dataset has been chosen by which covers the information about its biological activity. The in vitro biological activity data was reported as IC50. The IC50 values were Converted to pIC50. The dataset consists of some highly Active and inactive molecules, with very few molecules in Between. 21 molecules were randomly chosen for training set and 19 molecules were selected for test sets according QSAR Calculations.

Fig 1: The Chemical Structure Of HEQ-1 Ligand.

Fig 2: 3D structure of HEQ-1 ligand molecule

The preparation of ligands was carried out using the LigPrep 2.4 module of the Schrödinger Suite (2010). Initially, the 2D chemical structures of the compounds were imported and converted into accurate 3D geometries. Energy minimization was performed using the MacroModel molecular mechanics force field (MMFF) to ensure realistic representations of the ligands. Since ligands are inherently flexible, a diverse set of thermodynamically accessible conformations was generated. The conformational search was executed using MacroModel’s torsional sampling algorithm, followed by minimization with the MMFF94 force field under a distance-dependent dielectric solvent model. Ionization states at physiological pH (7.0 ± 0.5) were assigned using the Epik 2.1 module, which systematically adds or removes protons based on pKa predictions. For each ligand, a maximum of 100 conformers were generated with an initial 100-step pre-minimization and a 50-step post-minimization. Conformers within a relative energy window of 11.4 kcal/mol (50 kJ/mol) were retained, ensuring structural diversity while discarding high-energy states. Conformers with atomic deviations less than 2.0 Å were filtered for further analysis. Molecular properties such as molecular weight, hydrophobicity (logP), hydrogen bond donors/acceptors, solvent-accessible surface area, and polar surface area were computed using the QikProp 3.3 module to assess drug-likeness and ensure compliance with pharmacokinetic requirements.

Protein Preparation

The BCL-2 protein (PDB ID: 1G5M) crystal structure was retrieved from the Protein Data Bank and prepared using the Protein Preparation Wizard in Schrödinger Suite. Missing hydrogen atoms were added, and protonation states of amino acid residues such as Asp, Glu, Arg, Ser, and His were corrected to reflect physiological conditions. Missing side chains and loops were reconstructed using the Prime module. Energy minimization of the protein was performed with the Impact Refinement module using the OPLS-2005 force field, eliminating steric clashes and optimizing hydrogen-bonding networks. A receptor grid was generated around the BH3 binding groove, defined by the coordinates of the co-crystallized ligand. This grid provided the spatial constraints for docking simulations.

Pharmacophore Site Creation

Pharmacophore modeling was performed in PHASE (Schrödinger Suite). Each ligand was represented as a 3D array of features essential for protein–ligand interactions. The six pharmacophore features considered were:

• Hydrogen bond acceptor (A)
• Hydrogen bond donor (D)
• Hydrophobic group (H)
• Negatively ionizable group (N)
• Positively ionizable group (P)
• Aromatic ring ®

The active analogue approach was employed to identify common pharmacophoric patterns. Ligands were aligned based on their pharmacophore features, and common pharmacophores were selected using a tree-based partitioning algorithm, which clusters pharmacophores based on intersite distances.

 Identification of Common Pharmacophore Hypotheses

From the pool of generated pharmacophores, common hypotheses were identified. The final set of variants included AAHHR, AAHRR, AARRR, AHHRR, AHRRR, and HHRRR. Among these, the best hypotheses (AAHRR, AARRR, and AHRRR) were selected based on scoring metrics. Pharmacophore models were evaluated by their ability to differentiate active from inactive ligands. The alignment quality was assessed using the Root Mean Square Deviation (RMSD ≤ 1.2 Å), while overall performance was determined using the survival score, which combines geometric alignment, activity correlation, and hypothesis selectivity.

3D-QSAR Model Development:

3D-QSAR models were constructed based on the best pharmacophore hypotheses. The dataset was randomly divided into a training set (70%) and a test set (30%). Atom-based 3D-QSAR modeling was employed, as it provides more intuitive insights into structure–activity relationships than pharmacophore alignment alone. In atom-based QSAR, molecules are represented as van der Waals spheres, and each atom is classified into categories:

• D = hydrogen-bond donor
• H = hydrophobic/non-polar atom
• N = negatively ionizable
• P = positively ionizable
• W = electron-withdrawing (hydrogen bond acceptor)
• X = miscellaneous

These classifications were placed on a 3D cubic grid with 1.0 Å spacing, resulting in a binary matrix representation. The data was modeled using Partial Least Squares (PLS) regression, limiting the number of factors to one-third of the training set size. The best model was identified based on high R² (goodness of fit), Q² (predictive ability), and low RMSE values. Validation was performed by predicting the activities of the test set compounds. Robustness was confirmed through randomization tests and external validation.

Molecular Docking Studies

Molecular docking simulations were performed using the GLIDE (Grid-Based Ligand Docking with Energetics) module of the Schrödinger Suite (2010). Receptor grids were defined around the BH3 binding groove of BCL-2, and ligands were docked using Extra Precision (XP) mode to achieve accurate predictions.

The docking workflow consisted of:

1. Grid generation around active site residues.
2. Ligand docking using the GLIDE XP algorithm, which accounts for steric clashes, electrostatics, hydrogen bonding, hydrophobic contacts, and desolvation penalties.
3. Scoring using the GlideScore function

G Score = 0.065 \times E_{\text{vdW}} + 0.130 \times E_{\text{Coulomb}} + \text{H-bonding} + \text{Hydrophobic interactions} + \text{Penalty terms}. Docking poses were further minimized, and post-docking energy refinement was performed using the Prime MM-GBSA method, which estimates binding free energies by combining molecular mechanics with implicit solvation models.

ADMET and Drug-Likeness Evaluation

Docked ligands with favorable binding energies were subjected to ADMET (Absorption, Distribution, Metabolism, Excretion, and Toxicity) analysis using QikProp 3.3, Swiss ADME, and pkCSM. Evaluated parameters included:

• Lipinski’s Rule of Five compliance
• Gastrointestinal absorption and blood–brain barrier permeability
• CYP450 inhibition potential
• Predicted toxicity (mutagenicity, hepatotoxicity, carcinogenicity)

RESULTS AND DISCUSSION

3D-QSAR modeling was performed by dividing the dataset into a training set (21 compounds) and a test set (19 compounds). Pharmacophore analysis generated several hypotheses, among which the five-featured AAHRR hypothesis was selected, comprising one hydrophobic group, three hydrogen bond acceptors, and four aromatic ring features (Figure 3). Scoring and ranking of pharmacophores identified the best ligand hypothesis based on activity (-log10 IC??) and conformational energy. The large difference between scores of active and inactive molecules confirmed the discriminative ability of the model. 3D chemical structure alignment was carried out using the PHASE module (Figure 4). Statistical validation revealed a significant regression model (F = 62.5, P < 0.05), with low standard deviation (SD = 0.285) and RMSE (0.3211), indicating robustness. The cross-validated correlation coefficient (Q² = 0.5147) confirmed the predictive power of the model (Table 2). Visualization of QSAR cubes (Figure 5) showed blue regions corresponding to favorable features enhancing activity, while red regions indicated unfavorable structural contributions. Predicted versus observed activity scatter plots further supported the model reliability (Figure 6). Glide XP docking of the most active compound (HEQ-1) with BCL-2 protein revealed a stable binding interaction, particularly with ARG98 residues, achieving a docking score of –6.12 kcal/mol (Figure 7).

Fig 7: Docking interaction between Bcl-2 protein with HEQ-1 Ligand potential affinity of docking Glide Score (-6.12).