Translate this page into:
In silico identification of 2-oxo-1,3-thiazolidine derivatives as novel inhibitor candidate of class II histone deacetylase (HDAC) in cervical cancer treatment
⁎Corresponding author. usman@ui.ac.id (Usman S.F. Tambunan),
-
Received: ,
Accepted: ,
This article was originally published by Elsevier and was migrated to Scientific Scholar after the change of Publisher.
Peer review under responsibility of King Saud University.
Abstract
Cervical cancer ranks third among the most prevalent deadly cancer in women worldwide and ranks first in developing countries. It is caused by human papilloma virus (HPV) infection. Thus HDACs have become prominent inhibition target for cervical cancer treatment. In order to discover the new alternative HDAC inhibitors (HDACIs), we conducted a computer-aided drug discovery and development (CADDD) based on de novo approach. The compound library is based on 4-[(2-oxo-1,3-thiazolidin-3-yl)carbonyl] aniline. Screening of the best drug leads was evaluated from several parameters, such as docking and interaction analysis, pharmacology, in silico preclinical trial and molecular dynamics analysis. The inhibitory activity of these new designed ligands against Homo sapiens class II HDAC was determined by molecular docking simulation. Docking analysis yielded eight best ligands which have better binding affinity than the standards. Therefore, interaction analysis indicated that all best ligands performed coordination with Zn2+ cofactor in HDAC charge-relay system which are essential for HDAC inhibitory activities of these inhibitors. Pharmacology analysis and preclinical trials of these compounds including pharmacology properties, bioactivity, mutagenicity–carcinogenicity, absorption, distribution, metabolism, excretion and toxicity (ADMET) properties were done through in silico methods. Through this analysis, the best ligands meet Lipinski’s rule of five, have a better drug score than standards, and show good bioactivities, oral bioavailability and ADMET properties. All best ligands also have good synthetic accessibility and were proved to be new compounds that have never been synthesized before. Stability of HDAC–ligand complexes was also calculated through molecular dynamics (MD) simulation. Based on this simulation, complex of the best ligands with corresponding HDAC has a good stability based on RMSD (root mean square deviation) and interaction analysis. The study thus reveals eight best ligands (F, Ib14, O38, Kb17, Gd40, Aa50, Gc42 and Bb38) which have better binding affinity against human class II histone deacetylase (HDAC) through molecular docking, dynamics and interaction analysis. The best ligands were also found to have good bioactivities, oral bioavailability and ADMET properties through in silico pharmacology analysis and preclinical trial. These compounds were found to have a good synthetic accessibility; therefore they could be synthesized for further biological and clinical test.
Keywords
HPV
Cervical cancer
2-Oxo-1,3-thiazolidine
HDAC
Inhibitors
CADDD
1 Introduction
Cervical cancer ranks third as the most common deadly cancer in women worldwide and ranks first in the developing countries (Martin and O’Leary, 2011; Zagouri et al., 2012). In 2008, it is estimated that 529,000 cervical cancer cases were diagnosed and 274,000 of them died (Martin and O’Leary, 2011; Zagouri et al., 2012). In Indonesia, about 100 new cases occur in 100,000 populations and it is known that 70% of them are in the late stage (Tambunan and Wulandari, 2010). Approximately, 99.7% of worldwide cervical cancer cases are caused by human papillomavirus (HPV) infection (Lin et al., 2009; Steben and Duarte-Franco, 2007). Human papillomavirus (HPV) belongs to the family of simple Papillomaviridae virus. It is non-enveloped, has icosahedral capsid and has genetic material in the form of double-stranded DNA with the length around 7800–7900 base pairs and diameter of 55 nm (Paavonen, 2007). Currently, there are more than 120 HPV genotypes which have been identified and classified based on their risk in generating cervical cancer; there are low risk, probable-high risk and high-risk HPV (Steben and Duarte-Franco, 2007). There are 12 genotypes classified as low-risk HPV (6, 11, 40, 42, 43, 44, 54, 61, 70, 72, 81, and CP6108), three genotypes as probable-high risk HPV (26, 53, and 66) and 15 as high risk HPV (16, 18, 31, 33, 35, 39, 45, 51, 52, 56, 58, 59, 68, 73, and 82). More than 70% of cervical cancer cases in the world are caused by HPV 16 and 18 infection, while 20% of cases are caused by types 31, 33, 35, 45, 52, and 58 (Kanodia et al., 2008; Palmer et al., 2009).
Till date, there is no specific treatment against HPV infection in cervical cancer. Moreover, on a more serious level, surgery and physico-chemo-radiotherapy are also limited (Lin et al., 2009; Tambunan and Wulandari, 2010). One of the treatments currently used for cervical cancer is chemotherapy; it inhibits the HPV activity and suppresses the growth of cervical cancer cell. Cervical cancer treatment currently aims to inhibit histone deacetylases (HDACs) (Lin et al., 2009). HDACs (E.C 3.5.1.98) are non-redundant chromatin modifying enzymes that are often over-expressed in cancer caused by HPV oncoprotein activity, especially E6 and E7 (Szalmás and Kónya, 2009; Trottier and Burchell, 2009). HDACs play an important role in the process of deacetylation or acetyl group release within histone protein tails, that cause histones to be tightly twisted with DNA, thereby inhibiting transcription of tumor suppressor genes and causing tumors to grow into cancer cells (Tambunan and Parikesit, 2012; Tambunan et al., 2011). Thus, HDAC inhibition has become a potential and effective strategy for inducing growth arrest, differentiation and apoptosis of cancer cells in cervical cancer therapy (Rajak et al., 2012).
Currently, variety of compounds has been synthesized and reported to have inhibitory activity against histone deacetylase (HDAC) and anti-proliferation activity against human cancer cells both from in vitro and in vivo evaluation (Fournel et al., 2008). Although there are many classes of HDAC inhibitors (HDACIs), the most potent discovered so far is trichostatin A (TSA) (Nair et al., 2012) and vorinostat or suberoylanilide hydroxamic acid (SAHA); SAHA is the most widely used and belongs to hydroxamic acid group (Tambunan et al., 2011). However, TSA, SAHA and other HDACIs have a clear disadvantage since their production is costly and profound concerns remain about their toxicity, non-specificity and side effects (Nair et al., 2012). Hence, there is a need for the discovery of alternative HDACIs.
Based on pharmacophore modeling, a good HDAC inhibitor has at least three sides/regions: the attachment side of the Zn2+ cofactor/HDAC active site (Zn2+ binding group/ZBG), hydrophobic cap (CAP) and linker containing connecting unit (CU) with electronegative groups (Mohan et al., 2011; Rossi et al., 2011). In this study, the novel HDAC inhibitors (HDACIs) were designed from 2-oxo-1,3-thiazolidine as Zn2+ binding group/ZBG, para-amino benzoic acid (PABA) as linker and acetamide as connecting unit (CU) on hydrophobic cap (CAP) (Marek et al., 2013). Furthermore, functional group variation of these compounds was conducted on hydrophobic cap region and inhibition studies of each varied compound were employed in silico through molecular docking and interaction analysis.
In addition, pharmacology analysis was also done for lead molecules which have better binding affinity against HDAC than standards (Vilar et al., 2008). This analysis includes pharmacological properties based on Lipinski’s rule of five, bioactivity, drug likeness and drug score (Tambunan et al., 2012). Furthermore, in silico preclinical trials including absorption, distribution, metabolism, excretion, and toxicity (ADMET) test, health effect probability and maximum recommended daily dose (MRDD) prediction were also performed in this study (Moroy et al., 2012; van de Waterbeemd and Gifford, 2003). The synthetic accessibility and novelty of the lead compounds were also predicted. This research is expected to obtain novel potential inhibitors of Homo sapiens class II HDAC which can be used as drugs for cervical cancer treatment.
2 Methods
2.1 HDAC sequences and structures searching
H. sapiens class II HDAC enzyme sequences searching was done at the National Center for Biotechnology Information (NCBI) site (http://www.ncbi.nlm.nih.gov). Sequences used in the research were the latest sequence and were supported by valid references. Three-dimensional structure of proteins was derived from the Protein Data Bank (PDB) at the Research Collaboratory for Structural Bioinformatics (RCSB) (http://www.rcsb.org/pdb). If it is not found, the 3D structure could be modeled by using the SWISS Model (Biozentrum Universität Basel) (http://www.swissmodel.expasy.org). All data were stored in.pdb format of structural enzymes. Sequence and structural data were then analyzed, making the enzymes used in this study valid and adequately represent the other enzyme sequences.
2.2 HDAC inhibitor (HDACI) ligands design
All ligands including standard SAHA (suberoylanilide hydroxamic acid), TSA (trichostatin A) and VPA (valproic acid) were designed by using ACDLabs ChemSketch 12.0 software (Advance Chemistry Development, Inc). The refinement and optimization of three-dimensional ligands structure were also performed. Ligand structural data were then stored in MDL molfiles .mol format, so they could be used in docking simulations and pharmacology analysis.
2.3 Molecular docking protocols
Initial preparation for enzymes and ligands was done before molecular docking simulation took place. All heteroatoms (except Zn2+ ion) and water molecules were removed from enzyme file. Polar hydrogen atoms were then added to the enzyme (3D-protonation) to correct oxidation and tautomeric state of amino acid residues. This process was also done to convert Zn0 (catalytically inactive cofactor) to Zn2+ ion (catalytically active cofactor). The enzyme was also minimized using AMBER99 force field. Similarly, ligands were washed/protonated and minimized with MMFF94 force field. Thus, enzymes and ligands were ready to undergo the simulation.
The simulation was performed using the Molecular Operating Environment (MOE) 2008.10 software (Chemical Computing Group, Inc). Docking between the prepared enzymes and ligands was done without solvent (gas phase). The placement, rescoring and refinement mode chosen after the validation of methods were Triangle Matcher, London dG and force field, respectively.
2.4 Docking analysis and interaction visualization
Docking simulations generate data in a MOE database .mdb format. Screening of ligands with better binding affinity against HDAC than the three standards was performed virtually from the lowest free energy of binding (ΔGbinding) value. The best ligands then went pharmacology analysis and molecular dynamics simulation. Visualization of enzyme–ligand complex interactions was also conducted in two- and three-dimensions.
2.5 Pharmacology analysis and preclinical trials
Pharmacology analysis and preclinical trial were conducted to determine the various pharmacological properties and bioactivities of the newly designed inhibitors. Pharmacology analysis in this study included pharmacological properties test based on Lipinski’s rule of five, bioactivity, drug likeness and drug score, whereas preclinical trials included ADMET (absorption, distribution, metabolism, excretion and toxicity) test, health-effect probability and maximum recommended daily dose (MRDD) prediction. These analyses were performed in silico using a variety of software such as ACD/I-Lab, FAF Drugs-2, LAZAR, Molinspiration, Osiris and Toxtree.
2.6 Synthetic accessibility and leads novelty prediction
Synthetic accessibility prediction is needed to assess whether the best designed drug leads are easy to synthesize in the wet laboratory experiment or not (Ertl and Schuffenhauer, 2009). This prediction may use CAESA (SimBioSys, Inc) software. The novelty of lead compounds was predicted on ChemSpider (http://www.chemspider.com) and ZINC compound database (http://zinc.docking.org).
2.7 Molecular dynamics simulation
Molecular dynamics (MD) simulation was performed to assess the stability of the enzyme HDAC with the best drug leads complexes. First, the complex system energy was minimized before MD simulation took place. AMBER99 force field and Generalized Born implicit solvation (GBIS) model were used for this. Initialization time was determined by 100 picoseconds (ps) dynamics at 300 K. The stability of the complex was determined by 5000 ps simulations after it was heated up to 310 K (temperature of normal human body/cervical cancer patients) for 20,000 ps. The simulation was ended with cooling/annealing stage to 1 K for 10 ps to obtain complex conformational structure which has the lowest energy. The Ramachandran plot of the most optimized ligand will be elucidated for determining the integrity of the protein structure.
3 Results and discussion
3.1 Determination of class II HDAC three-dimensional structure
Sequences search for the six H. sapiens class II HDAC (HDAC4, 5, 6, 7, 9, and 10) at the NCBI site produced a total of 65 types of enzymes and their isoforms. From those different HDAC enzyme sequences, one would be selected as a target in molecular docking simulation. Selection of enzyme sequences is based on several parameters such as the length of the corresponding sequence, the use of enzyme sequences in various studies (citation number) and sequences novelty level (newest sequences are the most valid for use). Based on these parameters, the most fulfilling enzyme sequences were obtained from UniProt Knowledge Base (UniProtKB)/SWISS-Prot database. The accession numbers of H. sapiens class II HDAC sequences used as targets are P56524.3 (HDAC4), Q9UQL6.2 (HDAC5), Q9UBN7.2 (HDAC6), Q8WUI4.2 (HDAC7), Q9UKV0.2 (HDAC9) and Q969S8.1 (HDAC10). Moreover, three-dimensional structure and crystal data of the HDACs could be searched and obtained from its accession code at the RCSB-PDB database. PDB IDs obtained are 2VQJ (HDAC4: P56524.3), 3PHD (HDAC6: Q9UBN7.2) and 3C0Y (HDAC7: Q8WUI4.2). The HDAC crystal structures were downloaded in .pdb format. Three other enzymes (HDAC5, HDAC9, and HDAC10) 3D structures were not obtained from PDB and needed to be modeled first. Three-dimensional structures modeling of proteins from their amino acid sequences in FASTA format could be done by using SWISS Model software (Arnold et al., 2006). Modeled enzyme structures were then stored in .pdb format to be applied in the MOE software for simulation. As a result of structure searching and modeling, the visualization of HDACs catalytic area is seen in Fig. 1.![The representative structure of HDACs catalytic area. (a) HDAC4 [PDB: 2VQJ]. (b) HDAC5 [UniProtKB/SWISS-Prot: Q9UQL6]. (c) HDAC6 [PDB: 3PHD]. (d) HDAC7 [PDB: 3C0Y]. (e) HDAC9 [UniProtKB/SWISS-Prot: Q9UKV0]. (f) HDAC10 [UniProtKB/SWISS-Prot: Q969S8].](/content/184/2019/12/2/img/10.1016_j.arabjc.2015.07.010-fig1.png)
In addition, for analytical purpose, it was necessary that the validity of the HDAC three-dimensional structures from SWISS Model be verified first; hence it was valid to employ it in the simulation. The structure validity could be determined from the QMEAN4 score and QMEAN Z-score values (Benkert et al., 2011). QMEAN4 score is based on the energy contribution to the protein structure. On the other hand, QMEAN Z-score is based on the structure comparison between modeled structures and experimental structures which use high resolution X-ray crystallography that was obtained from database. Modeled structures are valid to use if the QMEAN4 score value is between 0 and 1. They are also valid if QMEAN Z-score value is between −1 and 0, which indicates the best structure. Meanwhile, the value between −2 to −1 is acceptable, but if it is less than −3, the structure needs to be corrected. QMEAN4 scores for HDAC5, 9, and 10 are 0.622, 0.680, and 0.577 respectively. Based on QMEAN4 score, the three modeled HDAC structures are valid for use. Despite having a valid QMEAN4 score, HDAC10 modeled structure is less valid based on QMEAN Z-score with a value of −3.22. However, the structure still could be used because it has a good structural energy based on QMEAN4 score and the absence of protein structural data that is quite similar to HDAC10 in protein database.
3.2 HDAC multiple sequence alignment (MSA) and residue conservation
Multiple sequence alignment (MSA) was performed by using the Basic Local Alignment Search Tool (BLAST). It could be performed through the NCBI BLAST server (http://blast.ncbi.nlm.nih.gov/Blast.cgi). Based on the MSA results, the entire class II human HDAC enzymes belong to the HDAC arginase superfamily (Bottomley et al., 2008; Schuetz et al., 2008), except HDAC6 that belongs to the zinc-finger ubiquitin binding protein (ZnF-UBP) (Ouyang et al., 2012). According to these results, each enzyme sequence selected for docking simulation was able to represent corresponding sequences of the other enzymes. It could be inferred that similarity percentage between corresponding HDAC sequences and structures was also high, which varied between 99% and 100%. Therefore, the selection of sequences of enzymes used in docking simulations in this study was valid.
Furthermore, residue conservation was also conducted to analyze the conserved residue between class II HDAC. It could be performed by using Clustal-Omega at the EMBL-EBI server (http://www.ebi.ac.uk/Tools/msa/clustalo/). The results of HDAC residues conservation are shown in Fig. 2 which shows highly conserved residues between each H. sapiens class II HDAC.
3.3 Characterization of HDAC catalytically active residues
HDAC active site is needed to be determined first to begin the molecular docking simulation. For HDAC structures that are derived from PDB (HDAC4, HDAC6, and HDAC7), active residues could be known directly. Meanwhile, for the HDAC structures that are derived from the SWISS Model, active residues could be identified by interaction analysis between enzyme and inhibitor ligand from LigX interaction feature on MOE (Molecular Operating Environment) according to the template and site finder.
Generally, HDACs have two regions that are functionally important which are carboxy/C-terminal and conserved catalytic region of amino/N-terminal extension that have no similarity with the other proteins (Schuetz et al., 2008). All HDACs, except class III HDAC are zinc-dependent which needs Zn2+ ion as a cofactor. The class II HDAC active site forms a cylindrical/tubular cavity with Zn2+ ion cofactor at the bottom of the active site. In class II HDAC, a Zn2+ cofactor binds to charge-relay system that consists of two aspartate residues (Asp/D) and one histidine residue (His/H), except HDAC6 where the two aspartate residues are replaced by two cysteine residues (Cys/C). For example, visualization of cavity and charge-relay system position in HDAC5 is seen in Fig. 3. Based on LigX interactions analysis, charge-relay systems of each H. sapiens class II HDAC are seen in Table 1.
| Enzyme | Charge-relay system position |
|---|---|
| HDAC4 | Zn2+ ion, with Asp196, His198, and Asp290 residues |
| HDAC5 | Zn2+ ion, with Asp870, His872, and Asp964 residues |
| HDAC6 | Zn2+ ion, with Cys5, His7, and Cys78 residues |
| HDAC7 | Zn2+ ion, with Asp707, His709, and Asp801 residues |
| HDAC9 | Zn2+ ion, with Asp820, His822, and Asp914 residues |
| HDAC10 | Zn2+ ion, with Asp172, His174, and Asp265 residues |
3.4 HDAC inhibitor (HDACI) ligands design
In this study, HDAC inhibitors were designed by de novo approach from various pharmacophore fragments which are shown to have a good inhibitory activity against HDACs. De novo drug design is an approach where the structure of a chemical compound is formed from atoms or fragments corresponding to the protein target binding sites (Jain et al., 2011). Pharmacophore model for HDAC inhibitors consists of three main parts: zinc binding group (ZBG), linker and hydrophobic cap (CAP). The role of ZBG is to bind Zn2+ cofactor at enzyme active site that directly inhibits enzymatic activity. Linker domain resembles the substrate and is able to bind the cylindrical pocket of the HDAC active site. Meanwhile, CAP interacts with the surface and closes the cylindrical pocket of the HDAC active site. From this model, the compound series based on 4-[(2-oxo-1,3-thiazolidin-3-yl)carbonyl] aniline has been designed and modified to obtain the best drug leads for cervical cancer.
In this research, 2-oxo-1,3-thiazolidine fragment was selected for ZBG. It is a bioisosterical form of 1,3-thiazolidine-2,4-dione which has been shown to have good inhibitory activity against HDAC Zn2+ cofactor in liver cancer cell lines (Mohan et al., 2011). As a linker, para-aminobenzoic acid (PABA) was chosen as it has been shown to have strong interaction with hydrophobic residues on HDAC cylindrical pocket. Its modified compounds also appeared to have a good inhibitory activity against HDACs both in vitro and in vivo evaluation (Rajak et al., 2012). Then as a hydrophobic cap (CAP), acetamide group was chosen as connecting unit (CU) and a hydrophobic group –R that may interact with the HDAC cavity surface (Rajak et al., 2012). The general structure of the class II HDAC inhibitor in this study is seen in Fig. 4.
Modification target of designed HDAC inhibitors is the –R group on hydrophobic cap. Hydrophobic –R group substitutes include cyclohexane, cyclopentane, naphthalene, pyridine, pyrrole, furan, thiophene, pyrimidine, imidazole, oxazole, thiazole, triazole, quinoline, indole and polysubstituted benzene. In addition to these fragments, hydrophobic cap also varied by other groups such as halides, alkyl, aryl, hydroxyl, amine, ether and carbonyl. This design is expected to increase the hydrophobicity and binding affinity between ligand inhibitors with class II HDAC in human. As a result, HDAC inhibitor modifications yielded a total of 2688 ligands.
3.5 Docking analysis and virtual screening
Database containing free energy of binding of the complex from molecular docking simulation was screened based on the lowest ΔGbinding value. Most of the designed ligands have lower ΔGbinding value than SAHA, TSA, and VPA standards. But only eight best ligands were chosen as novel potential HDAC inhibitors. They were ligand Aa50 (–R = 2-propanoyl phenyl), Bb38 (–R = 3-carboxy cyclohexyl), F (–R = pyridin-2-yl), Gc42 (–R = 5-(hydroxyamino)-1H-pyrrol-2-yl]), Gd40 (–R = 3-(hydroxycarboxamide)-1H-pyrrol-2-yl), Ib14 (–R = 4-methoxythiophen-2-yl), Kb17 (–R = 5-acetyl-1H-imidazol-2-yl) and O38 (–R = 5-carboxyl-1,2,3-oxadiazol-4-yl). The structures of standards and eight best ligands are shown in Fig. 5. The best designed ligands free energy of binding calculation with corresponding class II HDAC results are shown in Table 2. For example, one of the best ligand that has strong interaction is Ib14 with HDAC5 (visualized in Fig. 6).
| Ligand | Gibbs free energy of binding, ΔGbinding (kcal/mol) | |||||
|---|---|---|---|---|---|---|
| HDAC4 | HDAC5 | HDAC6 | HDAC7 | HDAC9 | HDAC10 | |
| SAHA | −5.8131 | −12.2155 | −23.4231 | −6.3467 | −9.6285 | −8.2174 |
| TSA | −11.0343 | −7.7395 | −22.1517 | −6.9128 | −7.7594 | −9.4841 |
| VPA | −7.9939 | −13.5089 | −28.5147 | −6.5775 | −17.9370 | −1.0557 |
| Aa50 | −11.5372 | −13.3201 | −20.9017 | −14.6737 | −9.1748 | −13.9624 |
| Bb38 | −8.8182 | −14.7906 | −35.7428 | −9.2168 | −27.3834 | −7.6069 |
| F | −18.4812 | −9.8521 | −22.8164 | −6.5852 | −6.3970 | −9.4943 |
| Gc42 | −20.2563 | −12.4177 | −20.7936 | −8.4926 | −8.5369 | −10.2846 |
| Gd40 | −9.2588 | −7.9550 | −23.9998 | −7.3624 | −29.1305 | −9.2218 |
| Ib14 | −18.1408 | −18.1124 | −19.5423 | −7.8312 | −13.3199 | −13.4759 |
| Kb17 | −10.2072 | −12.6084 | −22.9797 | −15.1957 | −11.3457 | −9.4706 |
| O38 | −12.2219 | −17.7518 | −31.0753 | −11.0377 | −25.6119 | −10.9574 |
![Ligand Ib14 binding to HDAC5. (a) Electron density map of the Ib14 ligand (green mesh). (b) Electron density corresponding to the Ib14 ligand that is shown as sticks (N-(4-methoxythiophen-2-yl)-2-({4-[(2-oxo-1,3-thiazolidin-3-yl)carbonyl] phenyl}amino)acetamide). (c) Surface representation of the HDAC5–Ib14 ligand complex. The active site/cavity of HDAC5 is indicated by a dashed yellow circle. (d) Ib14 ligand (shown as sticks) entering the HDAC5 cavity (indicated by a dashed yellow circle) in 12 Å sphere. (e) Three-dimensional visualization of detailed interaction of Ib14 (green sticks) with HDAC5 (ribbon). Side chains from HDAC5 that make contact with Ib14 are labeled and gray colored. Hydrogen bond and electrostatic interaction are indicated by dashed lines.](/content/184/2019/12/2/img/10.1016_j.arabjc.2015.07.010-fig6.png)
Fig. 5(e) shows that the interaction between the ligands with enzyme followed the hypothesis. Carbonyl (C⚌O) group of the 2-oxo-1,3-thiazolidine fragment acted as ZBG that strongly interacts with Zn2+ cofactor on the charge-relay system (Asp870, His872, and Asp964) and reaches 100.0% interaction (data in Table 8). This was also the case with the other best ligands. Methoxythiophene fragment acted as a hydrophobic cap and interacted with Gly962 residue on HDAC cavity surface. PABA as a linker was also contributed to bind HDAC cylindrical pocket. Similar interactions were also shown by the other best ligands.
3.6 Pharmacology analysis
Initial pharmacology analysis in this study was pharmacological properties testing based on Lipinski’s rule of five. Lipinski’s rule of five is quite a good approach for drug lead design, especially in the hydrophobicity and oral bioavailability prediction (Morya et al., 2012; Oprea, 2002). This rule is based on the approach that the human body contains about 90% of water. Drug molecules should be hydrophobic enough to last longer in the body and not easily excreted, as this will enhance its bioactivity efficiency in the body. Modified Lipinski’s rule of five consists of the following: molecular weight (MW) is from 160 to 500 Da, octanol–water partition coefficient (log P) is between −0.4 and +5.6, molar refractivity is from 40 to 130 cm3, hydrogen donor is less than 5 atoms, hydrogen acceptor is less than 10 atoms, and topological polar surface area (TPSA) is not more than 140 Å2 (Ekman et al., 2005). These pharmacological properties of standards (SAHA, TSA, and VPA) and designed ligands could be tested by using Molinspiration and Osiris online software (Table 3).
| Ligand | Molecular descriptor | ||||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|
| MW | LogP | TPSA | n atoms | n ON | n OHNH | n viol | n rotb | Sol | DrugL | DrugS | |
| SAHA | 264.325 | 2.467 | 78.424 | 19.0 | 5 | 3 | 0 | 8 | -3.33 | -8.87 | 0.35 |
| TSA | 302.374 | 2.683 | 69.635 | 22.0 | 5 | 2 | 0 | 6 | -3.26 | 1.24 | 0.37 |
| VPA | 144.214 | 2.804 | 37.299 | 10.0 | 2 | 1 | 0 | 5 | -1.97 | -2.62 | 0.3 |
| Aa50 | 411.483 | 2.538 | 95.576 | 29.0 | 7 | 2 | 0 | 7 | -4.72 | 3.24 | 0.39 |
| Bb38 | 405.476 | 1.403 | 115.804 | 28.0 | 8 | 3 | 0 | 6 | -3.46 | 0.98 | 0.41 |
| F | 356.407 | 1.285 | 91.397 | 25.0 | 7 | 2 | 0 | 5 | -3.49 | 4.44 | 0.49 |
| Gc42 | 484.319 | 2.443 | 94.296 | 26.0 | 7 | 3 | 0 | 6 | -3.55 | 3.69 | 0.48 |
| Gd40 | 403.42 | 0.625 | 132.765 | 28.0 | 10 | 4 | 0 | 5 | -3.89 | 3.82 | 0.46 |
| Ib14 | 391.474 | 1.953 | 87.739 | 26.0 | 7 | 2 | 0 | 6 | -3.83 | 3.8 | 0.46 |
| Kb17 | 403.42 | 0.762 | 133.493 | 28.0 | 10 | 3 | 0 | 7 | -3.67 | 3.66 | 0.46 |
| O38 | 391.365 | 0.071 | 154.728 | 27.0 | 11 | 3 | 1 | 6 | -3.41 | 3.71 | 0.48 |
As seen in Table 3, all best ligands meet the Lipinski’s rule of five, except O38 ligand which has insignificant violations in TPSA and numbers of hydrogen acceptor. Despite having violations, O38 ligand could still be used as potential HDAC inhibitor for having a better drug likeness and drug score than standards. We also see that all best ligands have a lower solubility in water than the three standards. This result indicated that all best ligands have a good hydrophobicity, will be able to circulate longer and are not easily excreted out from the body in order to enhance their activities and interactions with target.
One of the pharmacology parameters that were also analyzed in this study is drug likeness. The drug likeness is a pharmacological parameter that indicates fragments and bioactivity similarity between designed ligands and already known drugs (Hartenfeller et al., 2012). The larger value means the higher similarity and probability that the particular molecule will be active. The drug likeness score is calculated from the combined score of several bioactivity parameters similarity such as GPCR (G-protein coupled receptor) ligand, ion channel modulator, kinase inhibitor, nuclear receptor ligand, protease inhibitor and enzyme inhibitor properties (Kubinyi, 2006). From Table 3, it is shown that all best ligands have a larger drug likeness value than standards, except Bb38 ligand which has a little lower drug likeness value than TSA standards. We may conclude from this result that all best designed ligands have high fragments and bioactivity similarity with existing drugs.
Meanwhile, the drug score is a combination of several pharmacological parameters such as log P (octanol–water partition coefficient), solubility, drug likeness, and toxicity parameters prediction (mutagenicity, tumorigenicity, irritant, and reproductive effect) within one useful practical value (Singh and Tomar, 2011). The higher drug score value means that the compounds have a better chance to be developed as drug leads. Table 3 shows that all eight best ligands have a better drug score than the standards and could be developed as novel drugs for cervical cancer treatment.
3.7 In silico preclinical trial
The earliest in silico preclinical trial of the designed inhibitor ligands in this study is toxicity prediction. Two applications were used in this study to predict the toxicological properties of the lead compounds such as mutagenicity and carcinogenicity. They are Toxtree v2.1.0 and LAZAR which have several differences to determine the toxicological properties of the ligands.
Toxtree determined the toxicological properties based on Benigni–Bossa rules. This rule could evaluate compound toxicity based on certain functional groups or chemical fragments that may induce mutagenicity or carcinogenicity (Benigni et al., 2008, 2007). The results of ligands toxicity prediction by using Toxtree are seen in Table 4. From the results, we could conclude that TSA standards and Gc42 ligand have potential to be genotoxic, carcinogenic and Bb38 ligand has potential to be mutagenic to Salmonella typhimurium TA100. Gd40, F, Kb17, Aa50, O38 and Ib14 ligands are in accordance with Benigni–Bossa rules because they do not have any tendency to become mutagenic or carcinogenic.
| Ligand | Negative for genotoxic carcinogenicity | Negative for nongenotoxic carcinogenicity | Potential S. typhimurium TA100 mutagen | Potential carcinogen |
|---|---|---|---|---|
| SAHA | Yes | Yes | No | No |
| TSA | No | Yes | Yes | No |
| VPA | Yes | Yes | No | No |
| Gc42 | No | Yes | No | No |
| Gd40 | Yes | Yes | No | No |
| Bb38 | Yes | Yes | Yes | No |
| F | Yes | Yes | No | No |
| Kb17 | Yes | Yes | No | No |
| Aa50 | Yes | Yes | No | No |
| O38 | Yes | Yes | No | No |
| Ib14 | Yes | Yes | No | No |
The next toxicity prediction is using LAZAR (Lazy Structure–Activity Relationship). This software could evaluate toxicological properties of the ligands based on chemical fragments similarity with already known toxic fragments from compounds database (Maunz and Helma, 2008). The results are shown in Table 5. Compounds mutagenicity was evaluated from DSS-Tox (Distributed Structure-Searchable Toxicity) database and S. typhimurium assay test based on Kazius–Bursi method. From this evaluation, we may see that almost all designed ligands have mutagenic properties based on DSS-Tox database and S. typhimurium assay test, except Bb38 and Aa50 ligands which have no mutagenic properties based on DSS-Tox DBS. Moreover, due to its mutagenic properties, it does not always lead to carcinogenicity; the ligands carcinogenic properties were also evaluated by animal assay test, which are hamster, mouse, and rat. We can see in Table 5 that all best ligands were not carcinogenic to those animals.
| Ligands | Mutagenicity | Carcinogenicity | ||||||
|---|---|---|---|---|---|---|---|---|
| DSS-Tox DBS | Salmonella Kazius–Bursi | Hamster | Mouse | Rat | Single cell call | Multiple cell call | ISSCAN v3a Canc | |
| SAHA | Yes | Yes | No | No | No | No | No | Yes |
| TSA | No | Yes | No | Yes | No | No | No | No |
| VPA | No | No | No | No | No | No | No | Yes |
| Gc42 | Yes | Yes | No | No | No | Yes | No | No |
| Gd40 | Yes | Yes | No | No | No | Yes | No | No |
| Bb38 | No | Yes | No | No | No | No | No | No |
| F | Yes | Yes | No | No | No | Yes | No | No |
| Kb17 | Yes | Yes | No | No | No | Yes | No | No |
| Aa50 | No | Yes | No | No | No | No | No | No |
| O38 | Yes | Yes | No | No | No | Yes | No | No |
| Ib14 | Yes | Yes | No | No | No | Yes | Yes | No |
Carcinogenic properties were also calculated by using ISSCAN (Instituto Superiore di Sanita Chemical Carcinogens) method which is developed by Benigni and Bossa from the Department of Environment and Primary Prevention, United States Environmental Protection Agency (US-EPA) (Benigni et al., 2007). From this evaluation, all best designed ligands have no carcinogenic properties, especially as cancer causative agents. In contrast, SAHA and VPA standards may cause cancers based on this method.
The last stage of in silico preclinical trial was ADMET test including absorption, health effect probability, and oral bioavailability prediction using ACD/I-Lab (Lagorce et al., 2008; van de Waterbeemd and Gifford, 2003). Absorption prediction yielded maximum passive absorption percentage from the contribution of transcellular and paracellular routes. Passive absorption is a process in which the molecules are absorbed through the intestinal mucosa or epithelial cells that do not use energy or ATP contribution along the process. Passive absorption is divided into two major sub-types based on the route it takes: transcellular and paracellular.
In transcellular route absorption, the compounds cross the cell membranes and enter the cell. Meanwhile, in paracellular route, the compounds cross the tight junction between cells (Pendekal and Tegginamat, 2013). Transcellular route was the favored pathway for most drugs, probably due to its larger accessible surface area. The paracellular permeability of a drug was size- and charge-dependent, which is decreased with increasing molecular size and negative charge (Agoram et al., 2001).
Table 6 shows that all best ligands, except Gd40 and O38, have a high maximum passive absorption percentage with contribution only (or dominantly, for Kb17 ligand) from transcellular route, the favored pathway for most drugs. We may conclude then that all best ligands except Gd40 and O38 have a good passive absorption, which is an essential parameter for drug development.
| Ligand | Max. passive absorption (Trans/Para) | Health effect probability | Oral bioavailability (%) | MRDD (mg/kg/day) | |||||
|---|---|---|---|---|---|---|---|---|---|
| Blood | Cardio vascular | Gastro intestinal | Kidney | Liver | Lungs | ||||
| SAHA | 100 (100/0) | 0.36 | 0.25 | 0.07 | 0.11 | 0.11 | 0.37 | >70 | 8.49 |
| TSA | 100 (100/0) | 0.59 | 0.5 | 0.51 | 0.27 | 0.52 | 0.63 | 30–70 | 10.05 |
| VPA | 100 (100/0) | 0.09 | 0.09 | 0.08 | 0.05 | 0.05 | 0.05 | >70 | 15.46 |
| Gc42 | 100 (100/0) | 0.51 | 0.11 | 0.78 | 0.54 | 0.74 | 0.11 | <30 | 4.55 |
| Gd40 | 20 (59/41) | 0.42 | 0.14 | 0.24 | 0.74 | 0.26 | 0.76 | <30 | 13.33 |
| Bb38 | 100 (100/0) | 0.85 | 0.04 | 0.5 | 0.42 | 0.56 | 0.13 | <30 | 1.51 |
| F | 100 (100/0) | 0.42 | 0.33 | 0.65 | 0.36 | 0.25 | 0.25 | 30–70 | 4.95 |
| Kb17 | 98 (99/1) | 0.61 | 0.43 | 0.23 | 0.51 | 0.7 | 0.18 | <30 | 4.52 |
| Aa50 | 100 (100/0) | 0.31 | 0.23 | 0.51 | 0.63 | 0.31 | 0.42 | <30 | 6.99 |
| O38 | 54 (99/1) | 0.47 | 0 | 0.16 | 0.37 | 0.7 | 0.1 | <30 | 11.88 |
| Ib14 | 100 (100/0) | 0.18 | 0.36 | 0.09 | 0.25 | 0.15 | 0.11 | 30–70 | 11.23 |
Health effect probability prediction was also conducted in this study by using ACD/I-Lab. The probability value ranges between 0 and 1; increasing probability reached 1 (one)-which means the compounds have a higher health effect to corresponding organs or organ systems. The threshold value of health effect probability is 0.85; the probability score below this value is still acceptable. We can see in Table 5 that all best ligands have a low health effect probability, except Bb38 ligand which has relatively high health effect probability that reaches threshold value on blood. On average, the ligand that has the lowest health effect probability is Ib14. Due to its health effect, ACD/I-Lab could also predict the maximum recommended daily dose (MRDD) if the designed ligands will be developed as cervical cancer drugs.
Based on the hydrophobicity of the designed ligands, the oral bioavailability could also be calculated using FAF-Drugs2 and ACD/I-Lab software. From FAF-Drugs2 prediction, all best ligands have a good oral bioavailability based on Veber’s and Egan’s rules (Lagorce et al., 2008). The best ligands that have the highest oral bioavailability are F and Ib14 with oral bioavailability that ranges between 30% and 70%.
3.8 Synthetic accessibility and ligands novelty prediction
The next step in drug discovery after computer-aided drug discovery and development (CADDD) is clinical test. For further clinical trial, the best designed ligands from in silico approach have to be synthesized first to determine their in vitro and in vivo bioactivity. For this reason, synthetic accessibility of in silico designed ligands then has to be assessed to evaluate the synthetic probability and feasibility of the ligands.
The bases of this prediction method are screening and informational analysis of the contribution of particular fragments and the complexity of the structures derived from database containing millions of compounds that are already synthesized. So from this method, the synthetic accessibility of drug lead compounds could be estimated from de novo synthesis approach (Ertl and Schuffenhauer, 2009). Software that could be used to predict the synthetic possibility by using this method is CAESA (Computer Assisted Estimation of Synthetic Accessibility). These tools can estimate the percentage of synthetic accessibility based on the existence of starting materials and structural or residual complexity (Ertl and Schuffenhauer, 2009). Synthetic accessibility prediction of the best ligands is presented in Table 7.
| Ligands | Residual complexity (%) | Starting materials (%) | Synthetic accessibility (%) |
|---|---|---|---|
| SAHA | 3 | 75 | 93 |
| TSA | 0 | 100 | 100 |
| VPA | 0 | 100 | 100 |
| Gc42 | 17 | 49 | 75 |
| Gd40 | 18 | 46 | 74 |
| Bb38 | 7 | 72 | 88 |
| F | 4 | 55 | 82 |
| Kb17 | 16 | 56 | 77 |
| Aa50 | 13 | 60 | 81 |
| O38 | 19 | 48 | 74 |
| Ib14 | 15 | 43 | 76 |
TSA and VPA standards have easy synthetic accessibility percentage that reaches 100%. It is because they have been synthesized and already listed on Acros, Aldrich, or Lancaster database. All best ligands have a low residual complexity and a high synthetic accessibility. Based on this result, ligands synthesis could be started from the starting materials that have been synthesized and registered in the database.
All designed inhibitor ligands are predicted to be compounds that have never been synthesized before. Novelty evaluation could be done by searching ligand chemical structures on ChemSpider and ZINC compounds database. Through identical and similar structures, tautomers, structural and geometrical isomers search, all designed ligands were not found on both databases. According to this result, it could be predicted that the designed inhibitors were new compounds and have never been synthesized before.
3.9 Drug screening based on docking, dynamics and drug scan analysis
The best inhibitor ligands could be determined based on docking and drug scan analysis (pharmacology analysis and preclinical trial). For HDAC4, HDAC5, HDAC6, HDAC7, HDAC9 and HDAC10, the best inhibitor ligands are Gc42, Ib14, O38, Kb17, Gd40, and Aa50 respectively due to their lowest free energy of binding compared with other ligands. Although Gc42 ligand is the best inhibitor for HDAC4, it is predicted to have genotoxic carcinogenicity based on Benigni–Bossa rules. It is also found to have a high health effect probability and low oral bioavailability. The other ligand that also has a good binding affinity to HDAC4 is F ligand. Nevertheless, based on drug scan, F ligand does not have any mutagenic–carcinogenic properties and has better pharmacological properties such as health effect probability and oral bioavailability than Gc42 ligand. Thus, F ligand was then chosen as the best inhibitor for HDAC4.
The other best ligands (Ib14, O38, Kb17, Gd40, and Aa50 for HDAC5, HDAC6, HDAC7, HDAC9, and HDAC10 respectively) are already in accordance with Lipinski’s rule of five; they have lowest mutagenic–carcinogenic properties, a good passive absorption, health effect probability score and oral bioavailability based on Veber’s and Egan’s rules.
The structure and property differences among the best ligands of respective H. sapiens class II HDAC are a consequence of characteristics differences among HDAC. Structural and amino acid residues differences in HDAC catalytic site cause tendency in interaction differences.
3.10 Molecular dynamics simulation
Protein does not have a fixed structure because of rapid conformational changes in the presence of solvent. This is so because the atoms movement in protein backbone causes changes of α-helix, β-sheet and loop conformation as protein domains (Levinson et al., 2006). During folding and unfolding process, the stability of the enzyme-inhibitor complex also changes to include inhibitor interaction with enzyme active site. Molecular dynamics simulation then aims to evaluate this change.
The selected best ligands to undergo this simulation according to the drug screening result were F with HDAC4, Ib14 with HDAC5, O38 with HDAC6, Kb17 with HDAC7, Gd40 with HDAC9, and Aa50 with HDAC10. Simulation was done at pico to nanoseconds time scale (10−12–10−9 s) as domain protein changes could be observed at this range. First, initialization time was calculated to determine the complex conformational energy and time needed for stability in the presence of implicit solvent. It is run for 100 ps at 300 K. At the end of the simulation, the entire ligands still interact with the enzyme.
MD simulation was then continued for 5000 ps. The heating for 10 ps from 300 K to 310 K was performed first and ended with cooling for 10 ps until 1 K. At the end of the simulation, all complexes were allowed to remain stable and all best ligands were still strongly interacted with the enzyme. The complex stability could be evaluated from RMSD (root mean square deviation) value along the simulation. The RMSD curve of all complexes is presented in Fig. 7. We can see in Fig. 7 that all complexes are stable in respective RMSD value. However, we took ligand interaction between HDAC4-F as a simulation proxy, as it is having the most stable RMSD value, and re-run the dynamics simulation for 20,000 ps. As shown in Fig. 8, the dynamics simulation is beginning to retain its stability starting at 4000 ps. Hence, the Ramachandran plot at Fig. 9 shows that the residues of the enzyme are still within the 15% threshold in the disallowed region.
Additionally, all best ligands also tend to have a better interaction with HDACs after molecular dynamics simulation takes place. The interaction data of the best ligands with corresponding HDACs in each MD simulation stage are presented in Table 8. According to these results it could be concluded that best ligand inhibitors have a strong interaction with enzyme. Consequently, the formed complexes remained stable despite being influenced by solvent.
| Complex | The best interactions between ligand and HDAC in molecular dynamics simulation | |
|---|---|---|
| HDAC4–F | Before simulation (docking) | Initialization 100 ps (300 K) |
| Glu329: 36.1%; H-acc; 2.58 Å | Zn-Asp196: 41.6%; ionic; 1.80 Å | |
| Zn-Asp196: 90.4%; ionic; 2.30Å | Zn-Asp196: 34.9%; ionic; 1.77 Å | |
| Zn-Asp196: 76.9%; ionic; 1.99 Å | Zn2+: 100.0%; ionic; 2.70 Å | |
| Zn-His198: 19.7%; ionic; 1.94 Å | ||
| Zn-Asp290: 21.4%; ionic; 1.91 Å | ||
| Zn2+: 100.0%; ionic; 2.69 Å | ||
| Equilibration – heating 10 ps | Simulation 5000 ps – cooling 10 ps | |
| Ser287: 24.5%; H-donor; 1.67 Å | Leu326: 13.8%; H-donor; 1.97 Å | |
| Ser287: 12.2%; H-acc; 2.98 Å | Zn-Asp196: 48.7%; ionic; 1.84 Å | |
| Zn-Asp196: 30.3%; ionic; 2.03 Å | Zn-Asp196: 48.6%; ionic; 1.84 Å | |
| Zn-His198: 23.8%; ionic; 1.79 Å | Zn-His198: 26.8%; ionic; 2.10 Å | |
| Zn2+: 100.0%; ionic; 2.15 Å | Zn-Asp290: 56.2%; ionic; 1.88 Å | |
| Zn-Asp290: 42.0%; ionic; 1.81 Å | ||
| Zn2+: 100.0%; ionic; 1.86 Å | ||
| HDAC5–Ib14 | Before simulation (docking) | Initialization 100 ps (300 K) |
| Ser960: 29.9%; H-acc; 2.73 Å | Zn-His872: 32.7%; ionic; 1.78 Å | |
| Asp964: 40.8%; H-acc; 2.80 Å | Zn-Asp964: 62.4%; ionic; 1.92 Å | |
| Zn-Asp870: 93.8%; ionic; 2.10 Å | Zn-Asp964: 45.3%; ionic; 1.82 Å | |
| Zn-Asp870: 86.8%; ionic; 2.35 Å | Zn2+: 100.0%; ionic; 1.97 Å | |
| Zn-His872: 35.0%; ionic; 1.96 Å | ||
| Zn-Asp964: 28.3%; ionic; 2.46 Å | ||
| Zn2+: 100.0%; ionic; 2.46 Å | ||
| Equilibration – heating 10 ps | Simulation 5000 ps – cooling 10 ps | |
| Ala961: 12.6%; H-acc; 3.17 Å | Gly831: 19.7%; H-donor; 1.79 Å | |
| His1006: 31.6%; H-acc; 2.83 Å | Ser960: 18.6%; H-donor; 1.58 Å | |
| Zn-Asp870: 50.3%; ionic; 1.85 Å | Gly831: 41.3%; H-acc; 2.85 Å | |
| Zn-Asp870: 41.8%; ionic; 1.80 Å | Asn846: 33.0%; H-acc; 2.59 Å | |
| Zn-His872: 28.0%; ionic; 1.77 Å | Ser960: 20.9%; H-acc; 3.09 Å | |
| Zn-Asp964: 11.1%; ionic; 1.73 Å | Zn-Asp870: 19.6%; ionic; 2.18 Å | |
| Zn2+: 100.0%; ionic; 1.86 Å | Zn-His872: 26.4%; ionic; 1.80 Å | |
| Zn-Asp964: 17.1%; ionic; 1.86 Å | ||
| Zn2+: 100.0%; ionic; 2.18 Å | ||
| HDAC6–O38 | Before simulation (docking) | Initialization 100 ps (300 K) |
| Ala10: 51.7%; H-acc; 2.82 Å | Zn-His7: 21.3%; ionic; 1.72 Å | |
| Zn-His7: 62.4%; ionic; 1.98 Å | Zn-Cys78: 15.1%; ionic; 2.04 Å | |
| Zn-Cys78: 11.6%; ionic; 2.38 Å | Zn2+: 100.0%; ionic; 1.84 Å | |
| Zn2+: 100.0%; ionic; 2.18 Å | ||
| Zn2+: 100.0%; ionic; 2.80 Å | ||
| Equilibration – heating 10 ps | Simulation 5000 ps – cooling 10 ps | |
| Zn-His7: 33.6%; ionic; 1.78 Å | His7: 45.3%; H-donor; 1.36 Å | |
| Zn-Cys78: 13.1%; ionic; 1.96 Å | Zn-His7: 34.0%; ionic; 1.79 Å | |
| Zn2+: 100.0%; ionic; 1.84 Å | Zn-Cys78: 88.7%; ionic; 2.51 Å | |
| Zn2+: 100.0%; ionic; 2.07 Å | ||
| HDAC7–Kb17 | Before simulation (docking) | Initialization 100 ps (300 K) |
| His669: 36.3%; H-donor; 1.79 Å | Asp707: 12.6%; H-donor; 2.19 Å | |
| Ser728: 59.0%; H-acc; 2.50 Å | Leu729: 21.7%; H-acc; 2.99 Å | |
| Zn-Asp707: 93.9%; ionic; 2.27 Å | Zn-Asp801: 12.1%; ionic; 1.92 Å | |
| Zn-Asp707: 86.9%; ionic; 2.04 Å | Zn2+: 100.0%; ionic; 1.74 Å | |
| Zn-His709: 55.6%; ionic; 1.93 Å | ||
| Zn-Asp801: 27.3%; ionic; 1.92 Å | ||
| Zn2+: 100.0%; ionic; 2.60 Å | ||
| Equilibration – heating 10 ps | Simulation 5000 ps – cooling 10 ps | |
| His670: 53.2%; H-acc; 2.60 Å | Tyr726: 22.1%; H-donor; 2.88 Å | |
| Zn-Asp707: 92.8%; ionic; 2.10 Å | Tyr726: 22.1%; H-acc; 2.88 Å | |
| Zn-Asp707: 88.1%; ionic; 2.06 Å | Leu729: 29.3%; H-acc; 2.86 Å | |
| Zn-Asp801: 87.1%; ionic; 2.34 Å | Zn-Asp707: 60.7%; ionic; 1.90 Å | |
| Zn-Asp801: 83.7%; ionic; 2.03 Å | Zn-Asp707: 60.2%; ionic; 1.90 Å | |
| Zn2+: 100.0%; ionic; 2.07 Å | Zn-His709: 23.5%; ionic; 1.72 Å | |
| Zn-Asp801: 13.8%; ionic; 1.73 Å | ||
| Zn2+: 100.0%; ionic; 1.92 Å | ||
| HDAC9–Gd40 | Before simulation (docking) | Initialization 100 ps (300 K) |
| Gly955: 66.0%; H-acc; 2.69 Å | Zn-Asp820: 85.2%; ionic; 2.03 Å | |
| Zn-Asp820: 59.6%; ionic; 2.12 Å | Zn-Asp820: 80.7%; ionic; 2.01 Å | |
| Zn-His822: 10.1%; ionic; 3.15 Å | Zn-Asp914: 48.5%; ionic; 2.04 Å | |
| Zn-Asp914: 39.0%; ionic; 1.87 Å | Zn2+: 100.0%; ionic; 2.09 Å | |
| Zn-Asp914: 10.7%; ionic; 3.06 Å | ||
| Zn2+: 100.0%; ionic; 2.42 Å | ||
| Equilibration – heating 10 ps | Simulation 5000 ps – cooling 10 ps | |
| His782: 20.5%; H-acc; 2.77 Å | Gly955: 66.0%; H-acc; 2.69 Å | |
| Zn-Asp820: 94.2%; ionic; 2.10 Å | Zn-Asp820: 59.6%; ionic; 2.12 Å | |
| Zn-Asp820: 77.0%; ionic; 1.98 Å | Zn-His822: 10.1%; ionic; 3.15 Å | |
| Zn-Asp914: 17.1%; ionic; 1.90 Å | Zn-Asp914: 39.0%; ionic; 1.87 Å | |
| Zn2+: 100.0%; ionic; 2.06 Å | Zn-Asp914: 10.7%; ionic; 3.06 Å | |
| Zn2+: 100.0%; ionic; 2.42 Å | ||
| HDAC10–Aa50 | Before simulation (docking) | Initialization 100 ps (300 K) |
| Zn-Asp172: 79.5%; ionic; 2.03 Å | Zn-Asp172: 68.5%; ionic; 1.94 Å | |
| Zn-Asp172: 31.9%; ionic; 2.63 Å | Zn-Asp172: 54.6%; ionic; 1.87 Å | |
| Zn-His174: 11.7%; ionic; 2.92 Å | Zn-His174: 30.2%; ionic; 1.78 Å | |
| Zn-Asp265: 21.6%; ionic; 1.89 Å | Zn2+: 100.0%; ionic; 1.89 Å | |
| Zn2+: 100.0%; ionic; 2.22 Å | ||
| Equilibration – heating 10 ps | Simulation 5000 ps – cooling 10 ps | |
| Zn-Asp172: 13.3%; ionic; 1.75 Å | Zn-Asp172: 78.3%; ionic; 2.00 Å | |
| Zn-His174: 23.2%; ionic; 1.76 Å | Zn-Asp172: 62.3%; ionic; 2.47 Å | |
| Zn2+: 100.0%; ionic; 1.90 Å | Zn-His174: 40.0%; ionic; 1.95 Å | |
| Zn-Asp265: 20.4%; ionic; 1.92 Å | ||
| Zn2+: 100.0%; ionic; 2.25 Å | ||
4 Conclusion
In this study, novel potential inhibitors for cervical cancer targeting H. sapiens class II histone deacetylase (HDAC) were designed by de novo approach. Based on pharmacophore model, HDAC inhibitors (HDACIs) consist of three main parts: zinc binding group (ZBG), linker and hydrophobic cap (CAP). From this model, the compound series based on 4-[(2-oxo-1,3-thiazolidin-3-yl)carbonyl]aniline has been designed and modified. In silico studies have been performed for testing pharmacological properties and bioactivities of designed ligands.
The best ligands selection from molecular docking simulation was based on free energy of binding (ΔGbinding) and molecular interaction with enzyme active site. The selection yielded eight best ligands which have better binding affinity (lower free energy of binding) than SAHA, TSA, and VPA standards. Through interaction analysis of enzyme–ligand complex, all best ligands have interactions which followed hypothesis. Fragment 2-oxo-1,3-thiazolidine as zinc binding group (ZBG) forms a covalent coordination bond with the Zn2+ cofactor in HDAC charge-relay system with interaction that reached 100.0%.
Pharmacology analysis was then performed to assess bioactivity and pharmacological properties. All eight best ligands meet Lipinski’s rule of five and have drug likeness and drug score higher than standards. They also have a lower solubility in water than standards. Hence, it could be concluded that the best ligands are not easily excreted to enhance their bioactivity.
In silico preclinical trial was also conducted to evaluate the ADMET properties of the ligands including toxicological properties, passive absorption, health effect probability, MRDD prediction and oral bioavailability prediction. Based on the results, on average, all best ligands have a low toxicity and health effect probability on blood, cardiovascular system, gastrointestinal system, kidney, liver and lugs. They also have a high passive absorption and good oral bioavailability based on Veber’s and Egan’s rules.
Molecular dynamics (MD) simulations were also conducted to predict the stability of enzyme–ligand complexes on solvent effects. Based on this simulation, complexes formed by HDAC and corresponding best ligands have a good stability. All best ligands still strongly interacted with HDAC, mainly fragment 2-oxo-1,3-thiazolidine as zinc binding group (ZBG). The entire ligands also have a good synthetic accessibility predicted by CAESA. It could be inferred through these assessments that all designed ligands are potential to be novel HDAC inhibitors (HDACIs) and cervical cancer drug leads.
From the drug screening based on docking, dynamics and drug scan analysis, the best inhibitor ligands for HDAC4, HDAC5, HDAC6, HDAC7, HDAC9, and HDAC10 are F, Ib14, O38, Kb17, Gd40, and Aa50 respectively. The best ligands from drug screening then could be developed as novel potent drugs for cervical cancer targeting H. sapiens class II histone deacetylase (HDAC).
It can be deduced that the combination of in silico drug design based on de novo lead molecules design approach, molecular docking, pharmacology analysis, virtual screening, and manual interpretation can be more efficient to discover new effective inhibitors and drug leads that may have a great impact for future experimental studies. Computational or in silico methods play the important roles in drug discovery and development and can improve the efficiency of the novel HDACIs. This method may limit the use of chemicals for synthesis and biological testing, thus greatly reducing resource requirements than conventional methods. This study is expected to produce more potent HDAC inhibitors as drugs for cervical cancer treatment.
Acknowledgments
This research is supported by Hibah PUPT BOPTN Ditjen Dikti 2015 No. 0542/UN2.R12/HKP.05.00/2015. Usman Sumo Friend Tambunan and Arli Aditya Parikesit supervised the research and gave critical suggestions to the whole experiments; Abi Sofyan Ghifari and Cipta Prio Satryianto worked on the technical and experimental details, prepared the manuscript, and re-verified the data; Arli Aditya Parikesit added more advanced materials to the manuscript. All authors read and approved the final manuscript.
References
- Adv. Drug Deliv. Rev.. 2001;50(Suppl 1):S41-S67.
- Bioinformatics. 2006;22:195-201.
- Benigni, R., Bossa, C., Jeliazkova, N., Netzeva, T.I., Worth, A.P., 2008. EUR.
- Environ. Mol. Mutagen.. 2007;48:754-771.
- Bioinformatics. 2011;27:343-350.
- J. Biol. Chem.. 2008;283:26694-26704.
- J. Mol. Biol.. 2005;348:231-243.
- J. Cheminform.. 2009;1:8.
- Mol. Cancer Ther.. 2008;7:759-768.
- PLoS Comput. Biol.. 2012;8:e1002380.
- Int. J. Drug Des. Discov.. 2011;2:605-610.
- Int. J. Cancer. 2008;122:247-259.
- Success stories of computer-aided design. In: Ekins S., ed. Computer Applications in Pharmaceutical Research and Development. Hoboken, NJ, USA: John Wiley & Sons Inc; 2006.
- [Google Scholar]
- BMC Bioinform.. 2008;9:396.
- PLoS Biol.. 2006;4:e144.
- Clin. Cancer Res.. 2009;15:570-577.
- J. Med. Chem.. 2013;56:427-436.
- Best Pract. Res. Clin. Obstet. Gynaecol.. 2011;25:605-615.
- SAR QSAR Environ. Res.. 2008;19:413-431.
- Med. Chem. Res.. 2011;21:1156-1165.
- Drug Discov. Today. 2012;17:44-55.
- Clin. Proteomics. 2012;9:1.
- Comput. Biol. Med.. 2012;42:697-705.
- Molecules. 2002;7:51-62.
- J. Biol. Chem.. 2012;287:2317-2327.
- Int. J. Infect. Dis.. 2007;11(Suppl 2):S3-9.
- Exp. Mol. Pathol.. 2009;86:224-233.
- Hybrid drug delivery system for oropharyngeal, cervical and colorectal cancer – in vitro and in vivo evaluation. Saudi Pharm J.. 2013;21(2):177-186.
- [CrossRef] [Google Scholar]
- Eur. J. Med. Chem.. 2012;53:390-397.
- Bioorg. Med. Chem. Lett.. 2011;21:6767-6769.
- J. Biol. Chem.. 2008;283:11355-11363.
- Chem. Biol. Interface. 2011;1:209-227.
- Gynecol. Oncol.. 2007;107:S2-5.
- Semin. Cancer Biol.. 2009;19:144-152.
- In silico modification of suberoylanilide hydroxamic acid (SAHA) as potential inhibitor for class II histone deacetylase (HDAC) In: Ranganathan S., ed. BMC Bioinformatics. BioMed Central Ltd; 2011. p. :S23.
- [Google Scholar]
- Trends Bioinforma. 2012;5:25-46.
- HPV bioinformatics. In silico detection, drug design and prevention agent development. In: Rajkumar R., ed. Topics on Cervical Cancer with an Advocacy for Prevention. Rijeka, Croatia: Intech Publishing; 2012. p. :237-252.
- [Google Scholar]
- BMC Bioinform.. 2010;11(Suppl 7):S16.
- Public Health Genomics. 2009;12:291-307.
- Curr. Top. Med. Chem.. 2008;8:1555-1572.
- Nat. Rev. Drug Discov.. 2003;2:192-204.
- Gynecol. Oncol.. 2012;126:291-303.