Identifying selective agonists targeting LXRβ from terpene compounds of alismatis rhizoma

Hyperlipidemia is thought of as an important contributor to coronary disease, diabetes, and fatty liver. Liver X receptor β (LXRβ) was considered as a validated target for hyperlipidemia therapy due to its role in regulating cholesterol homeostasis and immunity. However, many current drugs applied in clinics are not selectively targeting LXRβ, and they can also activate LXRα which activates SREBP-1c that worked as an activator of lipogenic genes. Therefore, exploiting agonists selectively targeting LXRβ is urgent. Here, computational tools were used to screen potential agonists selectively targeting LXRβ from 112 terpenes of alismatis rhizoma. Firstly, a structural analysis between selective and nonselective agonists was used to explore key residues of selective binding with LXRβ. Our data indicated that Phe271, Ser278, Met312, His435, and Trp457 were important to compounds binding with LXRβ, suggesting that engaging ligand interaction with these residues may provide directions for the development of ligands with improved selective profiles. Then, ADMET analysis, molecular docking, MD simulations, and calculation of binding free energy and its decomposition were executed to screen the agonists whose bioactivity was favorable from 112 terpenes of alismatis rhizoma. We found that two triterpenes 16-hydroxy-alisol B 23-acetate and alisol M 23-acetate showed favorable ADMET properties and high binding affinity against LXRβ. These compounds could be considered as promising selective agonists targeting LXRβ. Our work provides an alternative strategy for screening agonists selectively targeting LXRβ from alismatis rhizoma for hyperlipidemia disease treatment.


Introduction
Hyperlipidemia is characterized by increased triglyceride (TG), total cholesterol (TC), and lower-density lipoprotein cholesterol (LDL-C) level and decreased high-density lipoprotein cholesterol (HDL-C) level in the blood [1]. Hyperlipidemia is considered an essential contributor to not only coronary arteries and coronary heart disease but also hypertension and diabetes. A 1% increase in plasma cholesterol level will enhance the risk of coronary events to almost 3% [2]. In addition, excessive accumulation of TG caused by hypercholesterolemia in the liver is responsible for nonalcoholic fatty liver disease [3]. According to the data model predicted by the World Health Organization (WHO), it will be up to 40% of all deaths that are related to cardiovascular disease by 2030 [4]. Hyperlipidemias are divided into six types by WHO based on whether levels of lipids and lipoproteins are increased [5]. The most common frequent hyperlipidemias are combined family hyperlipidemia and polygenic hypercholesterolemia [6]. Causes of hyperlipidemia are very complicated, which are relevant to environmental factors and genetic factors [5]. In recent years, due to people's diet and unhealthy lifestyle, the risk of hyperlipidemia has greatly increased [7].
One of the causes of hyperlipidemia is a high level of TC, which is tightly regulated by liver X receptor (LXR) [1,8]. LXR is a nuclear receptor with two isoforms, LXRα and LXRβ. LXRα is mainly distributed in the liver, intestines, adipose tissue, and macrophages while LXRβ can be found in various systemic tissues [9,10]. In vivo, LXRα and LXRβ all have the ability to promote cholesterol efflux by Chuanjiong Lin and Jianzong Li contributed equally to this work. modulating transcriptional activation of genes [11]. LXRα is the upstream activating gene of sterol response element binding protein-1c (SREBP-1c). SREBP-1c is a transcription factor to regulate the expression of various lipogenic genes, such as acetyl-CoA carboxylase (ACC) and fatty acid synthase (FAS), sequentially causing an undesirable increase of TC levels in the plasma and liver [12,13]. However, different from LXRα, LXRβ is considered to have less effect on the lipid synthesis pathway [14]. After forming permissive heterodimers with retinoid X receptor (RXR), LXR functions as lowering blood lipid by binding to LXR response element (LRE), which has been found in the target gene promoter region [15]. Moreover, LXR ligands induce the expression of the members of the ATP binding cassette family of lipid transporter genes, such as ABCA1, ABCG, and ABCA2. In addition, LXR ligands can also promote transcription of lipid synthesis genes, which are represented by SCD, ACC, and FASN [16].
Natural products are potential resources for the discovery of small molecules modulating target activity. Alismatis rhizoma, the dried tuber of Alisma orientale (Sam.) Juz. (Alismataceae), is a kind of traditional Chinese medicine (TCM) widely distributed in China, Russia, Japan, Mongolia, and North India [17]. In East Asia, alismatis rhizoma has been used to treat many diseases, such as hyperlipidemia, diabetes, hypertension, and urological diseases [18]. It's worth noting that the most effective constituents of alismatis rhizoma are attributed to terpenes [19]. Terpenes isolated from alismatis rhizoma have been proved to enhance transactivation for farnesoid X receptors (FXR), which belongs to the nuclear receptor superfamily along with LXR [20]. These terpenes are mainly composed of sesquiterpenes, diterpenes, and triterpenes [17]. Among them, protostane-type triterpenes are considered a major contributor to pharmacodynamic and bioactive compounds of alismatis rhizoma. Moreover, protostane-type triterpenes are also thought to be chemotaxonomic markers of alismatis rhizoma [21]. Therefore protostane-type triterpenes have drawn increasing attention in recent years because of their unique biological functions in treating chronic diseases, such as hyperlipidemia, diabetes, hypercholesterolemia, and cancer [22]. Alisol A, alisol B, and their acetylation and alisol C 23-acetate are major components of protostane-type triterpenes [23]. In addition, alisol B 23-acetate, alisol C 23-acetate, alisol A, and its acetylation have shown a superb effect to reduce lipid, and alisol B and alisol B 23-acetate are involved in autophagy, apoptosis, and endoplasmic reticulum stress in cancer cells [21].
The study aims to screen selective agonists targeting LXRβ from terpenes of alismatis rhizoma. Firstly, structures of LXRα and LXRβ were analyzed, and differences of protein sequences and residues which are highly stabilized in sitedirected mutation were evaluated. Then, we constructed a verification standard to verify agonists whether possessed LXRβ selectivity by analyzing the differences of the binding mode of developed agonists between LXRα and LXRβ. A dataset of 112 terpenes from alismatis rhizoma was constructed as the raw material of the study. SwissADME and ADMETlab were used to predict ADMET properties in order to eliminate compounds whose pharmacokinetic parameters were unfavorable. Then, molecular docking was executed by AutoDock Vina for the sake of preliminarily evaluating binding affinities of ligands with LXRα and LXRβ. In order to verify the reliability of the docking model, DecoyFinder 2.0 was applied to construct a set of decoy molecules for a group of ligands which were binding closely with LXRβ. The receiver operating characteristic (ROC) curves of the docking score of decoy molecules were plotted and area under curve (AUC) values were calculated. Moreover, a molecular dynamic simulation was conducted using Gromacs to test the dynamic binding mode of ligands and receptors during a specified simulation time. Furthermore, binding free energy and its decomposition were calculated by g_mmpbsa package to predict binding affinity between ligands and proteins. In addition, the binding model analysis was executed by Ligplot + and PyMOL. This study serves as an alternative strategy for screening selective agonists targeting hypolipidemic target from traditional Chinese medicine.

