E-pharmacophore filtering and molecular dynamics simulation studies in the discovery of potent drug-like molecules for chronic kidney disease

Chronic kidney disease (CKD) is a prominent health issue reported globally. The level of the vitamin D receptor (VDR) and cytochrome P450 enzyme 24-hydroxylase (CYP24A1) are crucial in the pathogenesis of secondary hyperparathyroidism (sHPT) in CKD. An elevated expression of the CYP24A1 leads to the deficiency of vitamin D and resistance to vitamin D therapy. Hence, VDR agonists and CYP24A1 antagonists are suggested to CKD patients for the management of biochemical complications. CTA-018 is a recently reported analog and acts as a potent CYP24A1 inhibitor. It inhibits CYP24A1 with an IC50 27 ± 6 nM, about 10 times more potentially than the non-selective inhibitor ketoconazole (253 ± 20 nM), and it is also been reported to induce the VDR expression. Thus, CTA-018 is under clinical trial among CKD patients. In this study, combined molecular docking and pharmacophore filtering were employed to identify compounds better than CTA-018. A huge set of 9127 compounds from Sweet Lead database were docked into the active site of VDR using Glide XP program. E-pharmacophore was developed from both the targets along with CTA-018. The compounds retrieved from the two different pharmacophore-based screening were re-docked into the active site of CYP24A1. The hits that bind well at both the active sites and matched with the pharmacophore models were considered as possible dual functional molecules against VDR and CYP24A1. Further, molecular dynamics simulation and subsequent energy decomposition analyses were also performed to study the role of specific amino acids in the active site of both VDR and CYP24A1.


Introduction
Chronic kidney disease (CKD) is a global public health problem, which affects 10% of the world population (Couser, Remuzzi, Mendis, & Tonelli, 2011). An impaired Vitamin D metabolism has been commonly observed among CKD patients (Cozzolino et al., 2006). Vitamin D is either synthesized in the skin or ingested through diet and it is transported to the liver in the form of 25-hydroxyvitamin D. This hydroxyvitamin D further hydroxylated by the enzyme 1α-hydroxylase in the kidney to yield 1α,25-dihydroxyvitamin D, which is the active form of vitamin D (Zhang et al., 2013). This metabolite is responsible for the effects of vitamin D on calcium and phosphorus metabolism, bone health, and the regulation of parathyroid function (Deluca, 2014). Vitamin D insufficiency is the major cause of CKD and hyperphosphatemia due to the phosphate elevation.
In recent years, different vitamin D analogs are used to mimic the function of 1α, 25-dihydroxyvitamin D (Hruska, Mathew, Lund, Qiu, & Pratt, 2008). But most of the drugs are ineffective because of protein-ligand binding affinity, disassociation of ligand and receptors, different metabolism of ligands (Adams & Hewison, 2008). Most of the drugs are designed against VDR drug target. The VDR therapies are controlled by the Cytochrome P450 enzyme 24-hydroxylase (CYP24A1)metabolizing agent for 1α, 25(OH) 2 D 3 , 25(OH)D 3 , and administered analogs by hydroxylation reaction (Chen, Sakaki, Yamamoto, & Kittaka, 2012).
Recent studies have explained that the elevation of phosphate leads to CYP24A1 over expression and subsequent vitamin D metabolism dysfunction (Zhang et al., 2013). Reports from uremic rats and humans have suggested that CYP24A1 expression causes vitamin D metabolism dysfunction. CYP24A1 has been known as a potent drug target for CKD (Petkovich & Jones, 2011). Posner et al. (2010) proposed a dual active compound named sulfone GHP-GH-16, 23-diene-25S02-1 (CTA018/MT2832). The CTA-018 compound mediates the VDR and inhibits the over expression of CYP24A1. This compound has a high binding affinity with VDR and also with CYP24A1. Persistent with these findings, it is also proclaimed that, CTA-018 is a selective drug candidate to treat CKD and hyperphosphatemia and for those who are undergoing regular hemodialysis (Posner et al., 2010).
In the VDR, two amino acids, namely His305 and His397 (His301 and His393 in rat) control the agonist and antagonist mechanism of VDR (Kudo et al., 2014). The activity of the VDR is based on these two amino acids. The CTA-018 was designed on the basis the selective binding mechanism of VDR and CYP24A1 (Posner et al., 2010).
Drug design and development are considered as a time and resource consuming process which involves various stages of screening. Computational chemistry is an ever growing field with the combination of computational power, chemical, and biological space that streamlines drug discovery, development, and organization. Structurebased drug design and ligand-based drug design are the two most important and commonly used computational approaches (Kumar, Hendriks, Janes, de Graaf, & Lauffenburger, 2006). Even though structure-based drug design holds a few drawbacks, it assists the pharmaceutical researchers in the development of new drugs in minimal time and significant level of confidence (Dror, Shulman-Peleg, Nussinov, & Wolfson, 2004). The molecular docking study is one of the successive ways to predict the binding orientation of the substrate with its receptor. The main aim of the docking study is to identify the appropriate binding pocket of protein, which can accommodate the correct conformations of ligands. It also calculates the binding affinities of a protein and ligand complex. More reliable protein-ligand complexes are predicted by fast and inexpensive docking protocols with molecular simulation techniques (Kitchen, Decornez, Furr, & Bajorath, 2004). The molecular docking technique includes an exhaustive conformational space search in a limited time, whereas the molecular dynamics (MD) method allows both the protein and ligand to be flexible and this induces the newly introduced ligands around the receptor binding site (Kitchen et al., 2004;Lin, Perryman, Schames, & McCammon, 2002). MD simulation can be used in various stages in computational chemistry viz., protein preparation before docking in order to optimize for structure flexibility, effect of solvent among proteinligand complex, calculation of binding free energies followed by an accurate ranking of the potential ligands (Huo, Wang, Cieplak, Kollman, & Kuntz, 2002;Reddy et al., 2013;Schames et al., 2004).
In this work, we attempted to find more potent dual active compound than CTA-018 using different in silico approaches such as molecular docking, E-pharmacophore generation, and MD simulation studies. Initially, the CTA-018 molecule was docked in the active site of both VDR and CYP24A1. Further, two different e-Pharmacophore models were generated using VDR/CTA-018 and CYP24A1/CTA-018 complexes. The Sweet Lead database compounds were screened to fit well in the active site of VDR. The screened compounds have undergone for the pharmacophore-based screening against VDR and CYP24A1 e-pharmacophores. Further, pharmacophore-based screening was performed against the Sweet Lead database which contained approved drugs and herbal medicines. Four different compounds were identified based on pharmacophore-based screening and molecular docking. Further, the contribution of individual amino acid-ligand interaction energies in total binding energy, i.e. amino acid decomposition analysis (ADA) (Hensen et al., 2004;Thangapandian, John, & Lee, 2012), was estimated for all the identified drug molecules. ADA was calculated on the basis of receptor structure and could be applied to different types of scaffolds. ADA efficiently predicted the effect of individual residues on ligand-receptor interactions and this information can be used as supporting information in the drug discovery process. Based on this available information, estimation of the optimum binding geometry could be more useful in choosing the best fragment(s) which leads to the improvement of ligand potency profiles. Finally, the conformational changes of these lead compounds in the active sites of VDR and CYP24A1 were analyzed by MD simulation studies. The findings from this study may provide meaningful clues for further designing of highly selective drug molecules to treat osteoporotic conditions among CKD patients.

