D 3 and its analogs with the human VDR receptor model

1,25-dihydroxyvitamin D(3) has quite significant anticancer properties, but its strong calcemic effect in principle excludes it as a potential anticancer drug. Currently, a lot of effort is being devoted to develop potent anticancer analogs of 1,25-dihydroxyvitamin D(3) that would not induce hypercalcemia during therapy. In this work, the free binding energy of the VDR receptor with 1,25-dihydroxyvitamin D(3) and its three potent analogs (EB 1089, KH 1060 and RO 25-9022) is calculated and compared with each other. With this approach, we could estimate the relative binding affinity of the most potent analog, RO 25-9022, and also revealed a quite distinct mechanism of its interaction with VDR.

Vitamins D are a group of lipid-soluble organic compounds causing a variety of physiological responses.The most important vitamins of the vitamin D group are vitamin D 2 (ergocalciferol) produced by plankton and fungi from ergosterol and vitamin D 3 (cholecalciferol) produced by animals and humans from 7-dehydrocholesterol.Both compounds are synthesized when an organism is exposed to UV radiation (Holick, 2008).Vitamin D 3 is biologically inactive; it undergoes hydroxylation at the C25 position and conversion to prohormone calcidiol (25-hydroxycholecalciferol) in the liver, and, subsequently, hydroxylation at the 1αC position and conversion in kidneys to calcitriol (1,25-dihydroxycholecalciferol or 1,25(OH) 2 D 3 ) (Guyton et al., 2003), a biologically active form of vitamin D 3 .The extent of biological activities of vitamin D is remarkable.The biological functions of vitamin D we are currently aware of include regulation of the intestinal absorption of calcium and phosphates (Meyer et al., 2006;Barthel et al., 2007), autoimmune disease prevention (Cantorna et al., 1996;Zella & DeLuca, 2003;Cantorna & Mahon, 2004), osteoclastogenesis and mineralization of bones (Fretz et al., 2006;Milat & Ng, 2009;Haussler et al., 2010), antiproliferative activity and prevention of prostate cancer (Zhuang & Burnstein, 1998;Boyle et al., 2001;Moffatt et al., 2001), breast cancer (James et al., 1996;Narvaez et al., 1997;Wu et al., 1997;Malinen et al., 2008), colon cancer (Evans et al., 1999;Gaschott et al., 2001), pancreatic cancer (Kawa et al., 1997;Schwartz et al., 2008;Chiang et al., 2009) and ovarian cancer (Li et al., 2004), mainly through the regulation of the expression of cyclin-de-pendent kinase inhibitors (p21 and p27), and it seems to be a good candidate for anticancer therapies (Nickeleit et al., 2007).1,25(OH) 2 D 3 bioactivity involves mainly binding with the vitamin D receptor (VDR) (Rochel et al., 2000).The VDR molecule has two binding sites for the 1,25(OH) 2 D 3 molecule: the main (genomic) binding site responsible for the regulation of transcription of various genes and indirectly for the biological effects mentioned above, and an additional (alternative) binding site responsible for the "rapid response", because certain specific biological effects can appear after 1-2 minutes.1,25(OH) 2 D 3 binding to the genomic site occurs in the VDR molecule localized in the nucleus, whereas binding to the alternative pocket usually occurs in the caveolae (Haussler et al., 2011).
Such a significant role of 1,25(OH) 2 D 3 in the human body drew attention of many research groups which initiated projects toward the ultimate goal: design and synthesis of new, potent and selective 1,25-dihydroxyvitamin D 3 analogs (Brown & Slatopolsky, 2008).It is especially promising as novel anticancer compounds are at hand, although the mechanism by which various 1,25(OH) 2 D 3 analogs elicit specific responses is yet unknown (Singarapu et al., 2011).
In this paper, the binding energies of 1,25(OH) 2 D 3 and selective and potent VDR agonists with significant anticancer properties (Leo Pharma EB 1089, Leo Pharma KH 1060 and Hoffmann-La Roche RO 25-9022, being one of the most bioactive 1,25(OH) 2 D 3 analogs in terms of the stimulation of cellular differentiation and inhibition of cellular proliferation) (Uskokovic et al., 2001;van den Bemd & Chang, 2002;Guyton et al., 2003;Fig. 1) to the genomic binding site of reduced VDR (redVDR) were assessed and compared with one another.

