Identifying new liver X receptor alpha modulators and distinguishing between agonists and antagonists by crystal ligand pocket screening

Liver X receptor alpha (LXRα) is one of the LXR subtypes belonging to the nuclear receptor superfamily. As a ligand-activated transcription factor, LXRα contains a highly conversed DNA-binding domain and ligand-binding domain (LBD), regulating the expression of genes through binding to the region of the so-called LXR response element of the downstream genes. LXRα participates in the regulation of lipid metabolism and inflammation, and also cancer. Nonalcoholic fatty liver disease is a common liver disease associated with LXRα [1,2]. LXRα also showed oncogenic properties in cancer cells, such as gastric and ovarian cancer [3]. Thus, LXRα can be an essential pharmaceutical target to alter the biological process for different diseases [4]. LXRα modulators, by their transcriptional activities, can be divided into agonists or antagonists. Depending on different expressions or organ heterogenetic of LXRα, the agonist in principle can relieve nonalcoholic fatty liver disease symptoms and the antagonist can be used to treat cancer. LXRα is an attractive target for its potential functions, and new LXRα modulators have a chance to be developed into new drugs [5,6]. To meet pharmaceutical demands, we need to identify more LXRα modulators and also to distinguish between the modulator types, especially for the early phase of drug discovery [7]. Here, we selected the appropriate crystal pockets as two binding site models for efficient molecular docking screening to identify novel LXRα modulators and differentiate types of modulators.

To identify new LXRα modulators, we chose a suitable binding site model to conduct efficient structure-based molecular docking screening. The molecular dockability [8,9] of the x-ray crystal ligand-binding pocket was examined by the self-docking of crystal complex ligands. If the ligand after self-docking can correctly find its crystal pose, it can indicate the dual facts that the native binding pocket and the methodology of molecular docking used are both reliable. The similarity between the best docking pose and experimental crystal pose is calculated using the root-mean-square deviation (RMSD) [10]. In primary molecular docking screening, we also utilized two restrictions to reduce the number of purchased library compounds: consensus scoring [11] and Lipinski’s rule [12]. The details of the methodology and procedure for efficient virtual screening are described in the Methods section. We screened 25 candidates and further analyzed them by LXRα transactivation luciferase assay [13]. 20 hits were classified as significantly either enhancing or reducing transactivation.

To differentiate types of modulators, we needed to choose another suitable binding site model to distinguish the agonism of the hits by structure-based molecular docking screening. According to the standard theory of nuclear receptor transactivation mechanisms, an agonist binding in the core of the LBD with contact to the AF2 function site induces a conformational change which can introduce a binding site for coactivators in a groove formed majorly by Helix12, leading to effective transactivation [14–16]. Molecular docking to see whether a hit can dock into the LBD and whether the hit can also show contact with LXRα Helix12 His421, which may denote its agonism. If not, lack of contact with His421 indicates antagonism. For this simulation, we used the largest binding pocket, which can accommodate the greatest variety of LXRα ligand structures. We found that the x-ray LXRα complex (Protein Data Bank [PDB] ID: 3FC6) has the largest ligand-binding pocket. The methods used to decipher the agonism of hits are described below. The docking poses of a hit tending to have the hydrogen bond binding to LXRα Helix12 His421 strongly imply that the compound is an agonist, while lack of binding implies an antagonist. The best agonist hit and the best antagonist hit in this study were demonstrated to qualitatively illustrate such features.

Appropriate crystal-binding pockets as models of ligand-binding sites are crucial for the success of molecular docking simulations. Each binding cavity of the x-ray crystal LXRα complex can reflect a snapshot of a true dynamic ligand-binding pocket. Instead of computationally altering the pocket residues, we kept it frozen with the two selected x-ray coordinates of the ligand-binding environment for hit screening and to distinguish agonism. These two approaches may provide more efficient and more correct results for structure-based molecular modeling for both purposes.

Computational & experimental methods

Molecular docking

LigandFit [17] is a molecular docking method for accurately docking ligands into protein-binding sites and also for high-throughput hit screening. All molecular docking experiments using Discovery Studio 2.5 [18] were performed on an Asus computer with an Intel Core i7 4.00 GHz processor, running Windows 7. The energy grids of set binding sites can be further calculated by three different force fields, Dreiding [19], consistent force field (CFF) [20] and PLP1 [18] as implemented in LigandFit. Before molecular docking, the ligand is removed from the active site, leaving a binding site model for molecular simulation. We added hydrogen atoms to LXRα receptors and their ligands according to the default values of the CHARMM force field [21]. DockScore [17] was then used to estimate the interaction energy between calculated ligand poses and the LXRα ligand-binding pocket. We used the same methodology and molecular docking procedure as a common base to inspect the quality of x-ray crystal-binding pockets [9]. Ligand self-docking is an important metric for evaluating docking methods [22,23] by correctly simulating the crystal ligand pose. We use self-docking to examine the intrinsic quality of x-ray binding sites since our docking method and procedure were the same. We emphasize the validation of the resulting RMSD for measuring the accuracy of the docking methodology and the quality of the x-ray ligand pocket. Ligand self-docking has the same procedure as mentioned [9].