Materials and methods
All the computational analyses were carried out on Centos 6.2 Linux platform in Intel core 2 duo processor RAM 4 GB.

Molecular docking studies
Since the unavailability of human CYP24A1 protein structure, the rat CYP24A1 was selected for the study. The crystal coordinates of rat VDR protein in complex with partial agonist 26-adamantyl-23-yne-19-norvitammin D were downloaded from PDB database. The sequence identity of human and rat VDR proteins shows 94% similarity in its sequence and the active site regions are similar in both the species. Superimposition of the human and rat VDR proteins are shown in Supplementary Figure 1. For the consistency of the experiment, we have selected both the proteins from rat species for this study. The VDR (pdb id -3VTC) (Kudo et al., 2014) and CYP24A1 (pdb id -3K9V) (Annalora et al., 2010) X-ray crystallographic structures were downloaded from the Protein Data Bank (Berman et al., 2000). Initially, the protein structures were prepared by protein preparation wizard module of Schrodinger software (Maestro, Schrodinger, LLC, New York, NY, 2015). The energy was minimized using OPLS_2005 force field (Jorgensen & Tirado-Rives, 1988). Glide docking module of the Schrodinger was used to simulate the binding mode of CTA-018 to VDR and CYP24A1 (Friesner et al., 2004). The CTA-018 molecule was prepared with the help of LigPrep module in Schrodinger (LigPrep, Schrödinger, LLC, New York, NY, 2015). This module generates the proper ionization, tautomers, stereochemistries, and ring conformations for the input molecule (CTA-018). Partial atomic charges were calculated using OPLS force field.
The CTA-018 molecule was docked into the catalytic site of both the proteins (i.e. VDR and CYP24A1) using Glide extra-precision (XP) mode with default protocols. Initially, the active site was defined within 5Å surrounding of the co-crystallized ligand. Further the grid was generated around the catalytic sites. The docking structures which correspond to the top ranking poses were selected for the e-pharmacophore generation (Phase, Schrödinger, LLC, New York, NY, 2015).