MatERIals anD MEtHODs
Crystallographic structures of the published VDR/ligand complexes are characterized by the almost identical conformation of the ligand binding domain (genomic site) (Rochel et al., 2000;Tocchini-Valentini et al., 2001;2004;Singarapu et al., 2011).The structures of the reduced VDR receptor and of 1,25(OH) 2 D 3 were constructed on the basis of the crystallographic structure of 1,25(OH) 2 D 3 /VDR (PDB entry: 1DB1).The remaining ligand structures were built on the basis of the 1,25(OH) 2 D 3 structure with a conserved trans configura-tion of the C6-C7 bond, conformation of the cyclohexane ring (A-ring), octahydroindene ring (C,D-rings) and the angle between the ring planes of 30° (Rochel et al., 2000;Haussler et al., 2011).The redVDR model was constructed using the following amino acids from the 1DB1 structure: Tyr143-Lys240, Ala267-Lys322, Leu393-Gly423 (Fig. 2).The addition of hydrogen atoms to ligands and redVDR was carried out using the Chimera program (Pettersen et al., 2004).Ionizable amino acid residues were protonated according to physiological (cellular) pH.The structures of ligands and redVDR were parameterized with Gasteiger atomic partial charges using Chimera.The molecular docking procedure was performed using Auto Dock software (Morris et al., 1998).The docking procedure was as follows: the ligands were docked inside a cuboid grid box with a spacing of 0.375 Å and dimensions of 19.5 × 19.5 × 14.25 Å which allowed sampling actual ligand-accessible regions.A Lamarckian genetic algorithm was used for 2.5 × 10 6 energy evaluations for each of 100 docking trials.To obtain the docked conformations of 1,25(OH) 2 D 3 and analogs, clustering at 0.5 Å was performed.The docking procedure used was a semi-flexible method which means that redVDR was treated as a rigid molecule and 1,25(OH) 2 D 3 and analogs had flexible torsion angles.During the docking procedure, the torsion angles C4-C3-O-H, C10-C1-O-H, C16-C17-C20-C22, C17-C20-C22-C23, C20-C22-C23-C24, C22-C23-C24-C25, C23-C24-C25-O and C24-C25-O-H were subject to change.In the calcitriol analogs, the torsion angles allowed to change were the same as in 1,25-dihydroxyvitamin D 3 with the exclusion of angles containing double bonds or a conjugated bonds system.Additionally, in EB 1089 and KH 1060, the rotational freedom of ethyl groups was taken into account.
Table 1 contains binding energies (Auto Dock score), RMS values and hydrogen bond parameters for the docked structures of 1,25(OH) 2 D 3 and analogs.
In the second part of our calculations a rescoring protocol based on a method developed by Fanfrlík et al. (2010) was employed as the results of calculations using this method were found to be consistent with known experimental data.Moreover, we used this protocol earlier to propose a novel binding mode of epothilone A to β-tubulin (Kamel & Kolinski, 2011).
To generate a ligand conformation that corresponds to the "global" potential energy minimum, a molecular dynamics (MD) method was used for 1,25(OH) 2 D 3 and the analogs.MD simulations were run using the AMBER force-field (Cornell et al. 1995) at a constant temperature of 1000 K using the Andersen thermostat.Equilibration time was 5 ps and production run was 10 ps.One-hundred resulting structures were generated for each ligand.Subsequently, the ligand/redVDR complexes and the redVDR molecule were geometrically optimized in a COSMO continuous solvent model (Klamt & Schuurmann, 1993) and in the gas phase using a PM6 semiempirical quantum mechanics method (Stewart, 2007) with dispersion and hydrogen bond corrections (PM6-DH2, Korth et al., 2010) with the MOZYME module (Stewart, 2009) using localized molecular orbitals implemented in MOPAC2009 (Stewart, 2008).The resulting structures from docking and MD were geometrically optimized using the PM6-DH2 method in continuous solvent.Apart from this methodology, an alternative approach was also attempted, which was a direct geometry optimization of ligand structures built on the basis of 1,25(OH) 2 D 3 .Such an approach was based on the fact that the backbone of all the ligands was almost identical,  In this part of the work, our aim was to calculate the free binding energy of 1,25(OH) 2 D 3 and analogs with the redVDR molecule.We adopted an approach based on the thermodynamic cycle proposed by Raha and Merz (2005): (1) (2) ΔG w bind is the estimated free binding energy in the water environment, ΔH f w (X) is the heat of formation in aqueous solution, where X stands for the free protein, free ligand or protein-ligand complex.Likewise, ΔH f (X) are the heats of formation in the gas phase of the free ligand, free protein or protein-ligand complex.ΔH f complex (X) is the enthalpy of the ligand or the protein molecule in the complex conformation.ΔG int is the interaction energy in the gas phase.Heats of desolvation and deformation are: X stands for the ligand or protein molecule.All heats of formation were calculated at 298 K. Entropy changes (ΔS) were calculated in vacuo using the AMBER forcefield.

