Design of potent telomerase inhibitors using ligand-based approaches and molecular dynamics simulations studies

doi:10.21203/rs.3.rs-4029957/v1

Telomerase is a well-recognised and a promising target for cancer therapy. In this study, we selected ligand-based approaches to design telomerase inhibitors for the development of potent anticancer agents for future cancer therapy. Till date no telomerase inhibitors have been clinically introduced. To investigate the chemical characteristics required for telomerase inhibitory activity, a ligand-based pharmacophore model of oxadiazole derivatives reported from the available literature was generated using the Schrodinger phase tool. The generated pharmacophore model displayed five features, two hydrophobic and three aromatic rings. This selected pharmacophore hypothesis is validated by screening a dataset of reported oxadiazole derivatives. The pharmacophore model was selected for virtual screening using ZINCPharmer against the ZINC database. The ZINC database molecules with pharmacophoric features similar to the selected pharmacophore model and good fitness score were taken for molecular docking studies. With the pkCSM and SwissADME tools we predicted the pharmacokinetic and toxicity of top ten ZINC database compounds based on docking score, binding interactions and identified two in-silicopotential compounds with good ADME and less toxicity. Then both the hit molecules were exposed to molecular dynamic simulation integrated with MM-PBSA binding free energy calculations using GROMACS tools. The MM-PBSA calculations exhibited that the free binding energy of selected protein-ligand complexes were found stable and stabilized with nonpolar and van der walls free energies. Our study suggests that ZINC82107047 and ZINC8839196 can be used as hit molecules for future biological screening and for discovery of safe and potent drugs as telomerase inhibitors for cancer therapy.

Telomerase inhibitors

Pharmacophore

Molecular docking

Virtual screening

MD simulation.

Ligand based pharmacophore model of oxadiazole derivatives was developed and used for screening of newer oxadiazole.
Virtual screening of ZINC database using coordinates of validated pharmacophore hypothesis.
Molecular docking and ADMET predictions have identifiedtwo hit molecules.
Molecular Dynamics integrated with MM-PBSA free binding energy calculations identified ZINC82107047as potential hit compound.

Cancer is a deadly disease, which involves several events such as irregulated cell growth, proliferation of cancer cells from one organ to other organ and inhibition of apoptosis due to abnormal gene expression. Telomeres are present at each end of chromosome and composed of repeated sequence of hexanucleotide TTAGGG and six proteins are called shelterin [1]. Telomeres mainly protect the terminal end of linear chromosomes by forming a cap like structure. Each cell division shortens telomere length, after reaching the Hayflicklimit, the telomeres shorten and induce cell senescence or apoptosis, and eventually lead to cell death [2].In cancerous cells, the telomere length (TL) maintenance process is initiated, leading to further cell proliferation, which is a hallmark of cancer. Telomere length is restored by activation of telomerase enzyme (telomere elongation enzyme) and alternative lengthening of telomeres (ALT) [3, 4].

Telomerase enzyme comprises of telomerase reverse transcriptase unit (TERT), telomerase RNA template (TER) with its binding domain (TRBD), dyskerin, NHP2, nop10, gar1, reptin and pontin [5]. In most of primary cancer, the telomerase enzyme is activated which leads to uncontrolled cell replication [6].

Literature review revealed that oncogenic factors activate telomerase enzyme and enhances cancer cell proliferation by restoring telomeres length [7]. Therefore, telomerase is a major factor in differentiating normal healthy cells and cancer cells [8]. Despite all scientific efforts no telomerase inhibitors have been introduced, only one telomerase inhibitors Imetelstat, has progressed to clinical trials. Therefore, telomerase is become a promising target for the development of novel and potentially tumour specific anticancer chemotherapeutics with less toxicity and better pharmacokinetic properties. The in-silico drug design has increasingly become a practicable approach in chemical and biological sciences through establishment of a statistical relationship between molecular features and activities of assorted compounds in developing newer telomerase inhibitors.

In our efforts to explore telomerase inhibitors for anticancer activity, ligand-based pharmacophore models of oxadiazole derivatives as telomerase inhibitors were developed using PHASE module of Schrodinger. The combination of pharmacophore based virtual screening and molecular docking were used to design newer or potent derivatives as telomerase inhibitors [9]. The docking accuracy of compounds were determined by molecular dynamic (MD) simulations and Molecular Mechanics-Poisson–Boltzmann Surface Area (MM-PBSA) binding energy calculations using GROMACS. Drug-like property of top scored molecules was identified by absorption, distribution, metabolism, excretion, and toxicity (ADMET) studies [10]. The whole study is depicted in a flow chart (Fig. 1). A promising computational strategy was used in this in-silico work and can be used for the design of novel derivatives as telomerase inhibitors.

The molecular modelling studies such as pharmacophore, virtual screening, docking were implemented using Schrodinger software [11]. The identified hit molecules were further extended for molecular dynamics integrated with MM-PBSA calculations using GROMACS. The ADMET predictions were performed using free wares like pkCSM and Swiss ADME.

2.1. Preparation of dataset compounds

Initially, the available dataset of oxadiazole derivatives as telomerase inhibitors was collected from published articles to develop Pharmacophore models. Thirty-nine compounds listed in Table 1 were submitted to Phase module of Schrodinger to generate pharmacophore hypothesis by ligand-based approach [12-14]. The chemical structures of oxadiazole derivatives were drawn in Maestro 2D sketcher [15]. The Table 1 reports the chemical structure of the oxadiazole derivatives along with biological activity taken as test set for generating pharmacophore model. Six compounds were identified from the literature to validate the generated hypothesis [16,17]. The biological activity (IC₅₀inµmoles/litre) value of all derivatives was converted to a negative logarithm of IC₅₀ i.e., pIC₅₀.

2.2. Ligand preparation

All the dataset molecules were subjected to ligand preparation using "Ligprep" module of Maestro v2.5 [18]. The 2D structure was converted into 3D structure by using clean-up wizard. The process of ligand preparation consists of several steps i.e., addition of hydrogen atoms, removal of counter ions, generation of conformers and energy minimization of ligands with OPLS_2005 force field [19].