Structure collection and preparation
Alismatis rhizoma has been proved to be an effective herb for treating hyperlipidemia, and terpenes were its main active components. Therefore, 112 triterpenes and sesquiterpenes from alismatis rhizoma were collected from databases and literature [24]. In order to identify key residues that play a selective role in ligand and LXRβ binding, 22 known selective or nonselective agonists were collected from different databases. The structures of LXRα (PDB code:5HJS, X-ray diffraction resolution with 1.72 Å) and LXRβ (PDB code:6S4T, X-ray diffraction resolution with 2.00 Å) were retrieved from the RCSB protein data bank [25]. open = 10; gap extend = 0.5; end gap penalty, false; end gap open = 10; and end gap extend = 1. Click "Submit" to run the alignment, and Smith-Waterman algorithm was applied during the process. The result was output in several seconds, and the length, similarity, gaps, and identity were shown as a list.

Alanine scanning
In order to estimate the stability effect of all possible mutations of each residue of LXR, FoldX 4.0 [28] was introduced to calculate the change of binding free energy of LXR between wild type and mutations. Each residue in mutations was changed into standard α-amino acids using the BuildModel function after applying the automatic RepairPDB function. The temperature was set to 298 K, and other options were set to default. The result for respective value was shown as ΔStability, which means binding free energy of the alanine mutations minus that of wild type (ΔStability = ΔG mut -ΔG wt ).

ADMET analysis
Suitable ADMET properties were the first hurdle for most drugs, and pharmacokinetic parameters (absorption, distribution, metabolism, excretion and toxicity, ADMET) were significant contributors to the discovery of drug candidates [29]. Therefore, SwissADME [30] and ADMETlab [31] were applied to predict ADMET properties. Here, compounds with high absorption and low toxicity and that met Lipinski's rule of five (RO5) were considered for further study. The details of RO5 are as follows: molecular weights (MW) ≤ 500; the calculated lipophilicity (log P) ≤ 5; the number of hydrogen bond acceptors (HBA) ≤ 10; the number of hydrogen bond donors (HBD) ≤ 5; and the number of rotatable bonds (RB) < 5 [32,33].

Molecular docking model validation
In order to identify the reliability of the docking model, the receiver operating characteristic (ROC) curves of docking functions were analyzed. DecoyFinder 2.0 was utilized to establish a set of decoy molecules for ligands which were obtained from the ChEMBL database according to affinities (K i < 3000) [34]. Parameters of DecoyFinder 2.0 were set as follows [35]: active ligand vs. decoy Tanimoto threshold ≤ 0.75, decoy vs. decoy Tanimoto threshold < 0.90, hydrogen bond acceptors ± 2, hydrogen bond donors ± 1, molecular weight ± 25 Da, rotational bonds ± 1, octanol-water partition coefficient (log P) ±1.00, maximum standard deviations away from active = 1.00, maximum decoys per active ligand = 36, and minimum decoys per active ligand = 36. Then, the ligands and decoys were regarded as the input set of molecular docking, and docking results were shown as the ROC curves, and the area under the curves (AUC) were calculated by the IBM SPSS program. The parameters of molecular docking will be described in the next paragraph.

Molecular docking
After verifying the druggability of molecules and the reliability of the docking model, molecular docking was performed using AutoDock Vina [36] as a tool of virtual screening. Firstly, the UCSF chimera molecular visualization application [37] was used to separate ligands and receptors and eliminate water on receptors. Then, AutoDock Tools (ADT, version 1.5.6) [38] was applied to add hydrogen atoms, Gasteiger charges, and atom types to receptors and the number of torsions for ligands. After that, ligands and receptors were output as pdbqt files which can be recognized by the AutoDock Vina program. A grid box was generated by ADT for AD4 atom types and centered on the ligand. The size of box was set as follows: size_x = 15 Å, size_y = 18.75 Å, and size_z = 15 Å. And the center of box was set as follows: center_x = 16.68, center_y = 10.492, and center_z = 15.764. The Vina score was calculated by AutoDock Vina.
LeDock score, Grid score, and Amber score were also used to test binding affinity between ligands and receptors. LeDock score was obtained from LeDock [39]. Grid score and Amber score were calculated from UCSF Dock6 [40]. First of all, the preparation of docking models of the three functions was the same. Water removal of protein and energy minimization of ligands should be executed by UCSF chimera [37]. Then, operation steps were different for the three functions. For LeDock, GetBox plugin in PyMOL program was introduced to set an enclosing box which is centered on ligands. Then, proteins and ligands files were input into LeDock to dock. The protein files were set as the input of LePro, and hydrogens were added automatically while the ligand files were input of LeDock. For Dock 6, the docking process was executed according to the docking user manual, and the Grid score and Amber score were applied.

Molecular dynamic simulations
In order to probe the flexibility and structural stability of screened compound candidates, we chose the output files of molecular docking and the water-free receptors as raw materials of molecular dynamic (MD) simulations. MD simulations were executed by Gromacs 5.1.4 [41] for a period of 20 ns. First of all, ligands and receptors were transformed into pdb2gmx format and mol2 format, respectively. To convert pdb files to pdb2gmx format, the AMBER99SB force field and the tip3p water model were applied to generate protein structure in the topology and gro files which contain protein sequences. The General Amber force field (GAFF) [42] was used for ligands, and partial charges of ligands were determined by the Antechamber program [43]. To check the missing force field parameters and amend the field files, the parmchk program was employed [43]. Antechamber Python parser interface (ACPYPE) tools were used to generate ligand structure in the topology and ligand sequences [44]. Then a box with a distance of 1.0 nm from the protein to the edge of boundary was defined. Na + and Cl − were added at a concentration of 0.1 mol·L −1 to ensure an electroneutral ionic environment, equivalent to one Na + and one Cl − for every 600 water molecules. In order to minimize the total energy of the system, adjusting wrong bonds and angles of stretching or compression, the energy minimization algorithm was introduced to stabilize the system [45]. After that, under the guidance of index files and with a limit of proteins and ligands, simulations of 100 ps NVT (constant number of particles, volume, and temperature) and 100 ps NPT (constant number of particles, pressure, and temperature) equilibrations were executed at constant temperature (300 K) and pressure (1 atm). MD simulations at 20 ns were performed with a timestep of 2 fs at constant pressure and temperature. The linear constraint solver (LINCS) algorithm was used to constrain all covalent bond lengths during the process [46], and the particle mesh Ewald (PME) method was used to calculate the long-range electrostatic interactions [47], while the SETTLE algorithm is applied for water molecules [48]. Finally, gmx_rmsd and gmx_rmsf in Gromacs 5.1.4 were employed for root-mean-square deviation (RMSD) and rootmean-square fluctuation (RMSF) analysis on the entire 20 ns track [41]. Then, xmgrace was used to draw RMSD and RMSF results in order to visualize the output of RMSD and RMSF analysis.

