INTRODUCTION

The growth and division of cancer cells are usually faster than normal cells. Chemotherapy is an effective way to treat tumor cells. However, chemotherapeutic drugs are powerful and principally cause impairment to healthy cells¹ leading to subside in their usage. This predicament has created a medical urge to flourish effective antitumor agents with heightened safety outline.^2,3 As clinically proven cancer targets, histone deacetylase (HDAC) inhibitors have been established as a flourishing tactic for the progress of new anticancer agents.^2,3,4

Acetylation and deacetylation of histone proteins play an essential role in transcription and regulation of genes in eukaryotic organisms. The enzymes viz., histone acetyltransferase and HDACs, play an influential role behind this.^5,6 The imbalance of any of them may result hindrance in differentiation and proliferation of typical cells and conduct start of tumor cells. Overexpressed HDAC effectuates the eviction of acetyl groups from histones, leading to compression of chromatin and downregulation of tumor suppressor genes.^7,8,9,10 Cell cycle arrest, chemosensitization, apoptosis induction, and overexpression of tumor suppressors are some of the primary mechanisms controlled by HDAC inhibitors.¹¹ To date, 18 members are present in the mammalian HDAC family, which are classified into four classes; class I-IV, on the basis of their sequence homology with the yeast protein. Class I enclose 1, 2, 3, and 8 isoforms promoting cellular proliferation and hindering apoptosis. Class II is further classified into class IIa with isoforms 4, 5, 7, and 9 and class IIb consist of 6 and 10. Class I and II forbid the cellular differentiation. Some class II isoforms (HDACs 4, 6, 7, and 10) boost cellular migration and angiogenesis, which are two crucial means for cancer metastasis. Class IV is with a lone member of HDAC 11. Classes I, II and IV act by Zn²⁺ reliant mechanism, while class III shows homology with silent information regulator 2, needing NAD⁺ as cofactor for catalysis.^10,11,12

Five HDAC inhibitors have been approved to date for the treatment of different types of cancers. Vorinostat [suberoylanilide hydroxamic acid (SAHA)] and romidepsin (FK228) have been approved for the treatment of cutaneous T-cell lymphoma, while belinostat (PXD101) and tucidinostat (chidamide) (CS055) have been approved for peripheral T-cell lymphoma. Panobinostat (LBH589) finds application in the treatment of multiple myeloma. Besides this, several HDAC inhibitors  are currently under different phases of clinical trials, i.e., rocilinostat (multiple myeloma), entinostat (breast cancer), and tacedinaline (lung cancer) are in phase I, II, and III, respectively, as represented in Figure 1. HDAC inhibitors are also structurally classified as hydroxamic acids, benzamides, cyclic tetrapeptides, short chain fatty acids, electrophilic ketones, etc.^10,13

Remarkably, because these medications are mostly pan-HDAC inhibitors or target many HDAC isoforms, they have a lot of negative effects. Because of their low toxicity and limited off-target effects, isoform-selective HDAC inhibitors may provide therapeutic benefits. As a result, in recent years, investigations on HDAC inhibitors have focused on isoform- or class-specific inhibition.^14,15

Although all these types bear a resembling core structure comprising of three key components, i.e. (i) zinc binding group (ZBG) responsible for chelation of zinc ion at active site; (ii) cap group (a hydrophobic or aromatic or heteroaromatic moiety) accounts for interaction with residues of HDACs external pocket, and (iii) a linker (with optimal length) accounts for joining the ZBG and cap group. The latter two components, i.e. cap group and linker, are being employed for structural modification to obtain compounds with selective and optimum anticancer activity.¹³

Computer-aided molecular drug design plays an important role in design and discovery of novel chemical entities. The role of computational study of HDAC enzymes is evolving nowadays, with particular emphasis on molecular modeling for the development of HDAC inhibitors with enhanced selectivity and effectiveness. Generally, 3D quantitative structure-activity relationship (QSAR) studies are complemented by docking studies.¹⁶