Utility-based selection of x-ray binding sites

Binding site for molecular docking was defined by the occupied space of the cocrystal ligand, which was set by expanding the ligand hydrogen radius to 2.0 Å and heavy ligand atom to 2.5 Å (default settings) until reaching residues. We collected 19 x-ray crystal protein chains of human LXRα cocomplexed with ligands [14,24–29] deposited in the PDB (more than one decade) [30]. The same definition can also calculate the size of each binding site. The best x-ray binding site for new modulator screening was examined by the self-docking strategy, while the x-ray binding site for distinguishing agonism used the largest one.

New modulator screening & consensus scoring

To be as close to reality as possible, we set all the crystal protein residues as still. All crystal waters were also removed since the orientation of water is unknown. A maximum of ten poses was set to be saved after molecular docking. Then, for each docking pose, 11 scoring functions, LigScore1 [17], LigScore2 [17], DockScore [17], PMF [31], PMF04 [32], Ludi1 [33,34], Ludi2 [33,34], Ludi3 [35], PLP1 [18,36], PLP2 [18,36] and Jain [37], were used for consensus scoring. Consensus scoring was applied as follows: each library compound can have up to ten after-dock poses, and each pose has 11 scores. If a pose’s score in one of the 11 scoring functions was within the top 6% of total library docking poses, by definition, the pose gained one consensus point. An in-house Ruby (www.ruby-lang.org) programming script for docking-pose-based consensus scoring was used to calculate consensus points.

Lipinski’s rule of five

In drug discovery, desired molecules can fulfill particular criteria for further drug development. Lipinski’s rule of five (ROV) [12] is a traditional way to address the drug-likeness of the organic compounds. ROV rules include no more than ten hydrogen bond acceptors, no more than five hydrogen bond donors, molecular weight (mass) lower than 500 Da and partition coefficient (AlogP) smaller than five. As improvements for the drug-likeness of organic compounds, Veber and coworkers [38] showed that the polar surface area (PSA) [39] and the number of rotatable bonds (RBs) can be discriminating molecular descriptors. Two additional criteria can be PSA <140 Ų and RBs ≤10. Those compounds fulfilling Veber’s rules are expected to have good oral availability. To further lower time and cost, we wanted to identify the best hits successfully passing ROV and the two additional criteria. The related molecular properties were calculated by the Discovery Studio 2.5 module for molecular properties.

Compound library

The Maybridge Screening collection, with >57,000 compounds, was the compound library used for screening. All Maybridge compounds were optimized within Pybel [40,41] using the MMFF94 forcefield until the structure’s energy reached the RMSD gradient of 0.05 kcal/mol/Å [42]. The partial charges and the optimized structures of chemical compounds were automatically set and saved in mol file format for further hit screening of molecular docking.

LXRα transactivation luciferase assay

LXRα transactivation tests including Huh7 cells were used for the luciferase assay. Transfections were done in Huh7 cells using Lipofectamine Plus (Life Technologies, Inc., MD, USA) according to the manufacturer’s instructions. In each transactivation assay, the following plasmids were involved: a reporter containing five Gal4 binding sites (0.2 µg), an expression plasmid for LXRα LBD fusion protein (0.6 µg) and 60 ng SV40-Relina plasmid. Pure compounds were added 12 h after cell transfection was done for 1 day. Cells were collected and lysed for luciferase assays following the manufacturer’s instructions. L783483, being a pan-agonist for the nuclear receptors, was used as an agonist standard. The LXRα ligand transactivation activities of 25 screened Maybridge compounds were tested with 10 μM as the primary screening assay.

Results & discussions

Best x-ray binding site for molecular docking screening