REsults anD DIsCussIOn
The molecular docking of the 1,25-dihydroxyvitamin D 3 molecule its and analogs to redVDR generated complexes very similar to corresponding crystallographic structures.For the 1,25(OH) 2 D 3 /redVDR complex all hydrogen bonds present in the crystallographic structure were reconstructed (see Fig. 3A, Table 1).The hydrogen bonding pattern observed in the remaining crystallographic structures was also reconstructed (Fig. 3B).The most important contacts defining bioactivity are O1-H hydrogen bonds with Ser237 and Arg274, and O25-H with His397 and His305.An exception to the rule is the EB 1089/redVDR complex in which the O25•••H-N-His397 hydrogen bond was broken, but the polar contact was maintained (3.62 Å).This effect could be attributed to the high rigidity of the conjugated bonds system between C22-C23 and C24-C24a  and the longer side chain compared to 1,25(OH) 2 D 3 .
For the RO 25-9022/redVDR complex, the reconstruction of hydrogen bonds formed by O25-H with His397 and His305 and O1-H with Ser237 and Arg274 was rather unsurprising, as the RO 25-9022 side chain was of equal length to the 1,25(OH) 2 D 3 side chain and had only one double bond at C23, which did not impair its flexibility in a significant way.The binding pocket of the vitamin D receptor consists of many hydrophobic residues capable of interacting with 1,25(OH) 2 D 3 and analogs through van der Waals forces.Since the conformation and position of the backbone fragment inside VDR are virtually identical for all the compounds, with a slight exception of RO 25-9022 which lacks the methylene group at C10, we focused mainly on the side chain/protein contacts.The side chain/protein aliphatic van der Waals contacts are listed in Table 2.
The geometry optimization procedure at the PM6-DH2 theory level resulted in an improvement and increased number of contacts between the ligand side chains and the receptor molecule, as expected (see Table 3).The van der Waals contacts with Leu227, Ala231, Val234, Val300, His305, Leu309, His397, Leu414, Val418 and Phe422 observed in experimental structures (Tocchini-Valentini et al., 2004) were reconstructed.
It can be clearly seen from Table 3 that the number of aliphatic chain/protein contacts in PM6-DH2optimized complexes increased considerably in comparison with the number of contacts in experimental and docked structures (Tocchini-Valentini et al., 2004; Table 2).The main cause of this effect is probably the fact that the VDR crystallographic structures are always characterized by almost identical receptor molecule conformations, which in turn forces the almost identical ligand conformations.This is most probably due to the strong lateral contacts between VDR molecules and their packing that force the adoption of one, specific, geometrically favored conformation (Singarapu et al., 2011).The geometrically optimized structures are not part of a crystal structure, so they are not subjected to strong lateral interactions generated by neighboring VDR molecules and therefore they retain their significant number of degrees of freedom and can interact with ligand molecules as in the induced-fit model.In this way, the interactions between redVDR and a ligand molecule resemble those found in real biological systems.
The O25-H hydrogen bonds are among the most important interactions responsible for the ligand bioactivity, as noted above (Vaisanen et al., 2002).In the case of 1,25(OH) 2 D 3 and KH 1060 molecules, the length and angle of hydrogen bonds between O25-H and the two histidine imidazoles are very similar to each other, as clearly seen from Table 4.The length and rigidity of the side chain in the EB 1089 molecule places the O25-H hydroxyl group in a less favorable position, which weakens the O25−H•••N−His305 hydrogen bond and severely disrupts the O25•••H−N−His397 hydrogen bond (3.425Å; 160.2°).
The analysis of the RO 25-9022/redVDR complex O25-H hydroxyl group hydrogen bonding yielded qualitatively different results from those obtained for the compounds discussed above.During both in vacuo and continuous water model geometry optimization procedures using the PM6-DH2 method and the MOZYME function, the O25 hydrogen was transferred to the His305 imidazole nitrogen (NΕ atom).The cause of this effect is the extraordinary electron withdrawing property of the trifluoromethyl group (Swain et al., 1983).
Two simplified systems, each consisting of two histidine side chain fragments (from His305 and His397, respectively) and one of RO 25-9022 and one of the 1,25(OH) 2 D 3 side chain fragments (C23-C27) were built to obtain more quantitative results for this proton transfer effect.The systems were optimized using the PM6-DH2 method in the gas phase, in the continuous water environment and using a wide range of dielectric constants to simulate a protein interior environment (Schutz & Warshel, 2001).The analysis of geometry-optimized systems with the RO 25-9022 side chain fragment showed that proton transfer occurred for dielectric constant ≥ 4. Our preliminary calcula- tions using the Density Functional Theory with a small basis set (3-21G), a hybrid three-parameter functional (B3LYP) and COSMO solvation model carried out in the GAMESS program (Schmidt et al., 1993) showed similar results for dielectric constants ≥ 8.These findings correlate well with the dielectric constants applicable to protein interiors (Schutz & Warshel, 2001).It is obvious that such an effect was not observed in the system containing the 1,25(OH) 2 D 3 side chain fragment as the methyl groups at C25 have rather electron donating properties.The analysis of the free binding energy (or rather the "score" as a direct comparison of the results obtained from isothermal titration calorimetry and this method is not possible due to the differences in construction of the scoring function (Fanfrlík et al., 2010)) of 1,25(OH) 2 D 3 and analogs to redVDR yield-ed expected results consistent with experimental data.Our calculations clearly show that despite numerous side chain/protein contacts in the EB 1089/redVDR complex, its score has the least negative value that indicates the smallest binding affinity of this ligand (see Fig. 4), which is not very surprising given the weak hydrogen bond interactions of O25-H.1,25(OH) 2 D 3 has the highest binding affinity (score) of all of the compounds investigated, which is due to the very flexible side chain that facilitates a significant number of side chain/protein interactions and permits the formation of strong hydrogen bonds at O25-H.The high score achieved by the KH 1060 molecule could clearly be attributed to the length and flexibility of its side chain that enabled most side chain/protein contacts of all the molecules under investigation and provided a conformation necessary for strong hydrogen Panel A shows scores from the ligand geometry optimization procedure.Panel B shows scores from the molecular dynamics/ligand geometry optimization procedure.ΔH int -enthalpy of interaction in vacuo, -TΔS -entropy term, ΔG C -sum of desolvation and deformation enthalpies, -RTln(exp) -binding free energy with respect to the relative experimental binding constants.Free binding energy of 1,25-dihydroxyvitamin D 3 and its analogs bonding between O25-H and His397 and His305, respectively.The higher score obtained by RO 25-9022 could in turn be attributed to the proton transfer effect which certainly strengthens the hydrogen bonds formed and accepted by O25-H through electrostatic interactions.
A comparison of the results of our binding free energy calculations with the data provided by Tocchini-Valentini et al. (2004) showed very good correlation.In Tocchini-Valentini et al. (2004), relative values of ligand/VDR binding affinity were given.The binding of 1,25(OH) 2 D 3 was assumed as 100, whereas the relative binding affinities of KH 1060 and EB 1089 were determined as 90 and 80, respectively.Assigning the scores obtained from the simple ligand geometry optimization procedure (Fig. 4A) to the respective values given by Tocchini-Valentini et al. (2004), using the -RTln(experimental % of binding) equation produced relative binding affinities.The same approach was applied for the scores obtained from the MD/ligand geometry optimization procedure.On the basis of these results, we can estimate the relative binding affinity of RO 25-9022 to be about 92-93.