2.3. Pharmacophore hypothesis generation

The development of pharmacophore model based on the structural features of the ligand (Ligand-based) is one of the standard methods employed for the development of new hit molecules. In the Phase module, development of Common Pharmacophore Hypothesis (CPHs) was initiated with all the ligprep minimized molecules of selected dataset (“Schrödinger Release 2018-3: Phase, Schrödinger, LLC, New York, NY,” 2018). Pharmacophore hypothesis was run to find the common 4-5 pharmacophore features (default setting) that are similar in all ligands in spatial arrangement. Pharmacophoric features were generated for all the dataset molecules based upon the observed activity threshold of active and inactive molecules. As per the hypothesis, molecules having highest survival score were selected. Survival score is the quality of alignment of all active ligands [20]. The hypotheses was evaluated on the basis of the activity of ligand by using site, post-hoc, vector and volume scores and tabulated [21]. The vector score measures the angle formed by the two vector features such as acceptor, donor and aromatic ring aligned structure and a volume score measures the overlay of all ligands with the reference ligand [22]. Phase module of Schrodinger generates a pharmacophore hypothesis with 6 pharmacophore features based on the characteristic site such as A-hydrogen bond acceptor group, D-hydrogen bond donor group, H-hydrophobic group, R-ring aromaticity, P-positively ionizable group and N-negatively ionisable groups. Pharmacophore models were selected based on the alignment of pharmacophore features on the active ligand, with the maximum survival score.All the molecules were screened and parameters generated were recorded and tabulated.

2.4 Pharmacophore model validation

To validate the efficiency and selectivity of selected pharmacophore model, a dataset of oxadiazole molecules as telomerase inhibitors were screened on the generated pharmacophore hypothesis. The selected dataset included six molecules which were evaluated for telomerase enzyme inhibition assay by the same evaluation method. The Hypothetical screening of six compounds was correlated with experimental results for telomerase inhibitory activity.

2.5 Pharmacophore based virtual screening through ZINCPharmer webserver

Virtual screening of the ZINC database based on the pharmacophore features is a specific and appropriate method for the identification of new and potent lead molecules. For virtual screening we used the best pharmacophore hypothesis. The X, Y, Z co-ordinates value of all five pharmacophore features were put into ZINCPharmer which generated a ZINC database (ZINC database-zincpharmer.csb.pitt.edu). In the ZINC database screening, the screened molecules must match at least three pharmacophore features for the hypothesis with three or four pharmacophore features and for the hypothesis with five or more than five features, the screened molecules must match at least four features. The maestro v9.3 virtual screening workflow was employed for this study. Primarily, all the ZINC database compounds were subjected to ligprep module of Schrodinger. The ZINC database molecules with the best fitness scores were further screened using molecular docking studies in extra precision (XP) mode against telomerase to estimate ligand-protein binding interactions.

2.6 Molecular docking

To identify the possible interaction and conformation between the targeted protein and ZINC database, Grid-based Ligand Docking with Energetics (GLIDE) tool of Schrodinger Suite was used [23]. Crystallographic structure of telomerase protein was retrieved from PDB (Protein Data Bank ID: 5CQG). Among the list of telomerase protein codes, PDB: 5CQG containing co-crystallized ligand BIBR1532 was found to be a potent telomerase inhibitor and was selected for our studies. Protein was prepared by energy minimization and refinement with default settings of protein preparation wizard of Schrodinger. It includes removal of hetero group, unwanted chains and water molecules with less than 3 H-bonds to non-water molecules (Schrödinger 2018). Minimization of protein was specified at 0.30 Å [24]. For accurate docking a grid was generated at the active site residues of the protein which are responsible for binding interaction of the ligand to produce desired biological activity. Receptor grid was generated with default procedures using OPLS_2005 force field.

Molecular docking method used for screening of telomerase inhibitors was validated by re-docking the co-crystalized ligand BIBR1532 of the telomerase receptor (5CQG) [25]. Co-crystalized ligand BIBR1532 was separated from the receptor and submitted to the ligand preparation and run for Glide XP docking with default parameter. The XP output docked pose of BIBR1532 is superimposed over the reported binding pose of BIBR1532 with 5CQG and the calculated RMSD value.

Molecular docking of top ninety ZINC database molecules was carried out with telomerase protein (PDB ID: 5CQG) by using GlideXP (Extra Precision) scoring function with default parameters [26, 27]. The Glide XP docking predicts docking pose with higher accuracy and scores docking results based on the size, shape of ligand and receptor active site binding pocket and the binding affinity of the ligand with the protein. The binding interaction of ZINC molecules with amino acid residues of 5CQG were observed and all the docked compounds were ranked on the basis of docking score.

2.7. Predictions of Pharmacokinetic (ADME) and Toxicity (T) parameters

However, computational prediction of ADMET properties was not as specific as in vivo or in vitro evaluation, but it gave potential information to analyse the drug like properties of compounds. Pharmacokinetic property such as Absorption (A), Distribution (D), Metabolism (M), Excretion (E) and Toxicity (T) of ZINC database compounds were predicted by online screening tools such as pkCSM [28] (Product of University of Cambridge) and Swiss ADMET Prediction (Swiss Institute of Bioinformatics ©2018) [29]. The pkCSM is graph based structural signature tool, it maintains the balance between potency, safety and pharmacokinetic properties of the drug molecules. The SMILES format of ligand molecules was used for prediction of ADMET parameters.

Absorption parameters of orally administered drug molecules were calculated in terms of water solubility, Caco-2 permeability and intestinal absorption. Distribution of drug molecules was described by predicting the value of steady state volume of distributions (VDss), fraction of unbound drug (Fu) in plasma, permeability of Blood Brain Barrier (BBB) and Central Nervous System (CNS) distribution. Metabolic rate of the drug molecules affect efficacy and toxicity potential. Most of the drugs were metabolized by CytochromeP450 enzyme. The excretion parameter of drug molecules is expressed by calculating the value of total clearance and RENAL Organic Cation Transport2 (OCT2) substrate. Toxicity of drug molecules was measured in terms of AMES toxicity and oral rat acute toxicity (LD₅₀).