The quality of the binding sites was our first concern for highly accurate molecular docking virtual screening for LXRα. X-ray crystallography can resolve high-resolution structures at the atomic level and structural details, but the intrinsic resolution of x-ray crystallography is still not enough for accurate structure-based drug design. Theoretically, a 2 Å deviation from reality at a critical binding-site residue can lead to failures in drug design, since the length of a single bond is <2 Å. We can see that the average resolution of the collected LXRα PDB complexes is 2.42 Å (Table 1, resolution column) [14,24–29]. Such a departure from reality, if it happens in the binding site, can severely decrease the accuracy of molecular docking, and thus, the efficiency of virtual screening. Therefore, a further inspection of the ‘dockable nature’ [8,9] of the ligand pockets of each crystal complex is essential for molecular docking. To raise our confidence for the crystal-binding pockets of LXRα for molecular docking, we used an additional strategy, self-docking, in which a geometry-randomized crystal ligand docks back into the natural crystal-binding cavity. The degree of regaining the ligand-bound pose is initially one of the central criteria for the correctness of the docking method [43,44], but here we used it to check the reliability and suitability of crystal-binding sites for hit screening. The shape overlap between the x-ray pose and the best ligand-docking pose was calculated by RMSD. The 19 results of self-docking are listed by the energy grids of the three different force fields, Dreiding, CFF and PLP1, in Table 1, where the RMSD ranged from 0.37 to 2.57 Å, from 0.40 to 2.57 Å and from 0.47 to 12.66 Å, respectively, as shown in the three violin distribution graphs (Figure 1A). We can see that better crystal protein resolution did not ensure lower RMSD values after self-docking. The best docking pose of the 3IPU ligand (Dreiding) was similar to the native bound crystal pose, with an RMSD of 0.37 Å (Figure 1B). Thus, the 3IPU ligand-binding site was considered to be the best x-ray binding site for the new molecular docking modulator screening of LXRα. Lastly, we stress here that all the protein residue side chains were not allowed to move to remain close to reality as possible, and also that crystal waters were also removed since the orientation of water is unknown in x-ray crystallography. In this case, there is no water existing in this binding pocket. We believe that each x-ray binding site reflects a snapshot of dynamic ligand–pocket interactions. Thus, we left the selected x-ray coordinates unchanged, rather than computationally altering the ligand-binding surroundings of the pocket. So, in this study, we excluded the method of molecular dynamics.

Possible LXRα modulators by screening & consensus scoring

The best 25 hit compounds were purchased for in vitro reporter gene assays. The chemical structures are shown in Figure 2. Each hit was picked out by the following sequence: when the best x-ray binding site was located, we did molecular docking to each of >57,000 library compounds. This step generated about 167,000 docking poses with 11 scoring functions, of which 228 screening compounds had consensus points >7. Ligand–receptor binding may involve a variety of mathematical forms for noncovalent interactions, presenting differently in diverse ligand–receptor interactions [45–47]. As the literature shows [22,23], no single scoring function can handle all types of binding pockets equally well. To ensure finding the best hit, we did consensus scoring [45–47] for all docking poses. If the specific score for docking poses was within the top 6% of total docking poses it gained one consensus point. We took all docking poses as the total population rather than only the best docking poses from the library compound, and thus, the consensus scoring was based on docking poses.

Moreover, the selected compounds have different structures than the parent crystal ligand. Scaffold hopping [48] seemed to be easily achieved in this molecular docking screening. Note that we chose hits by the prediction of ligand–receptor interaction for our concept validation. However, to reduce the time and cost of the luciferase assay, we further considered the ROV, PSA and number of RBs for the possible drug development of these hits. The related ROV values are listed in Supplemental Table 1. At last, we purchased 25 compounds for the in vitro reporter gene assay.

Novel LXRα modulators by in vitro reporter gene assay

The best 25 hit compounds were evaluated by in vitro LXRα luciferase reporter gene assay. The luciferase activities for all hit structures are shown in Figure 3. Luciferase assays are a typical cell-based experimental method used to measure the transactivation activity of LXRα. The 25 hits and L783483 [49] were assayed in transfected Huh7 cells at 10 μM. The control was the normal luciferase activity of the system when no compound was added. L783483 is a pan-agonist for nuclear receptors, here used as an agonist standard for comparison. The outcome of luciferase chemiluminescence reveals the degree of LXRα transactivation, showing whether a hit activates or inhibits it. More than half of 20 hits show a significant influence on LXRα transactivation. Then we further classified the identified hits into five groups. In group 1, six agonists (1, 9, 11, 12, 22 and 24) showed higher LXRα transactivation activity of at least 20% more than control. Notably, among group 1, the four agonists (1, 9, 11 and 22) had LXRα transactivation higher than or similar to that of L783483. Group 2 had the 14 antagonists (2–8, 13, 14, 16, 19, 20, 23 and 25), showing the LXRα transactivation activity less than 80% of blank control. We also noticed that two antagonists (2 and 5) showed significant antagonism less than 20% of control and most antagonists show less than 50% of control. Group 3 had five hits (10, 15, 17, 18 and 21) showing no significant alteration of LXRα transactivation activity and were similar (≤ ±20%) to blank control. Overall, 80% of our purchased LXRα modulators demonstrated molecular docking of selected x-ray binding sites, suggesting that this method could provide an extremely efficient way of virtual screening.