COnClusIOns
The Auto Dock docking/PM6-DH2 rescoring procedure employed in this study generated results in good agreement with those obtained by Tocchini-Valentini et al. (2004), which clearly shows the usefulness of this methodology in assessing the free binding energy of various compounds.
The direct application of a semiempirical quantum mechanics method (PM6-DH2) allowed for more precise calculations of the heat of interactions than would be possible using the molecular mechanics methodology only and we could study subtle effects such as the proton transfer effect observed.
The RO 25-9022 molecule has structural modifications that prevent it from being metabolized, therefore increasing its concentration in target cells and thus improving its bioactivity (Uskokovic et al., 2001).The substitution of the C26 and C27 methyl with trifluoromethyl groups is one of the most significant side chain modifications that prevents hydroxylation at these positions and in this way significantly alters the susceptibility to metabolism, but as it turned out, it also facilitates the proton transfer effect.
The calculation of the free binding energy (score) of the RO 25-9022 molecule with redVDR allowed us to estimate its high affinity and attribute it to the proton transfer effect and thus to strong hydrogen bonding with His397 and His305.Such an interaction mechanism results not only in a high binding energy, but it can lead to qualitatively distinct VDR conformation and, together with the high metabolic stability of RO 25-9022 molecule, could be responsible for its high biological activity.

Figure 2 .
Figure 2. superposition of human vitamin D receptor (gray, transparent) and reduced vitamin D receptor models (magenta, opaque).Sphere representation of 1,25(OH) 2 D 3 bound to VDR with stick representation of interacting amino acids.

Figure 4 .
Figure 4. total score and contributing terms for ligands/redVDR complexes.Panel A shows scores from the ligand geometry optimization procedure.Panel B shows scores from the molecular dynamics/ligand geometry optimization procedure.ΔH int -enthalpy of interaction in vacuo, -TΔS -entropy term, ΔG C -sum of desolvation and deformation enthalpies, -RTln(exp) -binding free energy with respect to the relative experimental binding constants.

table 1 . Energy and structural parameters for docked conformations of 1,25(OH) 2 D 3 and analogs.
Auto Dock binding energies (score, E, kcal/mol), hydrogen bond parameters (Å, deg) and RMS values (Å) between corresponding experimental and docked structures.RMS1 was calculated on heavy atoms.RMS2 is RMS1 with positional correction (rotation and translation were taken into account).