Calculation of MM/PBSA binding free energy
In order to explore the binding affinity between ligands and proteins, the Molecular Mechanics Poisson-Boltzmann Surface Area (MM/PBSA) [49] was introduced to calculate the binding free energy of ligand-receptor complexes by the g_mmpbsa tool. The last 2-ns trajectory of MD results was extracted as a raw material of binding free energy. In this approach, binding free energy consists of the following components: Here, G bind represented the binding free energy. G MM was the molecular mechanic energy of gas phase, and G sol was the solvation free energy. G int was the bonded interaction between molecules, including bond, angle, and dihedral energy. G vdw was Van der Waals' force, and G ele was the electrostatic interaction. G sol was composed of G polar and G apolar . The former represented polar solvation energy for unbound protein, unbound ligand, and protein-ligand complex, while the latter was nonpolar solvation energy for them. The entropic contribution was not considered due to a huge amount of computation and low accuracy. ΔG int was taken as zero because the receptor-ligand complex was usually considered as an invariant conformation in calculation. Consequently, binding free energy (G bind ) was composed of Van der Waals' force (G vdw ), electrostatic interaction (G ele ), polar solvation energy (G polar ), and nonpolar solvation energy (G apolar ). Thus

Decomposition of MM/PBSA binding free energy
The binding free energy decomposition was performed to obtain the contribution of key residues for the sake of analyzing interactions between ligands and LXR through a dynamic way. g_mmpbsa tool in Gromacs 5.1.4 was used to extract a section of 2-ns trajectory from MD results and turned the trajectory into a total of 100 snapshots. Subsequently, decomposition of residues was calculated. Finally, a script written in Python was applied to extract detailed information about binding free energy decomposition in the form of ΔG MM , ΔG polar , and ΔG apolar .

Binding model analysis
The binding model analysis was performed to explore the key residues of LXR-ligand binding. The output files of molecular docking and pdb files of LXRβ were combined as the raw material of binding model analysis. Then, Ligplot + (version 2.2) [50] was introduced to analyze the binding pattern in 2D. 3D LXR-ligand interaction diagram was analyzed and illustrated by PyMOL (version 1.8).

Structural analysis of LXR and construction of validation standard
In order to explore the structure of LXRα and LXRβ and search residues with significant differences between them, alanine scanning and pairwise sequence alignment were used to evaluate the stability effect of residues after site-directed mutation and differences of protein sequences, respectively. We also selected selective and nonselective agonists from different databases and discovered differences of the binding mode between them. We regarded residues which bind to all selective agonists but bind to nonselective agonists poorly as LXRβ selective residues. The residues were used as a standard to verify the selectivity of screened compounds.

Comparison of protein sequences of LXRα and LXRβ
In order to analyze the differences of protein sequences between LXRα and LXRβ, pairwise sequence alignment of LXRα and LXRβ was executed by APIs tool using Smithwaterman algorithm. The pairwise sequence alignment results are illustrated in Fig. 1. The full length of LXR attained 62.5% sequence identity and 72.7% sequence similarity. Ligand binding domain (LBD) shared 88.5% sequence similarity and 78.3% sequence identity, which was consistent with the previous research [51]. In addition, the full length of LXR attained 9.2% gaps, and there were no gaps in LBD sequences, which means that the LBD possessed a higher degree of homology and gaps were focus on non-LBD regions.

Stability landscape of all possible mutations of LXR
In order to evaluate the stability effect of all possible mutations of different residues, alanine scanning on LXRα and LXRβ was executed by calculating the change of free energy of mutations. The results were shown in Fig. 2. Here, ΔStability between 0.5 kcal·mol −1 and −0.5 kcal·mol −1 means the effect of mutations was neutral; ΔStability between 0.5 kcal·mol −1 and 1 kcal·mol −1 means the effect of mutations was slightly destabilizing, while values of −0.5 kcal·mol −1 to −1 kcal·mol −1 means the effect of mutations was slightly stabilizing; ΔStability between 1 kcal· mol −1 and 2 kcal·mol −1 means the effect of mutations was destabilizing, while values of −1 kcal·mol −1 to −2 kcal· mol −1 means the effect of mutations was stabilizing. ΔStability over 2 kcal·mol −1 means the effect of mutations was highly destabilizing, while ΔStability below −2 kcal· mol −1 means the effect of mutations was highly stabilizing. In LXRα, 178 of the 208 residues became unstable, accounting for as much as 85.58%. Of the 217 mutational residues, 181 destabilized accounting for 83.41% in LXRβ. It suggested that most residues on LXRα and LXRβ attained high thermodynamic stability effects. As shown in Fig. 3, the color of different residues from red to blue means the values of ΔStability changed from low to high. Residues whose ΔStability values were below −2 kcal·mol −1 were marked as spheres. The ΔStability values that were less than −2 kcal· mol −1 were exhibited in Table S1. The Lys 273 , Asp 318, Lys 410 , Ser 418 , Lys 433 , and Pro 437 were shared by LXRα and LXRβ, while the residue numbers in LXRβ are Lys 287 , Asp 332 , Lys 424 , Ser 432 , Lys 447 , and Pro 451 . The ΔStability of the residues of LXRα was slightly higher than that of the LXRβ, where the average of ΔStability of the former was −5.60 kcal·mol −1 and the latter was −5.72 kcal·mol −1 . The Gly 212 , Asp 270 , Asp 285 , Asn 347 , Asp 374 , Asp 400 , Asp 444 , and His 446 were the specific residues to LXRα. The minimum of ΔStability of the residues was −8.50 kcal·mol −1 , and other residues also possessed considerably low values, which means that the residues possessed an unstable structure before mutating. The Lys 240 , Lys 248 , Gly 296 , Lys 331 , Lys 337 , Asp 366 , Asn 377 , Pro 398 , Pro 412 , and Ser 454 were specific mutational residues to LXRβ. The average of ΔStability of the residues was −3.29 kcal·mol −1 , and the minimum was only −5.98 kcal· mol −1 . In addition, most of the residues except Lys 248 were located on helix regions, which are agonist conformation found in some developed LXRβ agonists especially helix 3 (residues from Arg 261 to Val 289 ) and helix 12 (residues from Pro 451 to Trp 457 ) [52,53]. Coincidentally, Ser 454 was located on helix 12, indicating it attained favorable binding activity with some developed agonists such as T0901317 and GW3965. Based on this, Ser 454 may become a potential residue for screening LXRβ.