2.8. Molecular dynamics simulations and binding free energy analysis

Molecular dynamics simulation (MD) and binding free energy calculations of two ZINC database compounds ZINC82107047, ZINC8839196 and standard drug BIBR1532 were carried out using GROMACS version 2019 simulation package [30]. The system preparation was done by adding appropriate amount of sodium (Na⁺) and chlorine (Cl⁻) counter ions to neutralize. The missing hydrogen in the crystal structure was added. All the systems were solvated into an SCP water box from all directions [31-33]. Topology of all ligands was generated from PRODRUG web server [34] and the protein parameters were generated using groomos5a47 force field. System was first vacuum-minimized for 1500 steps using the steepest descent algorithm. Structures were then solvated with water extended simple point charge (SPCE) model in a cubic periodic box. Complex systems were further maintained in a suitable environment having a salt concentration of 0.15 M and energy minimization was carried out for 50,000 steps. The next stage was equilibration of system, which was done in two steps. First step was constant number of atoms, volume and temperature (NVT) equilibration for 1000ps (1ns) at 310 k; and the second constant number of atom, pressure and temperature (NPT) equilibration for 1000ps (1ns) at 1 bar pressure. Each resultant structure from the NPT equilibration phase was subjected to a final production run in the NPT ensemble for a simulation time of 50 ns. Root mean square deviation (RMSD) and root mean square fluctuation (RMSF) of the protein were calculated using gmxrms and gmxrmsf tools respectively [35]. The gmx gyrate and gmxsasa tools were used to calculate the radius of gyration (Rg) and solvent accessible surface area (SASA) respectively. The MM/PBSA approach was employed to understand the binding free energy (ΔG binding) of the inhibitors with the selected complex throughout the simulation time. The GROMACS utility g_mmpbsa was employed to estimate the binding free energy [36]. To obtain an accurate result, we computed ΔG for the last 20 ns with dt 1000 frames [37].

3.1. Ligand based pharmacophore model

Ligand based pharmacophore design for the oxadiazole derivatives accomplished pharmacophoric hypothesis containing five features HHRRR_1 two hydrophobic groups(H) and three aromatic rings (R) with highest survival score of 6.013 (Supplementary InformationS1). The pharmacophoric features of the best pharmacophore model were aligned on the most potent compounds of the selected data series, the distance between two pharmacophoric features and the angle was calculated which are shown in Fig. 2a and Fig. 2b. All the compounds were ranked in terms of fitness score, align score, vector score and volume score along with number of matched hypothesis features. Compounds with highest fitness score might have complementary features required for the potential inhibition of telomerase.

3.2. Pharmacophore model validation

The reliability of pharmacophore model HHRRR_1was checked by screening the model with validation dataset compounds. Screened dataset with aligned features and good fitness score are shown in Table 2. Compounds with more telomerase inhibitory activity showed more fitness score compared to less active compounds. Among the two compounds, one compound showed the alignment with all five HHRRR pharmacophoric features with a fitness score of 2.097 and showed the good correlation between the telomerase inhibitory activity and fitness score. Five compounds showed alignment with four HRRR pharmacophoric features with a fitness score in range 1.715 to 1.380. Compounds with four feature alignment showed more fitness score for highly active compounds and less fitness score for less active compounds. One compound showed the alignment with five HHRRR pharmacophoric features with good correlation between telomerase inhibitory activity and fitness score. Inactive compound has no common features like the developed pharmacophore model and therefore showed a less fitness score. These interpretations represent the reliability of developed pharmacophore hypothesis for further screening of database to obtain lead molecules as a potent telomerase inhibitor. All the compounds of validation set were subjected to molecular docking studies to explore the binding affinity with protein (Supplementary Information S2).

Table 2

Validation of pharmacophore model
S. No.	pIC₅₀	No. of sites matched	Fitness score	Align score	Vector score	Volume score
1.	5.638	5-HHRRR	2.097	0.512	0.815	0.708
2.	5.553	4-HRRR	1.715	0.787	0.834	0.618
3.	5.509	4-HRRR	1.628	0.975	0.873	0.609
4.	5.319	4-HRRR	1.452	0.798	0.610	0.587
5.	5.260	4-HRRR	1.416	1.246	0.861	0.586
6.	4.967	4-HRRR	1.380	1.291	0.862	0.579

3.2.1. Virtual screening of ZINC database through ZINCPharmer webserver

For the virtual screening of the ZINC database, pharmacophore hypothesis HHRRR_1was selected (Supplementary Information S3). Dataset molecules were obtained by overlapping their pharmacophoric (chemical group) features over corresponding pharmacophore features of HHRRR_1 pharmacophore hypothesis. Compound with good fitness score might have complementary characteristic required for potential inhibition of telomerase protein. The structure of the best hits molecules after virtual screening i.e., ZINC84512574, ZINC88339196, ZINC85237790, ZINC20540819 and ZINC74124901 are shown in Fig. 3.

Out of five features, four features were observed commonly in all the molecules. The fitness score of top 90 Zinc database molecules was observed in the range of 2.3–1.81. The align pose of few Zinc database compoundsZINC84512574, ZINC88339196, ZINC85237790, ZINC20540819 and ZINC74124901over the best HHRRR_1 pharmacophore hypothesis with good fitness score is shown in Fig. 4. After screening, 90 ZINC database compounds with required pharmacophoric features and good fitness score were subjected to XP docking (PDB: 5CQG).

3.2. Docking studies for hit identification of telomerase inhibitors