Various studies have been reported by many scientists for studying and developing HDAC inhibitors using computational tools and techniques.¹⁶ Kim et al.¹⁷ synthesized many δ-lactam-based HDAC inhibitors containing modified cap groups. Hamblett et al.¹⁸, employed MS-275 as the lead moiety and modified it around pyridine ring for designing novel HDAC inhibitors with enhanced class I selectivity. Estiu et al.¹⁹ studied the structural basis for the selectivity of class II-selective HDAC inhibitors SAHA, turacin, and NK308 using molecular dynamics simulations approach. Huhtiniemi et al.²⁰ disclosed a relative modeling of human SIRT1. Xie et al.²¹ reported a QSAR study on HDACi for the identification of structural features responsible for anticancer activity. Chen et al.²² selected around 30 known HDAC inhibitors for designing 3D QSAR pharmacophore model to recognize critical ligand features for HDAC inhibition activity. Ragno et al.²³ accomplished 3D QSAR studies for their newly designed class II selective HDAC inhibitor (APHAs) against maize HD1-A and HD1-B with acceptable selectivity.

The upregulation of the HDAC enzyme has been linked to a variety of cancers, making it a possible therapeutic target. The research found possible inhibitors of human HDAC enzyme by screening several biologically active molecules from several databases.²⁴ Various bioinformatics tools can be used to screen prospective drugs before moving forward with wet-lab research.²⁵ As a result, computer-aided drug design has been shown to be extremely useful in lowering drug development costs and risks, while also increasing speed and accuracy of drug discovery.²⁶ Molecular docking, binding mode and energy, and hydrogen bond interactions aid in the identification of a possible inhibitor in the active site of the HDAC target protein from a dataset. In other words, the ligand-receptor interaction is predicted by molecular docking.²⁷ In addition, QSAR model assesses the biological activity of experimental data by comparing it with chemical descriptors of known training set substances. The application domain and appropriate validation approaches were used to determine the reliability and robustness of the developed QSAR models.²⁸ Nonetheless, the study discusses the creation of an atom-based 3D QSAR model that specifies molecular level comprehension and structure–activity relationship regions for a dataset of chemicals. The created QSAR model considers essential pharmacophoric properties such as average shape, hydrophobic/non-polar areas, electrostatic (positive ionic and negative ionic), electron withdrawing, and electron donating for their respective positive and negative coefficient patterns. Different metrics of QSAR models from partial least squares (PLS) statistical analysis, such as Q², R², standard deviation (SD), stability, F, and root mean squared error values, also show that the model has strong predictive capacity. As a result, the research above aims to provide useful information for designing innovative and effective HDAC inhibitors using computational and bioinformatic approaches.

MATERIALS AND METHODS

The present work comprises computer-aided drug design studies. Thus, it does not require ethics committee approval and patient consent for its accomplishment.

The computational analysis was performed employing the Schrodinger suite (Maestro v 9.3, LLC, New York) including protein prep wizard, ligprep, grid generation, glide XP dock, and 3D QSAR model designing.

Biological dataset

The data resources were collected from the research papers.^29,30,31,32 The literature review clearly exhibited that the heterocyclic linker in hydroxamic acid-based HDAC inhibitors improves activity by facilitating ligand receptor binding. The selected compounds have similar skeletons and biological assay methods. A data set of 57 compounds was chosen for the study along with inhibitory concentration 50 (IC₅₀) values in µM against human carcinoma cancer cell lines as shown in Table 1. IC₅₀ value was used as a dependent variable in the QSAR study. IC₅₀ values of all compounds for different pharmacophore studies were changed into negative logarithm of IC₅₀ (pIC₅₀).³³ These data are critical for constructing good 3D QSAR models for investigating structure-activity relationships.

Protein preparation

The protein preparation wizard in Maestro v. 9.3 was practiced to organize the receptor in order for docking studies.³⁴ The binding region of HDAC inhibitor was initially studied by complexed crystal structure of SAHA (proto type HDAC inhibitor) with HDAC protein (PDB ID: 1ZZ1).³⁵ This task was carried out in three steps, (i) importing the protein from PDB followed by processing to fix its structure, (ii) reviewing chemical correctness of structure and its modification by adding missing hydrogen atoms and neutralizing the remotely situated side chain from binding sites, (iii) refining the orientation of optimized H-bound groups and geometric minimizing the structure by OPLS_2005 force field by facilitating the realignment of hydroxyl groups of side chains.^35,36