Exploring specific key residues binding with LXRβ
In order to analyze the differences underlying the binding mode of agonists between LXRα and LXRβ, 22 known selective or nonselective agonists were docked to LXR, and the results were shown in Table 1. The binding affinity could be calculated by docking score, and we proposed that the selectivity of agonists between LXR could be reflected by ΔScore, the docking score of agonists against LXRβ minus that of LXRα. The lower value of ΔScore indicated the higher selectivity against LXRβ. As can be seen in Table 1, BDBM27173, BDBM27174, and BDBM50300572 [54] showed favorable binding affinity against LXRβ and poor binding affinity against LXRα, and the ΔScore of them were -8.0 kcal·mol −1 , −6.1 kcal·mol −1 , and -3.8 kcal·mol −1 , respectively. These compounds are analogs of GW3965 [55], a tertiary-amine LXR agonist. GW3965 used to be considered as a nonselective and nonsteroidal agonist for LXR. However, it has been demonstrated that GW3965 was a potent full agonist that possessed LXRβ specificity (EC50 values are 176 nM and 15 nM, respectively) [56]. Therefore, GW3965  a LXRα and residues whose ΔStability is less than −2 kcal·mol −1 . b LXRβ and residues whose ΔStability is below −2 kcal·mol −1 was classified as a selective agonist. RGX-104 [57], a small molecule LXR agonist, obtained −5.6 kcal·mol −1 of ΔScore. WYE672 as a LXR agonist showed potent binding affinity to LXRβ (IC 50 = 53 nM). It had a little binding affinity for LXRα (IC 50 > 1.0 μM) [58], and its ΔScore was −3.6 kcal· mol −1 . CHEMBL4215604 is a 2-thienyl analog exhibiting potent LXRβ selective agonistic activity (EC 50 = 0.559 μM) [ 5 9 ] , a n d i t s Δ S c o r e w a s − 3 . 1 k c a l · m o l − 1 . N-Acylthiadiazoline [60] is a selective LXRβ agonist with poor pharmacokinetic property, and its ΔScore calculated here was − 2. 8 kc al· mo l − 1 . 4 -( 3-Bi ar yl) qu i n ol i ne s ulf on e (BDBM50317731) is a high-affinity LXRβ agonist with modest binding selectivity over LXRα, and its ΔScore calculated here was −2.3 kcal·mol −1 [10].
Nonselective agonists, such as AZ876 [61], LXR623 [62], and T0901317 [63], attained ΔScore higher than −2 kcal· mol −1 . It has been indicated that the experimental binding free energy ΔG exp approaches to -RT lnIC 50 ; thus for every difference in binding free energy of 1 kcal·mol −1 , the difference in activity is about 5.9 times. Thereby we suggested that a value of −2 kcal·mol −1 was the threshold of whether a compound possesses selectivity. An agonist would possess selectivity against LXRβ if the ΔScore is lower than −2 kcal·mol −1 and a compound with ΔScore over −2 kcal·mol −1 is more likely to be a nonselective agonist. The physiological ligands for LXR are oxysterols, including 24(S)-hydroxyl cholesterol (cpd.9) and 24(S),25-epoxycholesterol (cpd.8); they did not show observed selectivity between LXRα and LXRβ compared with other synthetic ligands. Meanwhile, the affinity against LXR of physiological ligands calculated here was also more unfavorable than synthetic ligands. However, despite many efforts that have been made to discover effective and selective LXRβ agonists, there is rarely an approved drug that has been reported.
Furthermore, we investigated the differences of the binding mode of selective and nonselective agonists targeting LXR using 2D and 3D ligand-receptor interaction diagram. Among 22 known agonists collected from literature, 9 agonists are selective agonists of LXRβ, and the rests are nonselective agonists. As shown in Tables 2, 29 key residues from LBD participated in ligand binding. Residues Leu 274 , Thr 316 , and Phe 329 participated in LXRβ binding for all selective and nonselective agonists, indicating these residues are very important for agonist binding. Residues Phe 271 , Ser 278 , Met 312 , and His 435 were involved in selective agonist binding, and they all participated in 9 selective agonists binding but partially occurred in the binding of nonselective agonists. It has been Compound with a unique database ID, which takes the form of a "BDBM," '"HEMBL," or "SCHEMBL" prefix immediately by an integer. a The Binding Database. b CHEMBL database. c SureChEMBL database The ΔScore of collected nonselective agonists was all over −2 kcal·mol −1 but not over 0 kcal·mol −1 , for results of LXRβ binding patterns with agonists that possess LXRα specificity would not be considered in this study indicated that the hydrophobic part of the agonist should be located in the hydrophobic pocket constituted by Phe 271 , Phe 340 , and Phe 349 to simulate the activity of LXR [64].
Here, we suggest that Phe 271 , Ser 278 , Met 312 , and His 435 might be very important for selective agonist binding with LXRβ. The binding pattern between agonists and LXRβ was predicted by Ligplot + , and the results of eight representative ligands were illustrated in Fig. 4 and Fig. 5. First of all, as mentioned above, Leu 274 , Thr 316 , and Phe 329 which were considered as important residues for agonists were found in all agonists. Glu 315 , a residue located on the helix region centered on Arg 319 , was the specific residue binding with selective agonists. It was worth noting that only selective agonists could bind with Trp 457 while Trp 457 was not a key residue in nonselective agonists-LXRβ binding. Trp 457 located on helix 12 which was an agonist conformation is found in some developed LXRβ agonists. Trp 457 had been also proved to form a close contact to His 435 to help it form polar interaction with ligands easier [52]. Therefore, a standard that we verify whether compounds possess LXRβ selectivity was constructed. Firstly, Trp 457 was selected to identify a compound whether possessed selectivity due to its unique properties. Then, His 435 was indispensable, for it was the bridge of interaction between Trp 457 and ligands. In addition, Phe 271 , Ser 278 , and Met 312 were involved in all selective agonist binding. Taken together, a compound which targets LXRβ selectively should

Identifying selective agonists targeting LXRβ
In order to evaluate pharmacokinetic properties and binding affinity of agonists from 112 terpenes of alismatis rhizoma, ADMET analysis and molecular docking were adopted to screen compounds. Binding free energy and MD simulations were introduced to further eliminate ligands whose binding affinity is unfavorable.

