|Home | About | Journals | Submit | Contact Us | Français|
Understanding and identifying the receptor subtype selectivity of a ligand is an important issue in the field of drug discovery. Using a combination of classical molecular mechanics and quantum mechanical calculations, this report assesses the receptor subtype selectivity for the human retinoid X receptor (hRXR) and retinoic acid receptor (hRAR) ligand‐binding domains (LBDs) complexed with retinoid ligands. The calculated energies show good correlation with the experimentally reported binding affinities. The technique proposed here is a promising method as it reveals the origin of the receptor subtype selectivity of selective ligands.
Retinoid X receptors (RXRs) and retinoic acid receptors (RARs) are class 1 and class 2 nuclear receptors (NRs), respectively 1. Both receptors exhibit α, β, and γ subtypes. RXRs and RARs form heterodimers, which control crucial biological events such as cell differentiation and proliferation, morphogenesis, and homeostasis. All‐trans retinoic acid (ATRA; Fig. Fig.1)1) is a natural ligand of RARs. Its geometrical isomer, 9‐cis retinoic acid (9cRA; Fig. Fig.1),1), is thought to be a natural ligand of RXRs, although 9cRA also strongly binds to RARs. Recently, molecular mechanics docking simulations were used to investigate the binding of ATRA to RXRs, and this study suggested that ATRA could act as an endogenous ligand of RXRs 2. Moreover, the helix H3 three‐point initial‐binding hypothesis of ligand in the ligand‐binding domains (LBDs) of the NR superfamily was proposed, and the driving forces behind ligand entry into the NR LBDs were discovered, leading to the successful use of molecular dynamics (MD) simulations to understand the structural transition between the apo‐form and holo‐form of both hRXRα and hRARγ LBDs 3. Furthermore, the local motifs that determine the canonical fold of NR LBDs were discovered, and the agonism and antagonism of the NR superfamily at the electron level were proposed 1.
In this study, the origin of receptor subtype ligand selectivity was investigated for hRXRs and hRARs with their natural ligands (9cRA and ATRA; Fig. Fig.1)1) and the synthetic ligands LGD1069 4 (Bexarotene) and Am80 5 (Tamibarotene) (Fig. (Fig.1).1). LGD1069 is a representative RXR‐selective ligand and Am80 is a representative RAR‐selective ligand. Although some studies analyzed the binding affinities of different compounds for only one receptor by using quantum mechanics calculations, investigation of the correlation among receptor subtypes with different compounds has thus far not been reported in the literature. Here, a promising method for understanding the origin of receptor subtype selectivity with these selective ligands is demonstrated. Specifically, second‐order Møller–Plesset perturbation (MP2) theory calculations, in conjunction with ONIOM method 6, were undertaken to evaluate the binding free energies (ΔG bind) between the receptors and ligands. In addition, MP2 theory calculations using the fragment molecular orbital (FMO) method 7 were performed on the entire receptor–ligand system to analyze the interaction energies between the receptors and ligands.
With the exception of the hRXRγ LBD, the three‐dimensional (3D) structures of the liganded LBDs of hRXRs and hRARs in the agonist conformation have been obtained using X‐ray crystal structure analysis. In this study, the initial structures were prepared using Homology Modeling Professional for hyperchem 8, 9 (each structure contains a coactivator fragment and water molecules), as described in a previous report 2. The 3D structure of hRXRγ LBD in the agonist conformation was prepared using a homology modeling technique, as previously described 2. Docking simulations (AMBER99 force field) using Docking Study with hyperchem 9, 10 were performed under biomacromolecule‐rigid and ligand‐flexible conditions, as previously described 2. With the exception of the water molecules conserved in the ligand‐binding site, water molecules as well as the coactivator fragment were removed during this stage. A Gaussian 11 job file for the most stable complex was automatically prepared using ONIOM Interface for Receptor 12 integrated into the Homology Modeling Professional for hyperchem software. The three‐layer ONIOM calculations (B3LYP/6‐31G*:AM1:AMBER) were carried out for all combinations of LBDs with the ligands; the ligand was defined as the high layer (B3LYP/6‐31G*), the amino acid residues and conserved water molecules positioned within 4 Å of the heavy atoms of the ligand were defined as the medium layer (AM1), and the remaining structures were defined as the low layer (AMBER). In these calculations, the structures of both high and medium layers were fully optimized, whereas only hydrogen atom positions were optimized for the low layer. The binding free energies (ΔG bind) in the gas phase at 298.15 K were obtained from the single‐point frequency analysis for the converged complex (using two‐layer ONIOM calculations; MP2/6‐31G:AMBER) and the isolated receptor (AMBER) and ligand (MP2/6‐31G). For the converged structures, FMO calculations were performed at the MP2 level of theory using the ABINIT‐MP program 13. The job files were prepared using BioStation Viewer 14. Interfragment interaction energies (IFIE) for the ligands were obtained using single‐point calculations at the MP2/6‐31G level of theory.
It has recently been reported that the interaction energies obtained from biomacromolecule‐rigid and ligand‐flexible docking simulations using a classical molecular mechanics force field showed excellent correlations with the experimental binding affinities of 9cRA and ATRA 2. The docking simulations for structurally different synthetic ligands have now been investigated under the same conditions (Table 1 and Fig. Fig.2).2). Since the RXR subtype‐selective ligand is not obtained experimentally, LGD1069 was chosen in this study. On the other hand, Am80 exhibits RARα and RARβ subtype selectivity. From the results, it is apparent that the correlations between the calculated interaction energies and experimental binding affinities are not as strong as previously reported, although there is still some level of agreement (Fig. (Fig.2;2; bottom). It seems that the structure of the ligand‐binding site is strongly dependent on the crystal structure with the natural ligands, that is, 9cRA and ATRA. Biomacromolecule‐ and ligand‐flexible docking simulations were also performed but the correlations between the calculated interaction energies and experimental binding affinities were not improved (see Table S1 and Fig. S1).
The structures of the ligand‐binding site and complexed ligand obtained from the docking simulations were further optimized using QM/MM ONIOM calculations at the B3LYP/6‐31G*:AM1:AMBER scheme. For the optimized structures, the single‐point two‐layer ONIOM calculations (MP2/6‐31G:AMBER) for the complexes and the single‐point calculations of the isolated receptors (AMBER) and ligands (MP2/6‐31G) were performed to obtain the ΔG bind(calc) values. The ΔG bind(calc) values in the gas phase showed good correlation, with a correlation coefficient (R) of 0.43 (Fig. (Fig.2;2; top left). As is apparent from the right‐hand side of Fig. Fig.22 and from Table 1, with the exception of hRXRγ and hRARα LBDs, the correlation coefficient for the LBDs was high (0.73). This underlines the suitability of the method used in this study (when only 9cRA and ATRA were plotted for the four receptors, R was 0.94, thus strengthening our previous conclusions 2). Only the 3D structure of hRXRγ LBD was prepared by homology modeling, and this may be a reason for the low correlation observed for hRXRγ LBD; however, the reason for the low correlation for hRARα LBD is unclear (it may be attributed to allosteric effect from the antagonistic heterodimer partner of the original structure, as shown in Table S2).
Single‐point calculations for the whole system were then performed using the FMO method at the MP2/6‐31G level of theory. Under these conditions, all SCF calculations converged and the resulting IFIE values were reasonable. The middle of Fig. Fig.22 shows the correlations between the obtained IFIE values and experimental ΔG bind values for the six receptor subtypes studied. The IFIE values showed poor correlation, with a correlation coefficient of 0.29 (Fig. (Fig.2;2; middle left). However, as is apparent from Table 1, with the exception of hRXRγ LBD, the correlation coefficient for the LBDs improved to 0.41 (at the middle right of Fig. Fig.2,2, R was 0.40 for the four receptors, with the exception of hRXRγ and hRARα LBDs). It seems that the geometrical optimization calculations at the MP2 level of theory with larger basis set will be needed for the FMO method for the whole system.
Table 2 shows the amino acid residues of the ligand‐binding sites of hRXRs and hRARs. The residues involved in hRXR ligand‐binding sites are identical regardless of subtype, underscoring the difficulty faced during the assessment of receptor subtype selectivity. In contrast, one (α and β), two (β and γ), or three (α and γ) residue(s) of the ligand‐binding sites of hRARs differ(s) from each other. The docking simulations support the experimental observation that Am80 does not bind to hRARγ as well as hRXRs (Table 1) and suggest that Am80 is able to recognize the structure of the ligand‐binding site. Although the sequences of RXR and RAR LBDs are moderately conserved, with approximately 30% identity, 9cRA and ATRA can bind to both receptors. For example, there are only three identical ligand‐binding residues in the LBDs of hRXRα and hRARα: Leu309 and Leu269, Ile310 and Ile270, and Arg316 and Arg276 (Table 2). The ability of ATRA and 9cRA to bind to both receptors despite the different residues in the binding sites could be attributed to the conformational flexibility of these naturally occurring retinoids 2, which is absent in the synthetic retinoids (Fig. (Fig.3).3). As shown in Fig. Fig.3,3, the structure of LGD1069 is more flexible than that of Am80; thus, LGD1069 can bind to both hRXR and hRAR LBDs, and its agonistic potency for hRXRs is in the range of at most one or two order(s) of magnitude greater than that for hRARs 15. Therefore, the receptor subtype selectivity of the ligands could not be accounted for by considering only the differences in the binding residues. Rather, it appears that the binding affinities of the ligands can instead be estimated from the ΔG bind(calc) values obtained from quantum mechanical MP2 calculations for the entire receptor–ligand system.
This study demonstrates that receptor subtype selectivity can be predicted from energies calculated for an entire receptor–ligand system using MP2 quantum mechanical calculations. Quantum mechanical geometrical optimization of the complexed ligand and the surrounding region was essential for obtaining reliable values. The conventional method, wherein the initial structure is prepared using classical molecular mechanics calculations such as geometry optimizations and MD equilibrations, was insufficient to evaluate the binding affinities.
In conclusion, the method described in this study proved useful for assessing the receptor subtype selectivity of ligands. The origin of the receptor subtype selectivity of ligands highlighted by these results remains the subject of further investigation.
MT, KS, and HK designed the experiments. MT performed the experiments, analyzed the data, and wrote the paper. All authors approved the manuscript.
Table S1. Experimental ΔG bind values and calculated interaction energies of naturally occurring and synthetic retinoids.
Table S2. Receptor information used in this study.
Fig. S1. Correlations between ΔG bind(exp) and interaction energies of the most stable complex obtained from the docking simulations (the biomacromolecule‐ and ligand‐flexible conditions) for the binding of the α, β, and γ subtypes of hRXR and hRAR LBDs with ATRA, 9cRA, Am80, and LGD1069.
MT is deeply grateful to Dr. Neil S. Ostlund of Hypercube, Inc. for his help during the development of Homology Modeling Professional for hyperchem and Docking Study with hyperchem. MT also thanks Gaussian, Inc. for supplying the Gaussian program.