Docking methodology was validated by redocking of co-crystallized ligand BIBR1532 on to the telomerase receptor and superimposition of this docked pose was done over the binding pose of co-crystalized ligand. The RMSD value of superimposition of both the poses was found to be 0.9623 Å (Fig. 5). After validation of docking methodology, docking studies of top ninety ZINC database molecules was carried out to identify binding interaction with amino acid residues of telomerase receptor. All the docked poses were analysed to interpret the binding interaction with telomerase protein (PDB: 5CQG) and compared with co-crystalized ligand (BIBR1532).The XP score and binding interaction of top five scored ZINC database molecules is shown in Table 3, these molecules showed some common interactions which are hydrogen bonding of the ligands with ILE590, ARG486 and Pi-Pi stacking with TYR551, PHE494 amino acid residues of telomerase protein. These molecules showed similar binding interactions with the receptor relative to the binding interaction of the reference ligand. The binding affinity of the ligand with the targeted protein is calculated as that docking score.

From the obtained result, ZINC database compounds such as ZINC82107047, ZINC84512574 and ZINC88339196was the best docked inXP (-10, -9.2 and − 9.4 kcal/mol) docking mode Fig. 6. The docking score of hit molecules ZINC84512574 and ZINC88339196 obtained were near to co-crystalized ligand BIBR1532 (-6.9 kcal/mol).The compound ZINC82107047, ZINC84512574 and ZINC88339196 interacted with amino acid residue PHE494 and TYR551, via pi-pi stacking bonding. The binding pocket residue of ZINC88339196 is the same as obtained from the binding of BIBR1532 ligand (Fig. 6).

Table 3

Docking studies results of top five compounds from ZINC database
S. No.	Title	XP score	Binding interaction
S. No.	Title	XP score	Interaction residues	Type of interaction
Standard	BIBR1532	-6.9	PHE494	Pi-Pi stacking
1.	ZINC82107047	-10	ILE590	Hydrogen bond
	ZINC82107047		PHE494, TYR551	Pi-Pi stacking
2.	ZINC84512574	-9.2	TYR551	Pi-Pi stacking
3.	ZINC88339196	-9.4	TYR551, PHE494	Pi-Pi stacking
4.	ZINC85237790	-7.9	TYR551	Pi-Pi stacking
4.	ZINC85237790	-7.9	ARG486	Hydrogen bond
6.	ZINC20540819	-6.8	TYR551	Pi-Pi stacking

3.3. In silico Prediction of pharmacokinetic and toxicity parameters

Pharmacokinetic and toxicity parameters prediction play a major role in the process of drug discovery to identify potent and safe drug molecules which can be further taken for pre-clinical evaluation and lead optimization. pkCSM and Swiss ADME are open-source software, we used this two software to perform in silico prediction of ADMET.

3.3.1. pkCSM

The pharmacokinetic and toxicity prediction approaches are based on the concept of distance-based graph signatures of different physicochemical properties of the compound chemical structure. The ADMET properties have a significant role in discovery of drug molecules to identify the balance between pharmacokinetic properties, potency and safety of compounds which provide the information to proceed for clinical trials. The pkCSM toll based on the cut-off scanning concept and can be considered as a reliable tool for ADMET prediction. It builds thirty descriptors which are: 7 descriptors of absorption, 4 descriptors of distribution, 7 descriptors of metabolism, 2 descriptors of excretion and 10 descriptors of toxicity. The ADMET properties were calculated for the standard drug molecule and in silico active compounds identified from above computational approaches (Supplementary Information S4-S8).

Molecules screened for ADMET calculations were suggested to have more solubility in buffer and less water (aqueous) solubility. Pharmacophore screened ZINC database compounds ZINC82107047 and ZINC88339196showed good water solubility − 3.607 and − 3.845, Caco-2 permeability0.916 and 1.176 with human intestinal absorption − 2.752 and 95.908 percent respectively (Table 4). Compound ZINC82107047 and ZINC88339196 showed 0.076 and 0.108 unbound drugs in plasma. Drug metabolism affects the efficiency and toxicity of drug molecules. Chemical compounds are mainly metabolised by CYP450 enzyme. The pkCSM predicted the metabolism by different isoforms of CYP450 such as CYP2D6, CYP3A4, CYP1A2, CYP2C19, CYP2C9, CYP2D6, and CYP3A4. Compounds identified from ZINC database can be optimized for formulation design. Standard drug BIRB1532 showed 100 percent human intestinal absorption with − 3.456 water solubility and Caco-2 permeability 0.876 (Table 4). All the screened compounds should be effectively metabolized and eliminated from body after required therapeutic outcome is obtained. Safety profiles of compounds were determined by measuring total clearance of drug. The predictive toxicity results using the pkCSM free wares indicated the absence of mutagenicity or carcinogenicity in both ZINC database compounds (ZINC82107047 and ZINC88339196) in Ames test. The Oral Rat Acute Toxicity (LD₅₀) values for ZINC82107047 and ZINC88339196 were found to be 2.923 and 2.331 mol/kg respectively which comparatively more than the BIBR1532 (reference compound). Higher value of LD50 indicates less toxicity. Among, the dataset compounds ZINC82107047 and ZINC88339196demonstrated good clearance compared to standard drug (Table 4).

Table 4

ADMET parameters prediction of ZINC screened and standard telomerase inhibitors compounds using pkCSM
Category	Parameters	BIBR1532	ZINC82107047	ZINC88339196
Absorption	Water solubility (log mol/L)	-3.698	-3.607	-3.845
	Caco2 permeability (log Papp in 10^− 6 cm/s)	0.435	0.916	1.176
	Intestinal absorption (% human)	99.418	-2.752	95.908
Distribution	VDss (human) (log L/kg)	-2.027	0.135	0.127
	Fraction unbound (human)	0	0.076	0.108
	BBB permeability (log BB)	-0.217	-0.703	-0.911
	CNS permeability (log PS)	-1.786	-2.403	-2.262
Metabolism	CYP2D6 substrate	No	No	No
	CYP3A4 substrate	No	Yes	Yes
	CYP1A2 inhibitor	Yes	Yes	Yes
	CYP2C19 inhibitor	No	Yes	Yes
	CYP2C9 inhibitor	Yes	Yes	Yes
	CYP2D6 inhibitor	No	No	No
	CYP3A4 inhibitor	No	Yes	No
Excretion	Total Clearance (log ml/min/kg)	0.444	0.113	0.612
Excretion	Renal OCT2 substrate	No	No	No
Toxicity	AMES toxicity	No	No	No
Toxicity	Oral Rat Acute Toxicity (LD₅₀mol/kg)	2.173	2.923	2.331