Pharmacokinetic and physicochemical analysis of collected terpenes
In order to understand the physiochemical and pharmacokinetic properties of 112 collected terpenes, the ADMET (absorption, distribution, metabolism, excretion, and toxicity) properties were predicted. RO5 which signifies the number beyond Lipinski's rule of five regions was calculated based on the result. The values of RO5 of all compounds were less than 3, and that of 97 terpenes were located on the recommended range. A total of 52 compounds showed an RO5 of 0, which means the compounds possessed potential druggability. Log S, log D, and log K p stand for aqueous solubility, distribution coefficient D at PH = 7.4, and skin permeability, respectively. All of them are important properties to evaluate the physicochemical properties of drugs. Caco-2 permeability was represented by log Papp. In our results, 97 of 112 terpenes showed that the values of log Papp were located on the recommended range. HIA means human intestinal absorption. Compounds whose HIA is more than or equal to 30% were marked as 1, while HIA less than 30% were marked as 0. Of the 112 terpenes, 69 showed that values of HIA were 1, which means that the body may have a high absorption for these compounds. BBB stands for blood-brain barrier, and too polar compounds will not cross through the blood-brain barrier [65]. 0 means that a compound cannot pass through the blood-brain barrier, while 1 represents a compound that can go through the blood-brain barrier. Results showed that none of the compounds can pass through the barrier. However, LXRβ was mainly distributed in livers, and therefore BBB has little effect on the treatment effect of the compounds we screened. DILI means drug-induced liver injury, and 0 represents that the drug will cause little damage to the liver function. Compounds we screened must make sure the values of DILI were 0, for the drugs will mainly affect the liver. CYP1A2, CYP2C19, CYP2C9, CYP2D6, and CYP3A4 inhibitors were the probability of a molecule being the inhibitor of cytochrome P450 families. The ADMET properties of the compounds involved in this study are presented in Table 3

Assessment of the reliability of molecular docking mode
In order to verify the reliability of the docking models, ROC curves were used to evaluate the performance of docking functions that distinguish active compounds and inactive compounds from a given dataset. A set of decoy molecules for ligands which were obtained from the ChEMBL database were constructed while they were input to dock with LXRβ with different score functions. Then, the ROC curves of each function were plotted in Fig. 6a, and the AUC values of LeDock, Grid, Amber, and Vina are 0.710, 0.729, 0.786, and 0.88,7 respectively. The range of AUC values is from 0 to 1. In addition, 0.5 is deemed as the threshold of whether a function has ability to distinguish [70]. AUC < 0.5 suggests that the test results are worse than random results while AUC > 0.5 shows that test is effective [71]. Based on this, the study chose an AUC value of 0.5 as a reference line, which was named Random in Fig. 6a. To compare the AUC values of different functions with a reference line, all of them demonstrated acceptable discrimination, and the Vina score function was the most outstanding of them. Then, redock with LXRβ was executed by AutoDock Vina for the sake of identifying the accuracy of the Vina function in the most straightforward way [72]. RMSD is thought to be an essential standard of evaluating how much structure is different from the original structure [73]. Furthermore, if the RMSD value between pre-and post-docking structure is below 2 Å, the docking model is considered to be great and can be applied for docking [74]. Herein, the RMSD value between redock structure and its initial structure of LXRβ using Vina function in the study was 1.9 Å, which met the requirement mentioned above. The result of redock and its initial structure is shown as Fig. 6b, and the binding pocket of LXRβ used in the study is displayed in Fig. 6c.

Screening potential LXRβ ligands using molecular docking
The Vina scores were applied to compute binding affinity between ligands and LXR using AutoDock Vina program. In order to screen ligands targeting LXRβ selectively, the ΔScore was calculated which means the values of Vina score of LXRβ minus that of LXRα. As mentioned above, we thought that compounds whose ΔScore values below −2 kcal·mol −1 possessed preferential binding affinity with LXRβ. Seven screened compound candidates which met the requirements were shown in Table 4. cpd.8 and cpd.9 which were thought to activate LXR endogenously in the liver and brain, respectively [75], were also presented as reference ligands in Table 4. Compared with cpd.8, all screened compound candidates showed higher binding affinities of LXRβ, and most of them except cpd.4 attained lower binding affinities of LXRα, suggesting they are more likely to bind with LXRβ. Moreover, the values of HIA of most the candidates except cpd.2 were 1, indicating most candidates The recommended range of properties referred from ADMETlab [31] and SwissADME [30]. a Average of octanol/water partition coefficient by different prediction, which consists of iLOGP [66], XLOGP3 [67], WLOGP [68], MLOGP, and SILICOS-IT. b Aqueous solubility, which was implemented from ESOL [69]. c Distribution coefficient D at PH = 7.4. d Caco-2 permeability. e 0, non-substrate; 1, substrate. f Human intestinal absorption. g Skin permeability, unit is cm/s. h Blood-brain barrier. i Plasma protein binding, significant with drugs that are highly protein-bound and have a low therapeutic index. j 0, non-inhibitor; 1, inhibitor. k hERG (human ether-a-go-go-related gene). Blockers; 0, non-blockers; 1, blockers. l Druginduced liver injury; 0, DILI negative; 1, DILI positive The results of redock, carbon skeleton shown as white is the original structure while carbon skeleton shown as green is the redocking structure. c The binding pocket of LXRβ possessed the potential to enter into the internal environment and activate target protein. Among all candidates, cpd.5 was the most potential ligand for its performance in pharmacokinetic analysis and a great gap in docking affinities with cpd.8. cpd.1 and cpd.4 were also considerable compounds for LXRβ agonists. The former attained a high value of ΔScore, which means that it highly selectively targets LXRβ. The latter possessed a high docking score of LXRβ. Most of all, the values of DILI of all screened compound candidates were 0, which means they all possessed essential properties of drugs targeting LXRβ that is low toxicity. The discussion about cpd.9 will not proceed because it was a ligand found in the brain.