Ligand preparation

Lig prep version 2.5 (Schrodinger, LLC, NY) was employed for constructing and processing the selected ligands.³⁷ Initially, the structures of all these ligands were drawn in ChemDraw Professional version 16.0 and saved in mol format. In lig prep, the ligands were picked from their mol files and proceeded through several steps, like generation of 3D structures from their 2D structures, removal of low energy conformers, formation of stereoisomers and ionization state of ligands, addition and elimination of hydrogen atoms and counter ions, respectively, and lastly energy minimization using OPLS_2005 force field. The ligands were geometrically optimized through Optimized Potentials Liquid Simulations 2005 (OPLS_2005) force field.³⁸ The partial atomic charges were figured out employing the OPLS_2005 force field.³⁹

Docking studies

Glide version 5.8 (Schrodinger, LLC, NY) molecular docking tool was used for docking studies.⁴⁰ An effective interaction of hydroxamate derivatives with the target protein (PDB ID: 1ZZ1) to estimate the potential response against tumor cells can be predicted by this study. 1ZZ1, the target protein, was acquired from the protein data bank (PDP)³⁹ and was prepared for the docking task by working on “protein preparation wizard” in Maestro version 9.3. The selected ligands were prepared by lig prep in Maestro version 9.3. Low energy conformers of ligands were screened. The grid was generated on the receptor protein by following receptor generation module in glide and finally the screened ligands were docked into the receptor grid containing protein exercising XP and SP docking approaches.^41,42,43,44

Energetic (e)-pharmacophore hypothesis generation

Both ligand and structure based techniques are combined in energetic (e)-pharmacophore approach. The e-pharmacophore script feature permitting docking post processing option in Maestro v. 9.3 was accomplished for e-pharmacophore hypothesis study.^5,33,45 The module uses energetic tenets of the XP glide scoring function for mapping and creating energy-adjusted pharmacophores, i.e., e-pharmacophores. Thereafter, phase version 3.4 (Schrodinger, LLC, NY) application generates pharmacophore sites using a default set of six chemical attributes: hydrogen bond acceptor (A), hydrogen bond donor (D), hydrophobic group (H), positive ionizable (P), negative ionizable (N), and aromatic ring (R). The glide XP energies of each atom were summed to constitute each pharmacophore site. These sites were then ordered as per their energy and the most affirmative site was picked for pharmacophore generation.

Pharmacophore hypothesis generation

Phase version 3.4 (Schrodinger, LLC, NY) was engaged for the pharmacophore model (hypotheses) generation.⁴⁶ It is a commonly used system-based method for recognizing common pharmacophores and developing 3D QSAR models. Pharmacophore modeling is a ligand-based method for identifying new lead moieties.

The process gets initiated with cleaning of all 57 selected ligands. The conformers of these ligands were created by a macromodel search approach in which maximum number of conformers was 1000 per structure and minimization steps as 100 was set as default. Conformers were minimized using OPLS_2005.⁴⁷ Later, sites were created for all the ligands depending on the values set for activity threshold, which progressively generated common pharmacophore hypothesis (CPHs). The CPHs are based upon the activity threshold of active and inactive molecules. A maximum of six feature sites was present in each hypotheses, namely, hydrogen bond donor (D), hydrogen bond acceptor (A), hydrophobic group (H), positively charged group (P), negatively charged group (N), and aromatic ring (R). These generated hypotheses were monitored on the basis of survival, survival-inactive, and post-hoc scores. The hypotheses possessing lowest relative conformational energy and highest adjusted survival score were selected for building the QSAR model.⁴⁸

3D QSAR model development