Compounds have good intestinal permeability, good distribution with poor BBB permeability, metabolized by CYP enzymes with good clearance and less toxicity. Particularly, our hypothesis of study to investigate the potential novel molecules as telomerase inhibitor has shown impressive results and looking forward to screen for further in vitro and in vivoanticancer activity.

3.4. MD Simulations for in silico potential compounds

MD simulation studies were performed to identify the effects of protein (enzyme) structure changes and flexibility on complex interaction profile [38]. Along with the APO form of target protein, the complexes of two hit compounds ZINC82107047, ZINC8839196 (which exhibited favourable pharmacokinetic properties) and the standard drug (BIBR1532) were further analysed to study the dynamics and stability of protein–ligand complexes. In the present study, BIBR1532 was considered as the positive control drug and all the results of hit compound complexes were comparatively analysed. The stability of trajectory, flexibility, affinity of small molecules with receptor, and extent of compactness and folding behaviour were examined by analysing various structural parameters such as RMSD, RMSF, Rg, Hydrogen bonds and SASA between the target protein and respective ligand. The average values of structural parameters were recorded in Table 5.

Table 5

Average value of structural parameters from MD simulation
Structural parameter	APO	BIBR1532	ZINC82107047	ZINC8839196
RMSD (nm)	0.352075489	0.41171398	0.375250013	1.4844438
RMSF (nm)	0.194306711	0.233705034	0.197651678	0.622217785
Rg (nm)	2.805296041	2.839143211	2.841996133	4.125471412
SASA (nm²)	287.3564827	289.6046681	292.2441298	568.6321918

3.4.2. Root mean square deviation (RMSD)

RMSD provides an insight into whether or not the system has equilibrated and attained stability over the simulation period. The RMSD of the backbone atoms of the target protein was calculated and plotted as a function of time to analyze the structural stability of the receptor on binding to the ligand (Fig. 7). RMSD was calculated for protein backbone atoms for all the four complex structures that converged during the 50 ns MD simulation. The average RMSD values were calculated for the entire simulation trajectories. The average values of RMSD of APO form of 5CQG protein, BIBR1532 (standard drug), ZINC82107047 and ZINC8839196 were 0.35, 0.41, 0.37, and 1.48 nm respectively. Binding of standard drug BIBR1532 has slightly increased average RMSD value to 0.41. Interestingly, complex with ZINC82107047 has demonstrated with better RMSD value of 0.37. Such low RMSD value clearly point towards stability of all the three complexes [39]. From this data, it can be further inferred that the complexes BIBR1532 (standard drug) and ZINC82107047 were comparatively more stable.

3.4.3. Root mean square fluctuation (RMSF)

The Root Mean Square Fluctuation was comparatively analysed for the ligand-bound complexes (ZINC82107047 and ZINC8839196) along with the APO form (5CQG-APO) and standard drug BIBR1532 complexes to examine the average residual fluctuations, motion and flexibility of amino acid residues of target protein on binding to ligands during the simulation time. The mobility and flexibility of receptor-ligand complex and receptor APO form were represented by average RMSF values for Cα atoms of the protein (Fig. 8) [39]. Average RMSF values were recorded in Table 6. Here, binding of standard drug has shown slight increase in fluctuations with average RMSF value of 0.233. Few residues in the range of 150–200 have exhibited fluctuations above 1 nm and other fluctuations were less than 0.5 nm and stable. Similar to RMSD values, fluctuations were very minimal for ZINC82107047 complex with value of 0.197. For ZINC8839196 complex, fluctuations were high compared to other complexes but found less than 1.25 nm. Overall, RMSF values indicated that ZINC82107047 complex was found stable.

3.4.4. Radius of gyration

The radius of gyration (Rg) is the root mean square distance of atoms from their rotational axis [40]. It is the structural parameter that gives an insight about the compactness, rigidity and folding behaviour of the receptor and its change with time during simulation. Lower and constant Rg values indicate compactness and a stably folded nature while high and jerky fluctuations reveal instability in folding. The Rg values for backbone atoms of the target protein were calculated and plotted against the simulation time (represented in Fig. 9). For the APO form, the Rg showed slight fluctuations till 10000 ps with the values ranging between 2.9 and 2.8 nm. Subsequently, a decline was observed and then, the values did not fluctuate till the end of simulation and exhibited a constant value of 2.7 nm; this indicates the stably folded nature of 5CQG-APO. In the case of standard, even though fluctuation were observed throughout the simulation, the variations were moderate (2.92–2.79 nm), which is considerable. These fluctuations neither affected the regular folding pattern nor the steady binding of the ligand. Further, In the case of ZINC82107047, the variation pattern of Rg values was similar to standard drug. The variation of backbone gyrate value was in the moderate range 2.72–2.82 nm. In case of ZINC8839196, more fluctuation was observed during the simulation. The overall analysis clearly signifies that the receptor attained a compact state during simulation and there were no abrupt fluctuations, which indicate the stably folded nature of the protein on binding to the ZINC82107047.

3.4.5. Solvent accessible surface area (SASA)

The solvent accessible surface area analysis (SASA) was performed to understand the solvent behaviour of the target protein on binding to small molecules and was compared with surface area changes of the APO protein (5CQG-APO). The binding of ligands to receptors definitely induces structural and conformational changes, leading to variations in the protein volume; indirectly, this gives an insight about stability of the complex during simulation. The SASA values were contributed by the hydrophobic residues and their exposure from the hydrophobic core region leads to decompression of the receptor, resulting in instability [41, 42]. Lower and minimal fluctuation in values were expected, which signifies the stabilization, compression and folding of the target protein during simulation. The SASA values were calculated for the entire simulation and plotted against time, as shown in Fig. 10.