RMSD analysis of MD simulations
MD simulations at 20 ns for LXR-ligand were performed via the Gromacs program in order to evaluate dynamic interactions between LXRβ and ligand candidates. The calculated results of MD simulations were output in the form of an MD trajectory. Then, the assessment of backbone RMSD values for each MD trajectory was executed to provide a holistic sight to observe the stability of the LXR-ligand system. As illustrated in Fig. 7, the backbone values of the LXRα-ligand system were between 0.10 and 0.30 nm while that of the LXRβligand system were from 0.075 to 0.275 nm. The reference ligand, cpd.8, rapidly reaches the platform at 2 ns in both the LXRα-ligand system and LXRβ-ligand system. Cpd.1 and cpd.6 showed drastic fluctuation of around 0.20 nm and 0.225 nm in LXRα-ligand system, and LXRα-cpd.7 system appeared stable at a minimum equilibrium of 0.18 nm. The LXRβ-cpd.6 system fluctuated steadily while LXRβ-cpd.1 system showed fluctuation as drastic as the LXRα-ligand system. LXRβ-cpd.2 system was the most stable LXRβ-ligand system, which fluctuated at 0.125 nm. Moreover, fluctuations of LXRα-cpd.2 system were on the mean level among all LXRα-ligand systems. LXRβ-cpd.4 was the most fluctuant LXRβ-ligand system, which kept fluctuating during the 20-ns simulation experiment. However, LXRα-cpd.4 system showed a stable equilibrium. cpd.3 and cpd.5 presented increasingly large unstable fluctuation of RMSD values in the LXRβ-ligand system whereas the level of RMSD was moderately low in the LXRα-ligand system. In order to screen ligands selectively targeting LXRβ, compounds which fluctuated steadily in the LXRβ-ligand system but waved drastically in the LXRα-ligand system were selected. It demonstrated that cpd.6 was the molecule that most fits the requirement and cpd.1 and cpd.2 can also be candidates.

Calculation of MM/PBSA binding free energy
Because Vina score only computes the binding free energy of the lowest-scoring conformation and exist empirical component, MM/PBSA binding energy calculation was introduced to the work in order to assess binding affinity accurately between ligands and proteins [36,49]. The results of calculation of LXRα and LXRβ were listed in Table 5 and Table 6, respectively. Here, ΔG bind means a comprehensive evaluation of binding free energy while ΔG vdw , ΔG ele , ΔG polar , and ΔG apolar are key elements of binding free energy. The lower binding free energy suggests the binding mode between ligand and receptor is more stable. As shown in Table 5, cpd.2 and reference ligand cpd.8 presented a similar level in ΔG bind , which were − 45.84 kcal·mol −1 and −45.46 kcal· mol −1 , respectively. ΔG bind of cpd.4 was −40.67 kcal·mol −1 , which means that interactions between cpd.4 and LXRα are the most unstable. cpd.7 was the most stable ligand that binds to LXRα. The rests were at the same and moderately stable level. As listed in Table 6, the reference compound, cpd.8, showed considerably high ΔG bind of −30.29 kcal·mol −1 , which indicated that ΔG bind was detrimental to binding to LXRβ. ΔG bind of cpd.6 was −56.17 kcal·mol −1 , which was the lowest value among all screened compound candidates, indicating cpd.6 possessed the highest binding affinity to LXRβ. The binding affinity with LXRβ of cpd.1 and cpd.7

RMSF analysis of the key amino acid residues
RMSF values for each MD trajectory were calculated in order to assess the flexibility of proteins while higher values indicated better flexibility. Figure 8 shows the results of RMSF values of the LXR-ligand system. Low RMSF values mean that the system is stable, and therefore residues of this region  may participate in the binding interactions of the LXR-ligand system. It is clear that all compounds shared similar RMSF profiles in both the LXRα-ligand system and LXRβ-ligand system. Compared with the LXRβ-ligand system, the LXRαligand system showed relatively higher RMSF values of average 0.1328 nm while that of the LXRβ-ligand system was at a mean level of 0.1145 nm, which suggested screened compounds are more likely to bind to LXRβ on the whole. In the LXRβ-ligand system, region Ser 254 -Leu 264 was the most flexible, indicating the residues contribute less to the binding of the system. RMSF values of region Ala 306 -Arg 319 and Gly 369 -Ile 377 were lowest, suggesting the region is more responsible for the binding of the system. Interestingly, the unique mutational residue Asn 377 in LXRβ was located in the region, which means that Asn 377 contributes a lot to the binding of LXRβ and it is a candidate screening target for LXRβ. Region Ala 306 -Arg 319 was all situated on helix region except Ser 307 , and residues of protein-binding pocket Met 312 , Leu 313 , Glu 315 , Thr 316 , and Arg 319 also lay on the region. Met 312 which has been proved to bind to LXRβ such as T-0901317 [52] is a part of the validation standard for selective targeting LXRβ binding residues.
The RMSF values of all binding system on Met 312 were at a low level, which indicates that binding between screened compounds and Met 312 was stable. Arg 319 was a vital residue which participated in forming hydrogen bonds with ligands. Thr 316 is thought to be an active site of forming hydrogen bonds with LXRβ [76]. cpd.1, cpd.2, and cpd.6 all possessed the ability to bind with Thr 316 on LXRβ in the form of hydrogen bonds.

Decomposition of MM/PBSA binding free energy
Binding free energy decomposition was applied to identify the contributions of each residue in the LXR-ligand binding model. The contributions of each residue of LXRα and LXRβ are depicted in Fig. 9 and Fig. 10, respectively. Low values of energy indicate an important role of residues in the LXRligand binding. As demonstrated in Fig. 9, most residues showed values which were negligible to the system. Only the LXRα-cpd.8 system had a moderately high positive value located on Arg 305 , which means that the residue contributes less to binding to LXRα. Phe 257 , Leu 260 , Ala 261 , Met 298 , Leu 331 , and Trp 443 (values of energy of selected residues are below −0.8 kcal·mol −1 ) showed relatively larger contributions of binding to LXRα. As illustrated in Fig. 10, most residues had little effect on LXRβ-ligand binding. The reference, LXRβ-cpd.8 system, showed two extremely high values of binding free energy of 8.21 kcal·mol −1 and 3.71 kcal·mol −1 , which were located on Glu 281 and Glu 315 , respectively. Energy of Arg 319 was high in LXRβ-cpd.1 system while the Phe 271 , Phe 329 , and Trp 457 showed the lowest values of binding free energy. Phe 271 and Phe 329 have been proved an important contributor to hydrophobic interactions between LXRβ and ligand, and Trp 457 can form a close contact to His 435 which is able to form polar interaction with ligand [52]. All of the four residues are within validation standard and indicated cpd.1 possessed LXRβ selectivity. Arg 319 is an inhibitive site which can form hydrogen bonds with ligands. Compared to the residues above, a smaller change in ΔG MM is a direct cause of the higher energy value of Arg 319 . ΔG bind of Phe 271 , Phe 329 , and Trp 457 were lower than the rest residues because of lower ΔG MM while ΔG polar were even higher than the rest residues. ΔG MM indicates potential energy in the vacuum which are composed of ΔG vdw and ΔG ele . It demonstrated that potential energy was a crucial factor that affects the binding between cpd.1 and key residues in LXRβ. Energy of Glu 281 and Phe 329 was at a maximum of 2.22 kcal·mol −1 and a minimum of −2.80 kcal·mol −1 of LXRβ-cpd.2 system, respectively. Glu 281 , whose energy was highest in both LXRβ-cpd.2 and LXRβ-cpd.8 system, showed an extraordinary ΔG polar of 2.48 kcal·mol −1 and 8.80 kcal·mol −1 , respectively. Glu 281 was considered as a hydrophilic site which is able to bind to the ligand in the form of hydrogen bonds [77]. Hydrogen bonds are usually considered as a set of dispersion force and intermolecular or interatomic electrostatic interaction [78], whose energy is represented by ΔG MM . However, ΔG MM of the two compounds were not low, and therefore