The dataset was efficiently segregated into training (38) and test (19) sets for analysis using a random and rational division method. In the phase module, pharmacophore and atom-based alignments are available to orient 3D structures of compounds. In this study, an atom-based QSAR model was used, which explained the better structure-activity relationship. Initially, the overall dataset was segregated into a modeling set (80% compounds) and an external evaluation set (20% compounds) employing random division approach. The modeling set was further sectioned into a training set (comprising of 80% of the modeling set) and a test set (comprising of 20% of the modeling set) again using the rational division method. The best fitting model was generated through a random division. The atom-based QSAR model was generated for ligands by selecting the best fit hypothesis with good scoring value, keeping 1Å as grid spacing and maximum PLS factor as 5. The QSAR results were later visualized, which ultimately helped in the optimization of the thrust structure of dataset.^49,50,51

Validation of the pharmacophore model

The primary goal of QSAR model is to estimate biological activity of novel molecules. Internally and externally, the model developed would sound. A training set and a test set were created from the data. With 38 compounds in the training set, atom-based 3D-QSAR models were created for hypotheses. By estimating the activities of 19 test molecules, top QSAR model was externally validated.

Statistical analysis

Statistical criteria such as squared correlation coefficient (R²), q² (R² for test set), SD of regression, Pearson’s correlation coefficient (pearson’s r), statistical significance (P), and variance ratio were used to internally validate the developed pharmacophore hypotheses (F). The anticipated pIC₅₀ was calculated using the 5^th PLS factor. An increase in the number of PLS factors has no effect on the model’s statistics or prediction ability.

RESULTS AND DISCUSSION

Molecular docking study

The result of docking studies shown in Figure 2 indicated the probable interaction of ligands containing hydroxamate groups with the receptor 1ZZ1. Compound 39 from the dataset exhibits maximum structural alignment with that of SAHA in protein 1ZZ1. 2D interaction diagram of compound 39 docked with 1ZZ1 revealed metallic bond interaction with zinc (Zn2451) of receptor 1ZZ1 and keto group of ligand, hydrogen bond formation between NH- moiety of hydroxamic acid group of the ligand with GLY151 and TYR312 amino acids, and hydrogen bond formation between NH- moiety of hydroxamic acid group of the ligand with GLY151 and TYR312 amino acids. Hydrogen bonding is also obvious between the NH- moiety of indole ring and the ASP98 amino acid. The hydrophobic interaction of PHE208 with five-membered ring B of indole and benzene ring of compound 39 was also found in the diagram.

Energetic (e)-pharmacophore study

An energetic (e)-pharmacophore study displayed its result for compound 45 with protein 1ZZ1. A maximum of seven pharmacophore attributes were taken as default, but five pharmacophore sites were recorded. The created hypothesis presented one hydrogen bond acceptor, two hydrogen bond donors, and two aromatic rings, as presented in Figure 3b. Its ranking order and scoring value as represented in Table 2 clearly indicate that the two aromatic rings R10 and R11 create a hydrophobic environment and acceptor A2 and donor D4 and D6 groups participate in hydrogen bonding.

Pharmacophore generation and 3D QSAR model analysis and validation

Phase version 3.4 presented the outcome for pharmacophore generation and atom-based 3D QSAR modelling. The activity threshold was maintained at a range of 7.6 to 6.9, which divided the dataset into active, moderately active and inactive range. The dataset was further partitioned into training (38 molecules) and test (19 molecules) sets based on structural features and the range of biological activity. The five features containing CPHs were selected based on their high survival score to define the entire binding arena of the molecule as displayed in Figure 3a. The suitable 10 CPHs (Table 3) representing good scores (survival-inactive) were considered for 3D QSAR model design using 5 PLS factors.

CPH AADRR.139 executed best statistical conclusion for PLS factor 5, revealing Q² (0.7142), R² (0.9877), SD (0.1049), F (531.1), P (1.627e-030), root mean square deviation (RMSD) (0.4435), stability (0.4939), and pearson-r (0.8478) (Table 4). The scatter plots of actual vs predicted activity for training and test set compounds were plotted (Figure 4). The results showing comparison of predicted activity with their actual experimental activity was studied and is mentioned in Table 5.