The APO protein showed moderate fluctuations (gradual decrease) initially, followed by immediate stabilization and this stabilized value persisted till termination. At the end, the protein surface shrunk slightly compared to its native state and the values were in the range of 275–293 nm². In the case of 5CQG in complex with standard drug and ZINC82107047, the protein surface area fluctuations were less and the variation pattern was similar in both complexes. Variation in both standard drug and ZINC82107047 was 323 − 275 nm² and 326 − 280 nm² respectively, indicating that the complexes were stable during simulation. The deflection of values in ZINC8839196 was in the range of 626 − 540 nm². The surface area of the receptor protein withZINC8839196 was higher than the APO form and this might be due to structural changes on binding of ligand in the binding site.

3.4.6. Hydrogen bonds

To examine the binding affinity of ligands with the target protein, MD trajectories were analysed to interpret the extent of hydrogen bond formation during the entire simulation, as shown in Fig. 11. The standard drug complex formed a reasonable number of H-bonds with the receptor protein with two hydrogen bonds and three bonds at very intervals. Interestingly, ZINC82107047 complex has formed a maximum of four bonds at several time frames and five to six bonds at very few frames indicating stronger affinity towards the target. For the ZINC8839196 complex, the ligand formed one to three hydrogen bonds throughout the simulation at few time intervals. The ZINC82107047 complex stabilized and this can be distinctly inferred through other structural parameters too.

3.5. MM-PBSA: binding free energy calculations

The Molecular Mechanics-Poisson–Boltzmann Surface Area continuum solvation is a widely accepted method for estimation of the total binding free energy of ligands complex with biological macromolecules [43]. The g_mmpbsa tool was used to calculate the binding free energy for complexes BIBR1532, 5CQG-ZINC82107047 and ZINC8839196 on binding to 5CQG for the entire simulation (50 ns) by importing and analyzing the respective MD trajectories. The binding free energy is the gross summation of several energy components such as polar solvation energy, SASA, non-polar solvation energy and the non-bonded interaction energies such as van der Waals and electrostatic energies [40]. All the calculated energy values were listed in Table 6.

All the three complexes, i.e., BIBR1532, ZINC82107047 and ZINC8839196 exhibited a high mean negative value, ΔGbind of -188.812 ± 24.371 kJ/mol, -161.878 ± 15.251 kJ/mol and − 203.734 ± 16.501 kJ/mol, indicating stronger interactions between the ligands and receptor. Further, the energy terms that sum up to give the binding free energy were disclosed for all the complexes. Among them, the van der Waals energy was the chief motive for the strong binding of ligands. The van der Waals energy exhibited by complexes 5CQG-BIBR1532, ZINC82107047, 5CQG-ZINC8839196 were − 234.136 ± 21.478 kJ/mol, -219.956 ± 17.706 kJ/mol and − 212.525 ± 17.272 kJ/mol respectively, indicating stronger intermolecular interactions between the ligands and protein thereby strengthening the binding affinity. However, SASA energy contributed equally more to the binding energy along electrostatic energy in all three complexes. The polar solvation energy is the only component that contributes positively to the binding energy and it was found higher with ZINC82107047 complex. Surprisingly, ZINC82107047 complex was the only one found stabilized better with electrostatic energy than SASA energy. The complete estimation analysis distinctly revealed that both the hit compounds ZINC82107047 has strong binding affinity towards the target 5CQG and was at par and better with the positive control. Further experimental investigation might help in understanding better about the structural features.

Table 6

MM-PBSA energy values of respective complexes from GROMACS
Energy terms in KJ/mol	BIBR1532	ZINC82107047	ZINC8839196
van der Waals	-234.136 ± 21.478	-219.956 ± 17.706	-212.525 ± 17.272
Electrostatic	-18.451 ± 7.786	-45.255 ± 10.131	-13.301 ± 5.531
Polar solvation	80.915 ± 14.127	121.360 ± 20.614	38.714 ± 8.157
SASA energy	-17.140 ± 1.136	-18.027 ± 1.220	-16.622 ± 1.104
Binding	-188.812 ± 24.371	-161.878 ± 15.251	-203.734 ± 16.501

Interestingly, both compound ZINC82107047 and ZINC8839196 share some common chemical features. Two five membered heterocyclic rings joined with propyl bridge and both the heterocyclic rings are substituted with aliphatic, aromatic and heterocyclic rings. First and second compounds have bicyclic rings as substitution. Study concludes that ZINC82107047 can be considered as hit compound and can be taken up for synthesis and experimental investigation as telomerase inhibitor to correlate the simulation results. Further computational and experimental studies can be carried out to understand the structure activity relationship (SAR) on length of alkyl chain between two oxadiazole rings and substitution of different aromatic and heterocyclic rings on to oxadiazole ring.

In this study we described the design of novel telomerase inhibitors with different in silico molecular modelling techniques.The ligand-based pharmacophore hypothesis of telomerase inhibitors was generated using Phase module of Maestro v9.3. The generated pharmacophore hypothesis HHRRR_1 with two hydrophobic and three aromatic ring pharmacophoric features were screened by oxadiazole dataset and the screened molecules exhibited good fitness score. The HHRRR_1 pharmacophore model was utilized for virtual screening of the ZINC database andtop molecules were selected. Pharmacokinetic and toxicity of hit molecules were predicted using pkCSM and SwissADME free wares for further refining of hit molecules. Additionally, MD simulation studyof ZINC database identified ZINC82107047 and ZINC8839196 as potential hit compounds for further mechanistic studies as telomerase inhibitors for cancer therapy. MM-PBSA rescoring method was used to calculate binding free energy of two complexes and standard drug-receptor complex. This study concludes that all the techniques used in this work are capable to predict hit molecules as telomerase inhibitors. It further suggests carrying out in vitro and in vivo biological evaluation of the predicted molecules as potential and safe telomerase inhibitors for future drug discovery.

Acknowledgments

The authors gratefully acknowledge Department of Science and Technology (DST), New Delhi for providing Woman Scientist (WOS-A) Project to Miss. Shalini Bajaj (SR/WOS-A/CS-98/2016), DST-SERB (CRG/2019/001452) for financial assistance to YCM and SVKM’s NMIMS SPPSPTM, School of Pharmacy and Technology Management, Mumbai and Faculty of Pharmacy, M.S. Ramaiah University of Applied Sciences, Bangalore for providing the necessary facilities to carry out this research work.