Binding model analysis of screened compounds
In order to investigate the differences of 2D and 3D LXRligand binding mode of screened compounds, LigPlot + [50] was introduced to identify key residues which interact with ligands in LXRβ. Key residues of all screened compound candidates were predicted, and they were shown in Table S2. The results showed that most of them attained plentiful binding residues at helix3 and helix12, indicating they possessed a high binding affinity with LXRβ. The results of screened compounds and reference ligand was plotted in Fig. 11. Circled residues were equivalent side chains, which were considered as key residues in binding free energy decomposition. The green line represented hydrogen bonds between ligands and residues, and residues that formed hydrogen bonds were highlighted particularly. As depicted in Fig. 11, Phe 271 , Thr 316 , Ser 278 , and His 435 formed hydrogen bonds with screened compounds. These residues were also considered very important for selective agonist binding with LXRβ in our study. Interestingly, Thr 316 , a residue found in both selective and nonselective agonists, presented in all screened compounds, indicating all the compounds may be potential agonists against LXRβ. His 435 was proved to be an important residue that forms hydrogen bonds to stabilize the protein.
A binding pocket [52] could be found in all screened ligands. However, the reference ligand cannot bind with His 435 . Another residue, Trp 457 , which was located on helix 12 and could form contact with His 435 [52], had the same case. Arg 319 and Glu 315 were located in protein binding pocket, and Arg 319 Fig. 10 Binding free energy decomposition of LXRβ-ligand model (error bar of residues with an absolute value of energy less than 0.1 kal·mol −1 was not shown). a LXRβ-cpd.8 system. b LXRβ-cpd.1 system. c LXRβcpd.2 system. d LXRβ-cpd.6 system was an inhibitive site which can form hydrogen bonds. However, results of energy decomposition demonstrated that the binding affinity between ligands and Arg 319 on LXRβ was not high. And the binding affinity between ligands and LXRβ of Glu 315 was average.
The interactions between LXRβ and ligands in 3D were illustrated in Fig. 12 to show the hydrogen bonds intuitively. More hydrogen bonds and shorter length of hydrogen bonds indicate larger interaction intensity between ligands and receptors. It is clear that the intensity of interaction between LXRβ and screened compounds was higher than that of reference ligand, which was reflected on the number and length of hydrogen bonds between LXRβ and ligands. All screened compounds could bind to Thr 316 , and the distance between ligands and the residues was the nearest among all residues in Fig. 12, indicating ligands possessed a high binding affinity with Thr 316 . Phe 271 was the residue which binds closest with cpd.1, suggesting the residue was important in LXRβ-cpd.1 binding. Both His 435 and Ser 278 which were located on protein binding pocket had an average binding affinity, which means that all screened ligands could bind with the two residues. However, reference ligand, which was a kind of nonselective agonists in our body, did not find the characters described above.