Energy-optimized pharmacophore mapping
Energy-optimized pharmacophore (E-pharmacophore) was obtained by mapping the energetic terms of the Glide XP scoring function onto atom centers. The Glide XP scoring terms were computed, and the energies were mapped onto atoms. Next, pharmacophoric sites were generated, and the Glide XP energies of the atoms that comprised each pharmacophore site were summed. The sites were then ranked based on these energies, and the most favorable sites were selected for the generation of pharmacophore hypothesis. These pharmacophores were then used as queries for pharmacophore-based screening. Pharmacophore sites were automatically generated from the docked complex using Phase module (Phase, Schrödinger, LLC, New York, NY, 2015) with the default set of six chemical features: hydrogen bond acceptor (A), hydrogen bond donor (D), hydrophobic regions (H), negatively charged group (N), positively charged group (P), and aromatic ring (R). Hydrogen bond acceptor sites were represented as vectors along with the hydrogen bond axis in accordance with the hybridization of the acceptor atom. Hydrogen bond donors were represented as projected points, located at the corresponding hydrogen bond acceptor positions in the binding site. Projected points allowed the possibility of structurally dissimilar active compounds forming hydrogen bonds to the same location, regardless of their point of origin and directionality. Each pharmacophoric site was first assigned with an energetic value equal to the sum of the Glide XP contributions of the atoms comprising the site. This allowed the sites to be quantified and ranked on the basis of these energetic terms. Glide XP descriptors include terms for hydrophobic enclosure, hydrophobically packed correlated hydrogen bonds, electrostatic rewards, π-π stack pairing, π ⋯ cation pairing, and other interactions. Sites where less than half of the heavy atoms contributed to the pharmacophore feature were excluded from the final hypothesis. Thus, if only two heavy atoms in a six-member ring exhibited energetic interactions, the ring was not considered a pharmacophore feature.