Funding: Department of Science and Technology (DST), New Delhi, under the scheme of WOS-A Project to Miss. Shalini Bajaj (SR/WOS-A/CS-98/2016) and DST-SERB (CRG/2019/001452) for financial assistance to YCM and SPPSPTM, School of Pharmacy and Technology Management, SVKM’s NMIMS, Mumbai.

Conflicts of interest/Competing interests: The authors declare no conflict of interests.

Availability of data and material: Data generated is provided in manuscript and supplementary information.

Code availability: Software programs used are cited at appropriate places. No codes are used to perform the study.

Authors' contributions: SB and MM designed the strategy and performed the computational work and wrote the manuscript. MYC supervised the project and contributed to the final version of the manuscript.

O’Sullivan RJ, Karlseder J (2010) Telomeres: protecting chromosomes against genome instability. Nat Rev Mol Cell Biol 11:171-181.
Hayflick L, Moorhead PS (1961) The serial cultivation of human diploid cell strains. Exp Cell Res 25:585-621
Nault JC et al. (2013) High frequency of telomerase reverse-transcriptase promoter somatic mutations in hepatocellular carcinoma and preneoplastic lesions. Nat Commun 4:2218.
Dilley RL, Greenberg RA (2015) Alternative telomere maintenance and cancer. Trends Cancer 1:145-156.
Dey A, Chakrabarti K (2018) Current perspectives of telomerase structure and function in eukaryotes with emerging views on telomerase in human parasites. Int J Mol Sci 19(2):pii:E333
Kim NW, Piatyszek MA, Prowse KR, Harley CB, West MD, Ho PL, et al. (1994) Specific association of human telomerase activity with immortal cells cancer. Sci 266(5193):2011–2015.
Ozturk MB et al. (2017) Current insights to regulation and role of telomerase in human diseases. Antioxidants 6:1-13.
Felsher DW et al. (2017) Oncogenes and the initiation and maintenance of tumorigenesis. In Molecular Basis of Human Cancer (2^nd ed) (Coleman, W.B. and Tsongalis, G.J., eds) 143–145, Springer.
Koes DR, Camacho, C. J., ZINCPharmer: pharmacophore search of the ZINC database. Nucleic Acids Research, 40, W409-W414.
Pires DEV, Blundell TL, Ascher D.B. (2015). pkCSM: Predicting small-molecule pharmacokinetic and toxicity properties using graph-based signatures. J Med Chem 58(9):4066-4072.
Schrödinger Release 2018-3: Maestro, Schrödinger, LLC, New York, NY, 2018. (2018).
Zhang F, Wang XL, Shi J, Wang SF, Yin Y, Yang YS, Zhang WM, Zhu HL (2014) Synthesis, molecular modeling and biological evaluation of N-benzylidene-2-((5-(pyridin-4-yl)-1,3,4-oxadiazol-2-yl)thio) acetohydrazide derivatives as potential anticancer agents. Bioorg Med Chem 22:468-477.
Zhang XM, Qiu M, Sun J, Zhang YB, Yang YS, Wang XL, Tang JF, Zhu HL (2011) Synthesis, biological evaluation, and molecular docking studies of 1,3,4-oxadiazole derivatives possessing 1,4-benzodioxan moiety as potential anticancer agents. Bioorg Med Chem 19:6518-6524
Zhang YB, Wang XL, Liu W, Yang YS, Tang JF, Zhu HL (2012) Design, synthesis and biological evaluation of heterocyclic azoles derivatives containing pyrazine moiety as potential telomerase inhibitors. Bioorg Med Chem 20:6356-6365.
Schrödinger Release 2018-3: LigPrep, Schrödinger, LLC, New York, NY. (2018).
Zheng QZ, Zhang XM, Xu Y, Cheng K, Jiao QC, Zhum HL (2010) Synthesis, biological evaluation, and molecular docking studies of 2-chloropyridine derivatives possessing 1,3,4-oxadiazole moiety as potential antitumor agents. Bioorg med chem 18:7836-7841.
Sun J, Zhu H, Yang ZM, Zhu HL (2013) Synthesis, molecular modeling and biological evaluation of 2-aminomethyl-5-(quinolin-2-yl)-1,3,4-oxadiazole-2(3H)-thione quinolone derivatives as novel anticancer agent. Euro J Med Chem 6023-28.
Dixon SL, Smondyrev AM, Knoll EH, Rao SN, Shaw DE, Friesner RA (2006) PHASE: a new engine for pharmacophore perception, 3D QSAR model development, and 3D database screening. Methodology and preliminary results. J Computer-Aided Molecular Design 20:647.
Kleinjung J, Fraternali F (2014) Design and application of implicit solvent models in biomolecular simulations. Current Opinion in Stru Bio 25:126.
Friesner RA, Banks JL, Murphy RB et al. (2004) Glide: a new approach for rapid, accurate docking and scoring. 1. Method and assessment of docking accuracy. J Med Chem 47:1739-1749.
Sallam AA, Houssen WE, Gissendanner CR, Orabi KY, Foudah AI, El Sayed KA (2013) Bioguided discovery and pharmacophore modeling of the mycotoxic indole diterpene alkaloids penitrems as breast cancer proliferation, migration, and invasion inhibitors. Med Chem Comm 4 (10):1360-1369.
Hall MD, Salam NK, Hellawell JL, Fales HM, Kensler CB, Ludwig JA, Szakáes G, Hibbs DE, Gottesman MM (2009) Synthesis, activity, and pharmacophore development for isatinbeta-thiosemicarbazones with selective activity toward multidrug-resistant cells. J Med Chem 52:3191-3204.
Murahari M, Prakash KV, Peters GJ, Mayur YC (2017) Acridone-pyrimidine hybrids- design, synthesis, cytotoxicity studies in resistant and sensitive cancer cells and molecular docking studies. Euro J Med Chem 139:961-981.
Schrödinger Release 2018-3: Glide, Schrödinger, LLC, New York, NY, 2018. (2018).
Bryan C, Rice C, Hoffman H, Harkisheimer M, Sweeney M., and Skordalakes, E. (2015) Structural Basis of Telomerase Inhibition by the Highly Speciﬁc BIBR1532. Structure 23:1934-1942,
Glide ver. 5.9, Schrödinger LLC, New York, NY (USA), 2010
Friesner RA, Murphy RB, Repasky MP, Frye LL, Greenwood JR, Halgren TA, Sanschagrin PC, Mainz DT (2006) Extra precision Glide: Docking and scoring incorporating a model of hydrophobic enclosure for protein-ligand complexes. J Med Chem 49:6177-6196.
Cutinho PF, Roy J, Anand A, Shankar R, Murahari M, Venkataramana CHS (2019) Design of metronidazole derivatives and flavonoids as potential non-nucleoside reverse transcriptase inhibitors using combined ligand and structure-based approaches. J BiomoleStruc Dynamics 1626-1648.https://doi.org/10.1080/07391102.2019.1614094.
Daina A, Michielin O, Zoete V (2017) SwissADME: a free web tool to evaluate pharmacokinetics, drug- likeness and medicinal chemistry friendliness of small molecules. Nature Publishing Group 1-13. https://doi.org/10.1038/srep42717.
Lindahl E., Abraham M.J., Berk H., Van Der Spoel D., GROMACS 2019.4 manual, GROMACS Doc. (2019).
Gangadharappa BS, Sharath R, Revanasiddappa PD, Chandramohan V, Balasubramaniam M, Vardhineni TP (2020) Structural insights of metallo-betalactamase revealed an effective way of inhibition of enzyme by natural inhibitors. J Biomol StructDyn 38: 3757-3771, https://doi.org/10.1080/07391102.2019.1667265.
Kumar B, Parasuraman P, Murthy TPK, Murahari M, Chandramohan V (2021) In silico screening of therapeutic potentials from Strychnosnux-vomica against the dimeric main protease (Mpro) structure of SARS-CoV-2. JBiomol StructDyn1-19.
Krishna S, Kumar SB, Murthy TPK, Murahari M (2021) Structure-based design approach of potential Bcl-2 inhibitors for cancer chemotherapy. Comput Bio Med 104455.
Schüttelkopf AW, Van Aalten DMF (2004) PRODRG: a tool for high-throughput crystallography of protein-ligand complexes, Acta Crystallogr Sect D BiolCrystallogr 60: 1355-1363, https://doi.org/10.1107/ S0907444904011679.
Thangavel M, Chandramohan V, Shankaraiah LH, Jayaraj RL, Poomani K, Magudeeswaran S, Govindasamy H, Vijayakumar R, Rangasamy B, Dharmar M, Namasivayam E (2020) Design and molecular dynamic investigations of 7,8-dihydroxyflavone derivatives as potential neuroprotective agents against alpha-synuclein. Sci Rep 10:1-10, https://doi.org/10.1038/s41598-020-57417-9.
Kumari R, Kumar R, Lynn A (2014) G-mmpbsa -A GROMACS tool for high-throughput MM-PBSA calculations. J Chem Inf Model 54:1951-1962, https://doi.org/10.1021/ci500020m.
Prasanth DSNBK, Murahari M, Chandramohan V, Panda SP, Atmakuri LR, Guntupalli C (2020) In silico identification of potential inhibitors from Cinnamon against main protease and spike glycoprotein of SARS CoV-2. J Biomol Struct Dyn 1-15, https://doi.org/10.1080/07391102.2020.1779129.
Aliebrahimi S, MontasserKouhsari S, Ostad SN, Arab SS, Karami L (2018) Identification of phytochemicals targeting c-met kinase domain using consensus docking and molecular dynamics simulation studies. Cell Biochem Biophys 76:135-145, https://doi.org/10.1007/s12013-017-0821-6.
Ghosh R, Chakraborty A, Biswas A, Chowdhuri S (2020) Evaluation of green tea polyphenols as novel corona virus (SARS CoV-2) main protease (Mpro) inhibitors–an in silico docking and molecular dynamics simulation study. J Biomol Struct Dyn1-13, https://doi.org/10.1080/ 07391102.2020.1779818.
Reed JC, ZhaH, Aime-Sempe C, Takayama S, Wang HG (1996) Structure-function analysis of bcl-2 family proteins: regulators of programmed cell death. Adv. Exp.Med. Biol 406:99-112, https://doi.org/10.1007/978-1-4899-0274-0_10.
Garg S, Anand A, Lamba Y, Roy A (2020) Molecular docking analysis of selected phytochemicals against SARS-CoV-2 Mpro receptor. Vegetos 33 (4):766-781, https://doi.org/10.1007/s42535-020-00162-1.
Udhaya Kumar S, Thirumal Kumar D, Mandal PD, Sankar S, Haldar R, Kamaraj B, Walter C, Jebaraj E, Siva R, George Priya Doss C, Zayed H (2020) Comprehensive in silico screening and molecular dynamics studies of missense mutations in Sjogren-Larsson syndrome associated with the ALDH3A2 gene. Adv Protein Chem Struct Biol 120:349-377, https://doi.org/10.1016/bs.apcsb.2019.11.004.
Genheden S, Ryde U (2015) The MM/PBSA and MM/GBSA methods to estimate ligand binding affinities. Expet Opin Drug Discov 10 (5):449-461, https://doi.org/10.1517/17460441.2015.1032936.
Mukherjee S, Dasgupta S, Adhikary T, Adhikari U, Panja SS (2020) Structural insight to hydroxychloroquine-3C-like proteinase complexation from SARS-CoV-2: inhibitor modelling study through molecular docking and MD-simulation study. J Biomol Struct Dyn 1-13, https://doi.org/10.1080/07391102.2020.1804458

No competing interests reported.

Design of potent telomerase inhibitors using ligand-based approaches and molecular dynamics simulations studies

Status:

Version 1

Abstract

Figures

Highlights

1. Introduction

2. Material and methods

3. Result and Discussion