Discussion
Previous studies showed that activation of LXR activity could reduce TC, which is responsible for hyperlipidemia [1,8]. Many LXRβ agonists have already been screened such as T0901317 and LXR623 [62,63]. However, the agonists are not selective, which can also activate LXRα. LXRα played an essential role in activating SREBP-1c protein, which was considered as a pivotal activator of lipogenic genes [9]. However, activation of LXRβ had less impact on lipogenic synthesis pathway [14]. Therefore, screening a compound selectively targeting LXRβ becomes urgent. Alismatis rhizoma, a kind of TCM to treat dysuria, edema, obesity, diabetes, hyperlipidemia, and hypertension, was introduced to be the raw material of the study [18]. Terpenes from alismatis rhizoma were collected, for they have been considered as active components to enhance transactivation for FXR which belongs to nuclear receptors superfamily along with LXR. Especially protostane-type triterpenes had been proved the major contributor to pharmacodynamic and bioactive compounds in alismatis rhizoma [21]. Then a series of virtual experiment was performed.
First of all, in order to understand the structural differences between LXRα and LXRβ and specific residues of binding of selective agonists to LXRβ, a structural analysis of LXRα and LXRβ was conducted, and a validation standard was constructed. Pairwise sequence alignment of full length and LBD region was performed to identify regions of similarity of LXRα and LXRβ. The results showed that LXRα and LXRβ shared the same length of LBD region, and the gaps of LBD were 0. The %similarity and %identity of LBD were higher than that of the full length, suggesting LBD possessed a higher degree of homology. Alanine scanning was applied to evaluate the stability effect of different residues in site-directed mutation. The ΔStability of LXRα and LXRβ was similar in general, with only some local differences. Residues whose ΔStability below −2 kcal·mol −1 were selected, and LXRα was more stable than LXRβ in mutations for ΔStability of the residues of LXRα that was slightly higher than that of LXRβ in homologous residues. Moreover, the Lys 240 , Lys 248 , Gly 296 , Lys 331 , Lys 337 , Asp 366 , Asn 377 , Pro 398 , Pro 412 , and Ser 454 were unique mutational residues to LXRβ. All of them was located on the helix region except Lys 248 . Subsequently, 22 known selective and nonselective agonists were selected in order to analyze key residues that affected binding mode of LXRβ and ligands. The ΔScore means that the docking score of agonists against LXRβ minus that of LXRα and whether values below −2 kcal·mol −1 was threshold that a compound possessed selectivity. Then, LigPlot + was introduced to predict potential residues that interacted with ligands. The results showed that Phe 271 , Ser 278 , Met 312 , and His 435 might be very important for selective binding of ligands due to their high engagement in selective agonists binding. Moreover, a comparation between four representative selective agonists and four nonselective agonists was executed. The results showed that Trp 457 was responsible for selective binding between ligands and LXRβ. Finally, a standard wherein we verify whether compounds possessed LXRβ selectivity was set up; that is, ligands which bind with Phe 271 , Ser 278 , Met 312 , His 435 , and Trp 457 in LXRβ may be potential selective agonists targeting LXRβ.
Then, a series means of computer-aided drug design were executed to evaluate agonists against LXRβ from 112 terpenes of alismatis rhizoma. Suitable ADMET properties were the first hurdle for most ligand candidates [29], and unfavorable ADMET properties would lead that drugs possess low bioavailability. RO5, DILI, and HIA were used to assess whether properties of ADMET are favorable. The results showed that most of compounds possessed favorable ADMET properties. The ADMET properties of screened compound candidates in the study were all favorable except the absorption of cpd.2 was not high. Subsequently, molecular docking was employed to screen ligand whose binding affinity of LXRβ was far higher than that of LXRα. High values of AUC of ROC curves of docking functions and high RMSD value between redock structure and its initial structure of LXRβ suggested that the virtual screening model was reliable for screening active ligands effectively from compound database. The results of molecular docking showed that most terpenes have favorable binding affinity with LXRβ. The ΔScore value of −2 kcal·mol −1 which was considered as the threshold of selective agonists and nonselective agonists is used to screen compounds targeting LXRβ, and seven compound candidates were labeled as cpd.1 to 7. In order to further explore and verify the result of molecular docking, RMSD value analysis of MD simulations and calculation of MM/PBSA binding energy were followed. The results of MD simulation were used to evaluate dynamic interactions between screened compound candidates and LXR. Compounds whose RMSD curve of MD trajectory was stable were considered to possess high binding affinity with LXR. In this study, molecules with a stable trajectory of LXRβ and fluctuating trajectory of LXRα were thought to possess high dynamic binding affinities to LXRβ. It demonstrated that cpd.6 was the candidate that most fits the requirement and dynamic binding affinity of cpd.1 and cpd.2 was also favorable. cpd.4 was dramatically fluctuant during the specified simulative time of the LXRβ-ligand system, which means that the molecule is averse to binding to LXRβ. An MD trajectory of 2 ns was extracted to calculate the MM/PBSA binding free energy. Compared to the result of docking, the result of binding free energy possesses distinct distribution patterns. Some compounds of high ΔScore values possessed binding free energy of LXRα higher than that of the LXRβ. Because the Vina score only computes the binding free energy of the lowest-scoring conformation and existing empirical component, we adopted computation of g_mmpbsa which combines MM energy of the gas phase and solvation free energy as the results of binding free energy to simulate ΔG bind in vivo [36,49]. Screened compound candidates whose ΔG bind of LXRβ were higher than that of the LXRα were eliminated, and cpd.1, cpd.2, cpd.4, and cpd.6 have remained. Taken together, we screened three compounds through a series of experiments, namely cpd.1, cpd.2, and cpd.6. The HIA properties of cpd.2 were unfavorable, and cpd.4 did not meet the standard for it did not pass the screening of MD simulations.
Later on, RMSF analysis of MD simulations, binding free energy decomposition, and binding model analysis were carried out to find the differences underlying the binding mode of agonists between LXRα and LXRβ. The standard we proposed earlier was used to verify whether screened compounds possessed LXRβ selectivity. The RMSF values of molecular dynamic simulations provided flexibility of proteins. RMSF values of region Ala 306 -Arg 319 and Gly 369 -Ile 377 were lowest, suggesting the regions are more responsible for binding between LXRβ and the screened compounds. It is worth noting that Met 312 , Thr 316 , and Asn 377 were on the region. Thr 316 is an active site of forming hydrogen bonds with LXRβ [76], and it can be found in both selective and nonselective agonists. Asn 377 is specific residues of LXRβ in the structural analysis of alanine mutations. It demonstrated that screened compounds might possess the ability of agonists and bind active residues on LXRβ. Met 312 was a hydrophobic site which interacts with ligand by hydrogen bonds [52], and it was one of residues of verification standard. RMSF values of Met 312 were one of the lowest on the entire residues, suggesting the three compounds might possess LXRβ selectivity. Decomposition of MM/PBSA binding energy was also executed to calculate contributions of each residue in the LXR-ligand binding model. The results showed that the residues of lowest binding free energy of LXRβ-cpd.1 system were Phe 271 , Phe 329 , and Trp 457 . The residues of the lowest binding free energy of LXRβ-cpd.2 system were Phe 329 , and those of LXRβ-cpd.6 system were Phe 271 , Leu 274 , and Trp 457 . Leu 274 and Phe 329 participated in the LXRβ binding for all selective and nonselective agonists, indicating all of screened compounds possessed potential to be agonists. Phe 271 has been proven to be an important contributor to hydrophobic interactions between LXRβ and ligand, and Trp 457 can form a close contact to His 435 which is able to form polar interaction with the ligand [52]. They are also the members of verification standard. As mentioned earlier, cpd.1 and cpd.6 could bind to all of them, indicating the two compounds might possess LXRβ selectivity. The binding mode between LXRβ and screened compounds were calculated by LigPlot + and PyMOL. The results showed that cpd.1 could form hydrogen bonds to Phe 271 , Ser 278 , Phe 329 , and His 435 . cpd.2 could form hydrogen bonds to Thr 316 and His 435 . cpd.8 could form hydrogen bonds to Phe 271 , Ser 278 , and Thr 316 . Ser 278 was one of the residues that we thought contribute more to selective agonist binding with LXRβ. His 435 , which was the center of a protein binding pocket and could form hydrogen bonds to agonists [52]

Conclusions
In summary, the present study provides an alternative strategy for screening selective agonists targeting LXRα and LXRβ from terpenes. First of all, a standard wherein we verify agonists whether possess LXRβ selectivity was constructed by analyzing differences of residues between developed selective and nonselective agonists. Then, a series of computational tools such as ADMET analysis, molecular docking, calculation of binding free energy, and MD simulations were applied to evaluate physicochemical and biological properties of 112 compounds. Moreover, key residues and the binding mode of screened compounds were analyzed to confirm whether they possessed selectivity. Our result indicated that 16-hydroxyalisol B 23-acetate and Alisol M 23-acetate were identified as the most potential agonists. Phe 271 , Ser 278 , Met 312 , His 435 , and Trp 457 were considered as the major contributor to the binding between LXRβ and screened compounds. Finally, favorable physicochemical properties of screened compounds and remarkable binding affinity to residues which are in favor of the binding in previous study indicate that the result of the study is reliable. Thus, the screening strategies used in the study may guide the development of LXRβ-agonists from TCM.
Code availability All codes used in the study are mentioned in this article or can be obtained from the official website of the software.
Authors' contributions Chuanfang Wu and Jinku Bao conceived and designed the experiments; Chuanjiong Lin performed the experiment, and Jianzong Li analyzed the data; and Chuanjiong Lin wrote this manuscript. All authors read and approved the final manuscript. Data availability All data and material are included in this article.

Declarations
Ethical approval N/A Consent to participate N/A

Consent for publication N/A
Conflict of interest The authors declare that they have no conflict of interest.