E-pharmacophore-based database screening and binding free energy calculations
For the e-pharmacophore-based screening approach, explicit matching was required for the most energetically favorable sites which have scored better than −1.0 kcalmol −1 . Multiple sites were included in case of more than one site had the top score. Screening molecules were required to match a minimum of three sites for a hypothesis with three or four sites and a minimum of four sites for a hypothesis with five or more sites. The distance matching tolerance was set to 2.0 Å as a balance between stringent and loose-fitting alignment. The generated pharmacophore hypothesis was screened against the compounds available in the Sweet Lead database (https://simtk.org/home/sweetlead). Database hits were ranked on the basis of fitness score, a measurement of how well the aligned ligand conformer matches with the hypothesis based on RMSD site matching, vector alignments, and volume terms. The fitness scoring function is an equally weighted composite of these three terms and ranges from 0 to 3, as implemented in the default database screening of Phase module. The ligands with the best fitness scores were docked into the binding sites of the VDR and CYP24A1. The free energy binding of selected drug-like molecules was scored on the basis of Prime/MM-GBSA method (Lyne, Lamb, & Saeh, 2006;Prime, 2015). The final best docking pose was selected for the calculation of binding free energy by Prime/MM-GBSA method. The docked complexes were minimized with local optimization feature of Prime. The OPLS-AA and GB/SA continuum solvent model was applied to calculate the binding energy. The Qikprop module was used to calculate the absorption, distribution, metabolism, and excretion (ADME) properties of the selected potent hits (QikProp, 2015).

Guner-Henry (GH) scoring method
The generated pharmacophore model was validated by applying GH scoring method. This method is used to calculate the accuracy of the hits and the recall of actives from a data-set which consists of known actives and known inactives. The quantification of model selectivity and coverage of activity space of database mining were successively quantified by the GH scoring methodology. It can also be used to evaluate the effectiveness of similarity search in a database containing both structural and biological activity data (Guner, 2000). The GH score ranges from 0 to 1. The 0 indicates the null model, whereas 1 indicates the ideal model (i.e. containing all of, and only the active ligands). The normal GH value is expected to be more than .7 (Yang & Shen, 2005). %Y, recall accounts for the percent yield of actives in a database and %A, precision accounts for the percent ratio of actives in the hit list. The GH scoring method can be applied to find best tolerance of RMS in fixed atomic positions. An optimum GH score was calculated for fixed position tolerance. The GH scoring method was applied to the potent drug molecule with their 600 conformers and for the decoy data-set molecules in order to validate the pharmacophore models. The following variables were computed for the GH analysis: (a) the percent active (%A) (the coverage of activity space of the model), (b) the percent yield (%Y) (A measure of selectivity of the model), (c) the enrichment factor (E), and (d) the G score. All these variables were calculated by the following formula, The total number of compounds in the drug database (D), the number of actives in the database (A), the number of actives retrieved by the model (H a ), the total number of hits retrieved from the model (H a ), and the total number of hits retrieved by the model (H t ).

Reverse validation
To validate the generated pharmacophore models and active sites of both the targets, we has reverted the complete strategy. Similar parameters were adopted like the actual process for molecular docking and pharmacophoric mapping.

MD simulation studies
The best screened protein ligand complexes were subjected to the molecular simulation studies. The GROMACS 4.5.5 package (Hess, Bekker, Berendsen, & Fraaije, 1997) was used for all the MD simulation studies with the GROMOS96 (53a6) force field (Oostenbrink, Villa, Mark, & Van Gunsteren, 2004). The topology file was created on the basis of PRODRG server (van Aalten et al., 1996). The dimension of the cubic box is 10.60 × 11.21 × 10.51 (all in Å) in case of VDR models and 11.24 × 12.54 × 10.58 (all in Å) in case of CYP24A1 models. The protein was solvated in the simple point charge (SPC) water model in a cubic cell. Proper counter ions were added to neutralize the system. 50,000 minimization steps were carried out with steepest descent method with the tolerance of 1000 kJ mol −1 nm −1 . Particle Mesh Ewald method (Kawata & Nagashima, 2001) was used to apply electrostatic and van der Waals (vdW) interactions with a cutoff distance of 1.0 nm for short-range neighbor list (rlist) and 1.0 nm for coulomb cut-off (rcoulomb) and the same (1.0) were fixed for vdW interactions (rvdW). The geometry and solvent orientations of the final minimized and SPC solvated system measured the realistic structure. LINCS and SETTLE algorithms (Thangapandian, John, Sakkiah, & Lee, 2010) were used to constrain the bond angles and geometry of the water molecules. System pressure was maintained with 300 K, 1 bar with the help of V-rescale. Further, 30 ns of molecular simulation had run with a time step of 2 fs. The atomic structural coordinates were saved in each 2 ps.

Binding free energy calculation studies
Binding free energies of all the compounds to VDR and CYP24A1 were calculated using MM-PBSA program in Gromacs described by Kumari and Kumar (2014). For each complex, 100 snapshots were taken from the last 10 ns of the MD trajectory with an interval of 100 ps. Before performing the calculations, periodicity was removed in each complex.

Dual-selectivity design strategy
Two different pharmacophore modeling approaches that are widely used in computer-aided drug discovery research (Nagamani, Kesavan, & Muthusamy, 2012).
The first approach is based on a ligand with experimentally known inhibition for a particular therapeutic target (Pharmacophore modeling), whereas the second one exploits the receptor-based information for designing E-pharmacophore models. Our strategy was based on the information extracted from the receptor-based approaches to identify the dual functioned molecules. E-pharmacophore modeling could be a valid one since it merges the strength of ligand-based (pharmacophore modeling) and receptor-based (molecular docking) approaches together. At first, the dual activity compound (CTA-018) was prepared and docked within the active site of the two target proteins (VDR and CYP24A1) to predict their binding affinity and molecular interactions. The second step was E-pharmacophore generation of both the target proteins using two different docked complexes (VDR-CTA018 and CYPA24A1-CTA018). In the third step, the drug-like database (Sweet Lead database) was collected and docked with VDR. The compounds which bound very well in the active site of VDR have chosen for the next step. The docking hits from the VDR protein were mapped over the derived pharmacophore models. The compounds which accommodated very well upon both the pharmacophoric sites were docked in the active site of the second protein (CYP24A1). Finally, retrieved compounds from the second set of docking were known as possible potent dual selectivity molecules. The flowchart explaining the dual strategy is shown in Figure 1.

Molecular docking of CTA-018
Molecular docking was carried out using Glide module of Schrodinger software. In this study, known dual activity compound (CTA-018) was docked separately in the catalytic domain of two different target structures, namely VDR (PDB id: 3VTC) and CYP24A1 (PDB id: 3K9 V). In VDR CTA-018 complex, we could find two different interactions. The oxygen atom of the hydroxyl group of Ser274 was well interacted with the hydrogen atom of the CTA-018 molecule (HO⋯OH, 1.887 Å). Another interaction was observed in His393. The hydrogen atom from the amino group (-NH) of His393 well interacts with the oxygen atom of (O=S=O) CTA-018 molecule (NH⋯O=S=O, 1.875 Å). In CYP24A1 CTA-018 complex, one interaction has been observed in the allosteric site. The oxygen atom of the carboxyl group (C=O) of Met246 interacted with an oxygen atom of -OH group in CTA-018 molecule.
E-pharmacophore generation and validation E-pharmacophore generation is the combined aspects of ligand-based and structure-based approaches to enhance enrichments rather than ligand information alone. The method presented here attempts to take a step beyond the simple contact scoring since it incorporates both the structural and energetic information using the scoring function of Glide XP. Six pharmacophoric features were predicted in the VDR and CYP24A1, but only five pharmacophoric sites in VDR and two pharmacophoric sites in CYP24A1 were chosen on the basis of site score. The sum of Glide XP docking of the atoms present in the active site of the protein was equal to the energetic

GH studies
The best hypothesis was validated by GH method. We computed the following variables: (a) %A is the percent-age of known active compounds retrieved from the database (precision), (b) H a is the number of actives in the hit (true positives), (c) A is the number of active compounds in the database, (d) %Y, the percentage of known actives in the hit list (recall), (e) H t is the number of hits retrieved, (f) D is the total number of compounds in the database, and (g) GH is the Guner-Hery score. The best models (ADHHH and DD) succeeded in retrieving 76% of the compounds, and an enrichment factor of 2.24 and a GH score of .75, which shown the good quality of the model generated after the analysis (Table 2).

Molecular docking of drug-like database for VDR binding
Sweet Lead is a database containing 9127 compounds from the world's approved medicines and isolates from traditional medicinal herbs. This database provides us the accurate structures for drug discovery efforts and chemo informatics analysis. Initially, Sweet Lead database compounds were screened against the VDR active site by molecular docking. Molecular docking is a computational approach which samples conformations of small compounds at protein-binding sites. The best  Notes: H t is number of hits retrieved, H a is number of actives in hit list, D is number of compounds in a database, %A is a ratio of actives retrieved in hit list, %Y is fractions of hits relative to size database (hit rate of selectivity), E is enrichment of active bin by model relative to random screening is GH is Guner-Henry score.
complement of the protein-binding sites was assessed by scoring function. The quality of the docking method is assessed by two main aspects: (i) docking accuracy and (ii) screening enrichment. In this study, the data-set compounds were screened against the VDR active site by Glide High Throughput Virtual Screening (HTVS). A total of 7412 compounds were screened through this process, further it was screened by standard precision (SP) and extra precision (XP) mode of docking. Finally 854 compounds were identified as possible potential hit molecules for VDR.

Pharmacophore-based screening
As mentioned in the previous section, two different pharmacophore models (ADHHH for VDR and DD for CYP24A1) were generated and screened against the Sweet Lead database. Initially, 854 screened compounds underwent for the first set of pharmacophore screening with ADHHH hypothesis. We retrieved the compounds which had a fitness score of more than 1.5. There were 24 compounds which fitted with this hypothesis. These 24 compounds were screened against the DD hypothesis with the same criteria. Finally, 12 compounds were survived with both the hypothesis. These 12 compounds were docked in the active site of CYP24A1 protein.

Reverse validation
The Sweet Lead database was first docked in the active site of CYP24A1 enzyme and 642 compounds were selected which had good glide score when compared to CTA-018 molecule. These compounds were mapped over both VDR and CYP24A1 pharmacophore models. Among these 48, compounds were selected based on the best fitting over the VDR and CYP24A1 pharmacophore models. Finally, these 48 compounds were docked in the active site of VDR. This reverse validation procedure also retrieved the previously identified four compounds against the VDR and CYP24A1protein.

MD simulations studies of the drug molecules with VDR and CYP24A1
The behavior of the protein-ligand complexes can be studied in detail through MD simulation studies. Native contacts in biological interactions are the major determinant to analyze the dynamic properties of the protein  motions. In this context, computational tools are more valuable for assessing protein flexibility and to reveal the dynamic pattern in complex dynamics (Bahar, Lezon, Yang, & Eyal, 2010;Leo-Macias, Lopez-Romero, Lupyan, Zerbino, & Ortiz, 2005). After initial analyses, we found that the trajectories were significantly stable and the produced simulation was reliable for further studies.

Structure flexibility
The behavior of the newly identified lead compounds in the active site of VDR and CYP24A1 were studied in detail with the aid of MD simulation studies. The root mean square deviation (RMSD) of backbone atoms was calculated with respect to the starting structure as a function of time. The RMSD plots of the VDR models (VDR CTA-018 , VDRSW-03611 , VDR SW-03852 , VDR SW-03716 , VDR SW-03939 ) are displayed in Figure 6(A) and the RMSD plots of the CYP24A1 models (CYP24A1 CTA-018 , CYP24A1 SW-03611 , CYP24A1 SW-03852 , CYP24A1 SW-03716 , CYP24A1 SW-03939 ) are displayed in Figure 6(B). Overall, the RMSD was not fluctuated significantly, showing an average RMSD of less than .5 nm for most of the simulation period. The first 1 ns of the initial run was critical for the equilibration and later all the protein ligand complexes attained stable protein-ligand complexes. In both the cases, the identified lead compounds displayed different types of fluctuation when compared to the wild type (VDR CTA-018 and CYP24A1 CTA-018 ). In VDR models, the RMSD curve shows stability (except compound SW-03716) after 1 ns of simulation with .03 nm deviations and it also shows that the average RMSD of the complex is .27 nm with a SD of .03 nm. In CYP24A1 models, the RMSD curve shows stability near 20-30 ns time period with .04 nm of deviation. The average RMSD of the CYP24A1 models were .30 nm with a SD of .05 nm in MD studies. The radius of gyration (Rg) of the peptide backbone as a function of time could also be used as a probe for structural flexibility. In VDR models, all the protein structures maintained a stable structure which was around 1.8 nm, but the compound SW-03716 and compound SW-03939 was maintained around 1.90 nm. After 25 ns, both the compounds (SW-03716 and SW-03939) maintained stable conformers around 1.90 nm. In case of CYP24A1 models, all the protein structures sustained a stable structure around 2.15 nm. But the compound SW-03611 was deviated more (2.30 nm), since the compound is less tightly packed when compared to the remaining compounds. However, the compound has bound well within the active site of CYP24A1. Both the results substantiated the protein backbone stability in the presence of drug molecules. The Rg's of the protein backbone is displayed in Figure 6(C) and (D). The overall structural properties indicate that the drug-like molecules are stable in both the protein active sites.

Ligand flexibility
The conformational changes and ligand flexibility in the active site of VDR and CYP24A1 proteins were also assessed through MD simulations. The RMSD plots of the ligand molecules are displayed in Figure 6(E) and (F), respectively. In VDR models, there is no significant change in the binding pattern of all the drug molecules.
All the ligands were fluctuated vigorously up to 5000 ps in VDR active site. After 3000 ps, all the drug molecules produced a stable conformation with the average of being .10-.25 nm. Among these compound SW-03852 produced a good s conformer when compared to the other drug molecules. In case of CYP24A1 models, all the drug molecules generated a stable conformation of (.21 nm) when compared to the wild one (.20 nm). Among these compounds, SW-03611, SW-03852, and SW-03939 produced stable conformer after 5 ns, but the compound SW-03716 fluctuated vigorously throughout the MD simulations (Figure 6(E) and (F)). The overall results suggest that there is not much deviation in the active site of VDR and CYP24A1 models and the RMSD has been constantly maintained during the whole simulation period.
In order to understand the nature of the binding of drug molecules in the active site of VDR and CYP24A1, the time dependence of hydrogen bonds between receptor and drug-like molecule was monitored during the whole simulation period. The results showed that the calculated intermolecular hydrogen bonds were well correlated with the biological activities of the lead compounds. All the drug molecules exhibited constant hydrogen bonding interactions throughout the simulation time period, which gives direct clues about its affinity towards the protein-ligand complexes (Figure 7(A) and (B)).
The solvent-accessible surface area (SASA) which shows the free energy of solvation (ΔG slov ) was also calculated for the identified drug molecules along with VDR and CYP24A1 enzyme (Figure 7(C) and (D)). The results show that binding of the lead compounds does not induce any significant change in SASA. In both VDR and CYP24A1 models, all the drug molecules including the best compound (CTA-018) exhibits a similar type of distribution throughout the simulation time. The VDR models show an average solvation energy of 550.65 nm with a SD of 1.6 nm, whereas CYP24A1 models show an average solvation energy of 128.35 nm with a SD of 15.23 nm.

Ligand binding site
The agonist and antagonist mechanism of VDR analogs was first defined by Kakuda et al. (2010) using X-ray crystallographic studies. Expression of VDR mutations (His301Phe, His393Phe) clearly exhibits the role of these amino acids in agonistic and antagonistic activity. Interestingly, the identified drug-like molecules were mainly interacted with the His393 amino acid (Figure 7(E)). In 2010, Posner et al., suggested that the binding pocket of CTA-018 is different from the azole-based compounds (e.g. ketoconazole) in CYP24A1 enzyme. To understand the dynamic behavior of His301 and His393 in VDR and binding mode of CYP24A1 enzyme, we analyzed the distance between the compounds and N, O, and H atoms of the key residues at binding site throughout the 30,000 ps simulations.
In VDR, the distance of His393 and the identified drug-like molecules was measured throughout the simulation time, since it plays an important role in the agonistic mechanism. Initially, in the VDR-CTA018 complex, a drastic fluctuation was observed and after 5 ns the distance was maintained at reasonably constant at .30 nm. The total distance of VDR and all the drug molecules shows an average deviation of .30 nm (SD of .03 nm). It is very close to the active site amino acids. All the drug molecules produced a good stable conformers (Figure 7(E)). This interaction led to correct helix 12 folding to induce VDR mediated expression. The presence of hydroxyl atom (-OH) in all the identified drug molecules played major role in the agonistic mechanism of the VDR protein. In CYP24A1 enzyme, we mainly calculated the distance between the HEM prosthetic group and the drug molecules. Posner et al. (2010) reported the CTA-018 molecule was bound in a different binding pocket unlike azole-based compounds. Thus, we calculated the distance between the HEM prosthetic group and the drug molecules. The total distance between HEM and drug molecules of the CYP24A1 shows an average deviation of 1.00 nm (SD of .03 nm). Our previous work suggested that CTA018 bound in different binding site when compared to ketoconazole (potent inhibitor of CYP24A1). The calculated distance between HEM and ketoconazole was 1.4 nm, whereas the calculated distance between HEM and CTA-018 is .9 nm. Interestingly, in the present work the calculated distance between HEM and the identified drug molecules was around 1.00 nm (Figure 7(F)). We speculate that the identified drug molecules are able to generate different binding pose when compared to ketoconazole.

Analysis of binding free energy
The strength of interaction between a receptor and its ligands can be calculated by MM-GBSA and MM-PBSA calculations using MD trajectories (Srivastava & Sastry, 2012). Kuhn, Gerber, Schulz-Gasch, and Stahl (2005) explained the application of MM-PBSA energy function using MD snapshots. Moreover, the binding affinities of the saquinavir (protease inhibitor) to the wild type and HIV-1 protease resistant variant were calculated using MM-PBSA technique by Coveney and co-workers (Stoica, Sadiq, & Coveney, 2008). Chen and his co-workers compared the QM/MM MD simulation and QM/MM GBSA calculation to study the accuracy of the binding affinity of serine proteinases (Chen, Wang, Zhang, Chen, & Zhu, 2015). Shan and his co-workers compared the binding free energies of C-domain-Ang II and N-domain-Ang II complexes using MM -PBSA method (Guan, Han, Zhang, Wang, & Shan, 2015). Other than hydrogen bond interactions, the contribution of van der Waals interactions and electrostatic interactions also can be studied in detail using PBSA method (Arba & Tjahjono, 2015;Wang, Yang, Shi, & Le, 2015). Here, we used MM-PBSA method to calculate the binding free energy of each set of protein-ligand complex in order to compare the binding affinity between the CTA-018 with the identified potent lead compounds. The binding free energies were calculated for the last 100 snapshots extracted from the last 10 ns of MD trajectory at an interval of 100 ps, since taking up very limited number of states from MD simulation is the best practical method. Table 1 displays the binding free energies of the identified lead compounds in the active site of VDR and CYP24A1 along with the best dual function molecule CTA-018. The corresponding binding free energies have ranged from -60.42 to -48.62 kcal/mol for VDR models and in CYP24A1 models it ranges from -90.54 to -70.42 kcal/mol. The results have revealed that the drug molecules bind strongly with VDR and CYP24A1 in accordance with the experimentally proved CTA-018. The overall results indicate that the current energy analysis performed by MM-PBSA method is reliable.
The van der Waals interactions (ΔE vdw ) favor the binding of the lead molecule in the active site of VDR and CYP24A1. For all the complexes the van der Waals energies are not very different. Thus, they are all having good hydrophobic contacts. Moreover, the entropy charges mainly depended on inhibitors associations and weaken the bindings. According to Table 4, the polar interactions of all the drug molecules with VDR and CYP24A1, including the electrostatic interactions with the polar interactions (ΔG ele + pol) vary from each other in the same set of protein. Overall, the results from the binding free energy calculations clearly exhibit the identified drug molecules are able to produce the similar-type interactions and energy when compared to the already reported compound CTA-018.

Amino acid decomposition analysis (ADA)
The contribution of the each amino acid in total binding free energy was computed by the evolution of Lennard-Jones (LJ) and coulombic interaction energies. The results of the ADA are shown in Figures 8 and 9, and Tables 5 and 6. Residues made a favorable contribution to ligand-receptor interactions via negative energy shifts.
The 2D interaction diagram generated from the Schrodinger software was used to detect the interacting residues in each case. The identified drug molecules shared similar binding pattern-like CTA-018 in both the cases (VDR & CYP24A1). Amino acids His301 and His393 in VDR protein and Arg128 and Glu329 amino acids in CYP24A1 protein play an important role in ligand binding activity. Residues constructing hydrophobic pocket in the proximity of His301 and His393 in VDR protein are almost involved in all types of interactions with ligand.
Initially, in the first analysis, we discussed the VDR agonistic nature of the identified drug molecules. In case of VDR CTA-018 scaffold, His393 was found to be the most significant residue in ligand-receptor interactions (ΔE = −8.28 kcal/mol). The hydrogen atom from the amino group of His393 interacted with an oxygen atom of the sulfonyl group (O=S=O). The hydrophobic contacts between these groups made it a favorable interaction. In addition to it, one more interaction was available with Ser274 (ΔE = −2.35 kcal/mol). The hydrogen atom from the hydroxyl group of the CTA-018 interacted with the oxygen atom of the Ser274. But overall His393 produced more energy when compared to Ser274.
In VDR SW-03611 scaffold, we found a single interaction with His393. The hydrogen atom from the hydroxyl group of His393 interacted with the nitrogen atom in the drug molecule. However, this interaction produced more binding energy when compared to the remaining drug molecules (ΔE = −7.61 kcal/mol). This His393 had maximum coulombic and LJ interaction energies when compared to the other identified drug molecules (-3.50 kcal/mol, -2.38 kal/mol, respectively).
In VDR SW-03716 , we found three different interactions in the VDR active site. But the His393 yields more energy (ΔE = −6.31 kcal/mol) when compared to the remaining two interactions (ΔE = −5.23 and −2.12 kcal/mol, respectively). The oxygen atom of the carboxyl group of Ser274 interacted with hydrogen atom of the hydroxyl group in the drug molecule. Another interaction was between hydrogen atom from the amino group of Tyr143 and the oxygen atom of the hydroxyl group.
In VDR SW-03852 binding scaffold, unfortunately the main interaction was available with His301. But the binding energy was comparably low when compared to all other drug molecules (ΔE = −3.32 kcal/mol). The hydrogen atom from the hydroxyl group interacted with the nitrogen atom of -CN group of His301 amino acid. Another interaction was available between the oxygen atom of the carboxy group of Ala299 and the hydrogen atom of hydroxyl group in the drug molecule. Moreover, the SW03852 had -OH group in its structure. Thus, it may interact with His393.
In VDR SW-03939 binding scaffold, the main interaction was observed in His393. The interaction energy of His393 was −7.29 kcal/mol. The interaction with Tyr143 had second better interaction (ΔE = −5.21 kcal/mol) in all VDR complexes. The oxygen atom of the carboxyl group (C=O) of Ser233 interacted with hydroxyl group of drug molecule (ΔE = −4.28 kcal/mol). The same type of interaction had also been observed in the Ser274 (ΔE = −4.32 kcal/mol). The hydrogen atom from the hydroxyl group interacted with the oxygen atom of the carboxyl group of Ser271 (ΔE = −6.25 kcal/mol).
In the second set, we analyzed the CYP24A1 antagonistic nature of the identified drug molecules in the active site of CYP24A1. Out of four identified drug molecules, three drug molecules had main interaction with Arg128. Thus, the Arg128 may be the most  significant residue in ligand receptor interactions. Interestingly, the Arg128 is placed in the catalytic core of the enzyme near to the HEM prosthetic group. However, the identified drug molecules may target the HEM group via Arg128 interaction unlike azole-based compounds.
In CYP24A1 CTA-018 complex, we found main interaction with Leu325 and Glu329. Both the amino acids had interaction in the same pattern. The oxygen atom of the carboxyl group of Leu325 (ΔE = −5.45 kcal/mol) and Glu329 (ΔE = −6.85 kcal/mol) interacted with the hydrogen atom of the CTA-018.
In CYP24A1 SW-03611 complex, we found three different interactions. The oxygen atom of the hydroxyl group of Leu129 and Met148 interacted with the oxygen atom of two different carboxyl group (C=O). The binding energies were −2.12 and −3.16 kcal/mol, respectively. One more interaction was also observed between the oxygen atom of Arg128 and hydrogen atom of amino group (−NH) of the drug molecule with the binding energy of −4.15 kcal/mol.
In CYP24A1 SW-03716 complex, we found one hydrogen bond interaction with the binding energy of −5.15 kcal/mol. The hydrogen atom from the hydroxyl group of Thr394 interacted with the oxygen atom of the hydroxyl group in the drug molecules.
In CYP24A1 SW-03852 binding scaffold, we found two different interactions. The oxygen atom of the carboxyl group of Arg128 interacted with the amino group (−NH) of drug molecule (ΔE = −5.45 kcal/mol). Another interaction was observed between the hydrogen atom of Leu129 and the oxygen atom in carboxyl group of the drug molecule (ΔE = −3.16 kcal/mol).
In CYP24A1 SW-03939 binding scaffold, three different interactions were observed in the catalytic site of the CYP24A1 enzyme. The oxygen atom of the hydroxyl group of Leu79 and His497 interacted with the oxygen atom of the carboxyl group of the drug molecule with the binding energy of −3.61 and −2.61 kcal/mol, respectively. Another π-π stack pairing interaction was observed between the drug molecule and Arg128 amino acid.
Arg128 played an important role in the context of CYP24A1 binding mode. The binding energy was reasonably high in Arg128 in all the complexes when compared to other hydrogen bonding interactions. Even though, the active sites of the newly identified drug molecules were entirely different, the binding pattern was similar to CTA-018 as we mentioned in the MD simulation studies. Thus, the identified drug molecules may have better inhibitory activity against CYP24A1 enzyme.
Overall, our in silico approaches supported the identification of drug molecules may act more potent than CTA-018 molecule. The ADME properties of the druglike molecules are also under acceptable range (Table 7). Moreover, the identified compounds are from Sweet Lead database which contains approved medicines and isolated from traditional medicinal plants. So, the identified drug molecules may have the high potency to treat osteoporosis among CKD patients without any side effects. Further, in vitro or in vivo analysis may lead to the identification of potent drug candidate among the selected compounds to manage CKD patients.

Conclusion
In summary, an approach combining molecular docking and E-pharmacophore filtering was applied to meet the critical challenges faced in designing efficient multi-target compounds to treat complex sHPT in CKD. Initially, Sweet Lead database was docked in the active site into one of the two targets (VDR) using docking program Glide, and the compounds with better binding characteristics were selected for further study. Common E-pharmacophore models were developed for both the targets with the experimentally known compound (CTA-018). Post processing of the docking results with pharmacophore filtering allowed us to bypass the difficulty with scoring functions that produced false positives from the molecular docking calculations. The selection of compounds to the required chemical features was confirmed by the pharmacophore filtering method. The compounds fitted on both the pharmacophoric models were docked into the active site of the second target CYP24A1. Finally, the compounds with better binding characteristics and with the required pharmacophoric features were selected as possible dual function compounds against VDR and CYP24A1 proteins. From the Sweet Lead database, we found four promising compounds which have the possible dual parallel function against the targets VDR and CYP24A1. These could be promising compounds for the management of CKD.

Supplementary material
The supplementary material for this paper is available online at http://dx.doi.org/10.1080/07391102.2015. 1111168. Predicted octanol/water partition coefficient -2.0-6.5. d Predicted aqueous solubility, log S. S in mol dm −3 is the concentration of the solute in a saturated solution that is in equilibrium with the crystalline solid -6.5-.5. e Predicted IC 50 value for blockage of HERG K+ channels. Concern below -5. f Predicted human oral absorption on 0 to 100% scale. >80% is high <25% is poor. g Lipinski rule of fivemaximum is 3.