Distinguishing agonism

According to the recognized concepts of nuclear receptor transactivation mechanisms, a ligand is necessary to bind to the AF2 function site for the activation of a specific nuclear receptor [14–16]. In LXRα (as well as other nuclear receptors), we can also see that a full agonist always exhibits a feature whereby the head of the ligand structure body has a direct hydrogen bond link with the AF2 Helix12 His421 in x-ray crystallography images. We used this criterion to see if our molecular docking simulation could differentiate the agonism of a virtual screening hit. We believe any computational alteration of the protein backbone structure may cause some amount of uncertainty. A specific binding site can be seen as a snapshot of the dynamic reality of LXRα complex. We tried to avoid unnecessary modifications to the original structure coordinates revealed by x-ray crystallography. To allow the maximum docking poses, we took the largest binding pocket, which may accommodate the greatest variety of LXRα ligand structures. The volumes of different binding sites, as well as the largest-size ligand-binding pocket (PDB ID: 3FC6), are listed in Table 1 (binding site columns) and illustrated in Figure 4A and B. Our aim was to see if hits could dock into the LBD and whether the hit furthered show a hydrogen bond link with 3FC6 LXRα AF2 His421 [14–16] or not. As a result of all simulation observations, the best antagonist (compound 5) and the best agonist (compound 11) in our study were taken as examples. The docking pose of compound 5 showed multiple definite binding interactions: a strong π–π interaction with Trp443 (head of the ligand) and hydrogen bonds with Asn225 and Arg305 (tail of ligand) as shown in Figure 5A. The docking pose of compound 11 also demonstrated multiple definite binding interactions: strong π–cation interactions with Arg305 (tail of ligand) and hydrogen bonds with Ser264 (body of ligand) and, most importantly, with His421 (head of the ligand) as shown in Figure 5B. By observing plentiful docking poses of nuclear receptor ligands including LXRα ligands, we can safely draw a qualitative conclusion that either a potent binding nuclear receptor agonist or antagonist tends to have multiple significant interactions in the binding site and that the distinguishing feature is that the full agonist must have the hydrogen bond binding to Helix12 His421. Furthermore, the particular structural features of the strong LXRα agonists are having strong hydrogen acceptors on the heads of agonists, which can fit spatially into the Helix12 His421 vicinity surroundings. For example, the trifluoromethyl group in Y1, the carbonyl group in Y9 and Y22 and the carboxylic acid in Y11 serve as such hydrogen acceptors. In the molecular docking simulation, the greater frequency of this hydrogen bond link in hit docking poses indicates that the hit is a full agonist (a strong binder with Helix12 His421). In contrast, in a strong binding antagonist hit it is hard to see any such hydrogen bond link in any docking pose (a strong binder with no Helix12 His421 link).

Conclusion

Modulators of LXRα are of high pharmacological interests as LXRα regulates fatty acid metabolism, inflammatory processes and even related to cancer. Based on the suitably selected experimental native binding sites of x-ray receptor complexes for molecular docking, we aim both to identify the new LXRα modulators and to recognize the distinguishable feature between agonists and antagonists. Our two approaches, ligand self-dock and largest-cavity-size searching, purposely located two appropriate ligand-binding site models to reach the following two aims: for identifying new modulators, after sorting all the docking poses of a chemical library of 57,000 compounds, out of 25 purchased compounds, in vitro reporter gene assay confirmed up to 20 agonists or antagonists; and for distinguishing types of the modulators, the results of docking pose indicates a feature that at least one pose, out of ten poses, having one hydrogen bond to LXRα Helix12 His421 denotes its agonism. Thus, the main structural features for the LXRα agonists are having strong hydrogen acceptors on the heads of agonists fitting to the Helix12 His421 vicinity surroundings. Overall, the utility-based selections of x-ray binding sites showed that molecular docking not only can screen new active modulators but also can differentiate between agonism and antagonism.

Future perspective

In the future, new LXRα modulators will have a chance to be developed into new drugs. Appropriate crystal-binding pockets as models of ligand-binding sites are important for the success of molecular docking simulations. Each binding cavity of the x-ray crystal LXRα complex can reflect only a snapshot of a true dynamic ligand-binding pocket. The two selected x-ray coordinates of the ligand-binding environment were used for hit screening and to distinguish agonism. These two approaches may provide more efficient and more correct results for structure-based molecular modeling for both purposes. The full agonist must have the hydrogen bond binding to Helix12 His421, is also an important judgment for virtual LXRα modulators development. We conclude that computationally assessed qualities of crystal pockets are an important checkpoint that can help discover structurally novel nuclear receptor modulators and cost less.