The statistics and predictive ability (q²) of the model did not improve with an increase in the digit of PLS factors. The regression is carried out by creating a series of models with progressively more PLS components. When the number of PLS factors is increased, the model’s accuracy improves until overfitting occurs. Although there is no limit to the number of PLS factors that can be added, in general, adding factors should be halted when the SD of the regression is roughly equivalent to the experimental error. This problem began to appear on models generated after PLS 5. At 5^th PLS factor with the smallest SD of regression, statistical measures like R² and q² were also high (0.9877 and 0.7142, respectively). As a result, for the construction of our atom-based three-dimensional quantitative structure-activity relationship model, 5^th SD of regression component was chosen.

3D QSAR model visualization

There are some essential features in the form of different colored cubes for each feature observed in QSAR visualization maps highlighting an active ligand-receptor interaction. These features indicate the type and position of the attachment of functional groups for showing specific pharmacological activity. They also throw light on toxicity statement of ligands. The ligand 40 from the dataset, more specific from the training set, was carefully chosen as the template molecule for improved understanding of the study. The QSAR model made between hypothesis AADRR.139 and compound 40 is visualized in Figure 5.

The substitution of hydrogen atom of NH- and its adjacent -CH₂ group by a hydrogen bond donating moiety increases the activity; similarly, replacement of hydrogen atom of OH of hydroxamic acid also displays an elevation in activity. Replacement of oxygen of hydroxamic acid by any H-bond donor leads to decrease in activity. The attachment of the hydrophobic group rather than H present at 4^th position of the indoline ring leads to rise in activity. In addition, the substitution of hydrogen atoms of the ethylene moiety present in the linker also escalates the activity. Attachment of the hydrophobic group and an electron withdrawing group in the phenyl ring of linker causes a decline in activity. Substitution of the hydrogen of indoline nitrogen with electron withdrawing moiety increases the activity.

The outcome of these computational studies clearly indicate that among all compounds, higher fitness value and docking score, lesser toxicity, and superior drug properties and more complimentary conformation as compared to the original ligand has been shown by compound 39. Thus, ligand 39 can be considered as a possible lead moiety for the development of newer HDAC inhibitors.

CONCLUSION

HDAC inhibitors, a newer addition to chemotherapy, have been found to play a crucial role in the treatment of cancer. HDAC inhibitors are scarcely available on the market. The computational study performed on hydroxamic acid-based derivatives of the dataset exhibited a convincing outcome.

The various computational studies such as pharmacophore and atom-based 3D-QSAR, molecular docking (XP and SP), and energetic-based pharmacophore mapping effectively established a correlation between the structure of ligands with their predicted biological activity. Both ligand- and structure-based pharmacophore mapping approaches in combination efficiently forecasted this correlation and would be helpful in the design and development of novel HDAC inhibitors as anticancer agents. Moreover, they may also help in designing novel ligands more accurately. The molecular docking study showed maximum structural similarity of compound 39 with that of reference HDAC inhibitor (SAHA). It revealed the crucial intermolecular interactions between the ligand moieties and amino acids in the target protein. The created 3D QSAR pharmacophore model exhibited exceptional regression coefficient standards for the training set, with Q²: 0.7142, R²: 0.9877, and low RMSD: 0.4435. It is expected that the findings of these investigations will be used to develop new structural analogs of substituted phenyl hydroxamide derivatives with anticancer activity.

Ethics

Ethics Committee Approval: Not applicable.

Informed Consent: Not applicable.

Peer-review: Externally peer-reviewed.

Authorship Contributions

Concept: E.S., H.R., Design: E.S., R.V., H.R., Data Collection or Processing: E.S., Analysis or Interpretation: E.S., V.K.P., R.V., H. R., Literature Search: E.S., J.P., Writing: E.S.

Conflict of Interest: No conflict of interest was declared by the authors.

Financial Disclosure: The first author is grateful to DST, New Delhi for providing financial assistance to the author through Inspire Fellowship Program (IF210097). The corresponding author is thankful to the Indian Council of Medical Research (ICMR, New Delhi) for providing financial support in the form of Extra-Mural research project (grant no: 58/33/2020/PHA/BMS).