Scaffold-Based Pan-Agonist Design for the PPARα, PPARβ and PPARγ Receptors

As important members of nuclear receptor superfamily, Peroxisome proliferator-activated receptors (PPAR) play essential roles in regulating cellular differentiation, development, metabolism, and tumorigenesis of higher organisms. The PPAR receptors have 3 identified subtypes: PPARα, PPARβ and PPARγ, all of which have been treated as attractive targets for developing drugs to treat type 2 diabetes. Due to the undesirable side-effects, many PPAR agonists including PPARα/γ and PPARβ/γ dual agonists are stopped by US FDA in the clinical trials. An alternative strategy is to design novel pan-agonist that can simultaneously activate PPARα, PPARβ and PPARγ. Under such an idea, in the current study we adopted the core hopping algorithm and glide docking procedure to generate 7 novel compounds based on a typical PPAR pan-agonist LY465608. It was observed by the docking procedures and molecular dynamics simulations that the compounds generated by the core hopping and glide docking not only possessed the similar functions as the original LY465608 compound to activate PPARα, PPARβ and PPARγ receptors, but also had more favorable conformation for binding to the PPAR receptors. The additional absorption, distribution, metabolism and excretion (ADME) predictions showed that the 7 compounds (especially Cpd#1) hold high potential to be novel lead compounds for the PPAR pan-agonist. Our findings can provide a new strategy or useful insights for designing the effective pan-agonists against the type 2 diabetes.


Introduction
Type 2 diabetes are characterized by hyperglycemia, insulin resistance and defects in insulin secretion, thus the patient with type 2 diabetes often suffer from symptoms of dyslipidemia, hypertension as well as obesity [1]. As reported, there are nearly 151 million individuals in the world affected type 2 diabetes, which will increase to 346 million at the end of 2012 [2]. As important members of nuclear receptor superfamily, Peroxisome proliferatoractivated receptors (PPAR) play essential roles in regulating cellular differentiation, development, metabolism, and tumorigenesis of higher organisms [3], thus these receptors have been considered as attractive targets for treating type 2 diabetes. PPAR receptors have 3 identified subtypes: PPARa, PPARb and PPARc [4]. Although the ligand-binding domains of the subtypes share 60%-70% sequence similarities, they have specific organization distributions and active functions. PPARa is expressed at high level in the brown adipose tissue, liver, kidney, heart, and skeletal muscle [5]. This PPAR subtype has significant functions in fatty acid metabolism and energy homeostasis [6] as well as modification of the high-density lipoprotein (HDL) circulations [7]. The activation of PPARa receptor is propitious to improve glucose tolerance so as to decrease the risks of developing atherosclerotic lesions [8]. Due to its crucial role in the fatty acid oxidation in the liver cells, PPARa can increase the energy consumption so as to achieve the controlling of energy metabolism and body weight [9]. Thus, this PPAR receptor can be used as an anti-obesity target, against which a number of clinical drugs (i.e., beclofibrate, bezafibrate, ciprofibrate, clofibrate, fenofibrate and gemfibrozil) have been approved by US Food and Drug Administration (FDA) to treat the dyslipidemia.
PPARb is widely distributed in various body tissues with responsibility for controlling blood lipid concentration and insulin sensitivity [10]. Experimental evidences show that this PPAR receptor functions as a regulator in fatty acid catabolism, energy balance and cholesterol reversion transportation [11]. In animal models of type 2 diabetes, this receptor is found to have the ability of improving insulin resistance and decreasing plasma glucose [12,13]. PPARc is expressed in adipose tissue, macrophages and vascular smooth muscles at high levels. The activation of this PPAR receptor can increase the differentiation of the fat cells, improve the storage of fatty acids and enhance insulin sensitivity in the adipose tissue, skeletal muscle and liver [6,14].
As attractive targets for type 2 diabetes, many PPAR agonists, including PPARa/c and PPARb/c dual agonists, have been designed and synthesized. Due to the undesirable side-effects, many PPAR agonists (i.e., muraglitazar, tesaglitazar, ragaglitazar, TAK559 and KRP297) are stopped by US FDA in the clinical trials [15,16,17,18]. An alternative strategy is to design novel agonist that can simultaneously activate PPARa, PPARb and PPARc with an aim to treat both insulin resistance and dyslipidemia [19]. In the current study, we employed molecular modeling technologies with a core hopping approach to screen the fragment database, with an aim of searching for novel PPAR panagonists to treat type 2 diabetes.

Initial Structures for PPARa, PPARb and PPARc Receptors
The initial structures of PPARa, PPARb and PPARc receptors were obtained from the Protein Data Bank (PDB) with PDB entries of 1k7l.pdb [20], 1gwx.pdb [6] and 1k74.pdb [20], respectively. The crystal structure of PPARb was released in 2000 with a resolution of 2.5 Å by Xu et al [6], who released the crystal structures of PPARa and PPARc one year later [20], with resolutions of 2.5 Å and 2.3 Å , respectively. Except for the polar hydrogen and heavy atoms, all the other atoms including nonpolar hydrogen atoms in these crystal structures were removed. Hydrogen atoms were subsequently added to the receptor structures of all the PPAR receptors based on the computational pKa values for each residue in the crystal structures. Finally the structural models for PPARa, PPARb and PPARc receptors were obtained after a steepest descent energy minimization with OPLS force field parameters [21].

Molecular Docking Procedure with a Core Hopping Approach
Molecular docking technology has been increasingly used in the course of drug research and development [22,23,24]. To improve the activity of a lead compound, we usually vary the side chains attached to a core part of the compound, because in many cases it is the side chains that bind to the receptors [25]. However, it also makes sense to vary the core structure to find novel compounds (scaffolds). In the present study, core hopping algorithm (Combi-Glid 2.5. Schrodinger LLC, New York, 2009) was adopted during the molecular docking procedure with the functions to perform both fragment-based replacing and docking. The core hopping strategy is to screen multiple scaffolds against a guiding structure, searching for alignments of potential attachment points on the scaffold with the attachment points on the guiding structure. Generally speaking, the core hopping process contains 4 major steps: The first step is to define some possible points at which the cores are attached to the scaffold. It is also the define combinations step from the Combinatorial Screening panel packaged in Schrodinger (www.schrodinger.com). The second step is to define the receptor grid file in the Receptor Preparation panel in Schrodinger. The third step is the core preparation with the Protocore Preparation module in Schrodinger to find the core attaching to the scaffold using the fragment database derived from ZINC database (http:// zinc.docking.org). The last step is to align and dock the entire molecular structures built up by the core and scaffold into PPAR receptors. The cores are sorted and filtered by goodness of alignment and then redocked into the receptor after attaching the scaffold, followed by using the docking scores to sort the final molecules.

Molecular Dynamics Simulation
The candidates obtained from the core hopping docking procedure were then subjected to a series of molecular dynamic simulations by the open software GROMACS 4.0 [26]. The PPAR receptors were parameterized by GROMOS 96-53a6 force field parameters [27], while the topology file, partial charges and force field parameters of the candidates were generated by the online tool PRODRG [28]. The simulation systems were subsequently inserted into 108 DPPC lipid bilayers, and solvated in a specific box with SPC water and a space of 9 Å around the solute. To neutralize the redundant charges of the simulation systems, 4 sodium ions were added into the simulation systems to randomly replace 4 water molecules. The neutralized systems were subjected to an energy minimization for about 3000 steps with the steepest descents approach. Finally, 10-ns molecular dynamic simulations were performed with constant temperature (300K), periodic boundary conditions and NVT ensembles.
During the molecular dynamics simulations, all bonds in the simulation systems were constrained by the Linear Constraint Solver (LINCS) algorithm [29,30], and the atom velocities for the start-up runs were obtained based on the Maxwell distribution at 300 K [31,32]. To maintain the simulation systems at a constant temperature and volume, the Berendsen thermostat with a coupling time of 0.1 ps and v-rescale scheme were applied [33]. The electrostatic interactions were treated by the particle mesh Ewald (PME) algorithm with interpolation order of 4 and a grid spacing of 0.12 nm [34,35]. The van der Waals interactions were treated by using a cutoff of 12 Å [36,37]. The integration step was set to 2 fs, and the coordinates were saved every 1 ps.

ADME Prediction
In the current study, QikProp was employed for predicting the absorption, distribution, metabolism and excretion (ADME) properties for all the candidates. QikProp is designed by Prof. William L. Jorgensen, and can predict physically significant descriptors (i.e., the partition coefficient, van der Waals surface area of polar nitrogen and oxygen atoms, and predicted aqueous solubility) and pharmaceutically relevant properties for small druglike molecules. This software also provides ranges for comparing a particular properties of molecules with those of 95% known drugs, and flags 30 types of reactive functional groups that may cause false positives in high-throughput screening assays.

ligand-binding Domains of the PPAR Receptors
Experimental evidences showed that the biological functions of PPAR receptors is regulated by the precise shape of their ligandbinding domain induced by the binding of ligands including a number of coactivators and corepressor proteins. These ligands can significantly simulate or inhibit PPAR receptor functions. According to the X-ray crystallography studies, the ligand-binding domains of PPARa, PPARb and PPARc receptors are very similar with the RMS deviations between Ca atoms less than 1 Å . Besides, they also share some common features ( Figure 1): i) composed of 12 a-helices arranged in an antiparallel helix sandwich, and a 4-stranded antiparallel b sheet; ii) Y-shaped hydrophobic ligand binding pocket with a volume of , 1300 cubic angstroms; and iii) a C-terminal helix (Helix 12 or AF2 helix) showing widely conformational variations in different crystals and playing essential roles in activation of PPAR receptors.
Addition to the common features mentioned above, the ligandbinding domains for PPARa, PPARb and PPARc receptors also have their unique features. For instance, it is significantly narrower in the region adjacent to the AF2 helix in PPARb, which is not suitable for the ligands containing larger polar heads. Compared with the PPARb ligand-binding domain, the one of PPARa is comparatively more lipophilic, while the one of PPARc is more hydrophilic. The largest variation among the three ligand-binding domains is found in the omega loop displaying a comparatively high RMS deviation and differing among a number of reported crystal structures.

Core Hopping and Drug Design
During the core hoping processes, LY465608 is selected as a guiding structure. This small molecule is firstly reported as a PPARa/c dual agonist, latter proved to excite PPARb as well. Owe to playing a positive role in keeping lipid and cholesterol homeostasis stability by reducing serum triglycerides and increasing HDL cholesterol, LY465608 has been treated as a typical PPAR pan-agonist to lower plasma glucose levels and improve insulin sensitivity. The structure of LY465608 can be divided into 3 major components ( Figure 2): a polar acidic head (Core A), a linker group (Core B) and a hydrophobic tial (Core C). The Core A contains an acidic head, and is reported to be of most importance for the ligand-binding and PPAR receptor activation. Thus, this part will be retained during the entire core hopping processes as described below.
The 1 st core hoppling operation was aimed at the Core B (see the red part of Figure 2) to generate 5 scaffolds, named Fragment B 1 to B 5 respectively. The 2 nd core hopping operation was aimed at the Core C (see the blue part of Figure 2) to generate 4 scaffolds, named Fragment C 1 to C 4 , respectively. Consequently, a total of 20 combinations of LY465608 derivatives were thus obtained. Subsequently, each of the 20 derivatives was docked into PPARa (1k7l.pdb), PPARb (1gwx.pdb) and PPARc receptors (1k74.pdb), respectively. The 20 derivative compounds were then ranked roughly according to their docking scores to PPARa, PPARb and PPARc receptors, respectively. The derivative compounds that employed stronger binding affinities than the original LY465608 with all the PPAR receptors were listed in Table 1 (the chemical structures of the derivative compounds mentioned above were included in Table S1). Of the top 7 derivative compounds, the Cpd#1 with Core A-Fragment B 1 -  Diagrammatically showing the core hopping procedure. The Original LY465608 structure contains 3 major components: 3 major components: a polar acidic head (Core A), a linker group (Core B) and a hydrophobic tial (Core C). Owe to forming significant hydrogen bonds with the ligand-binding domain, the Core A would be retained during the core hopping procedure. The 1 st core hopping operation was aimed at the Core B to generate 5 scaffolds named Fragment B 1 to B 5 respectively. The 2 nd core hopping operation was aimed at the Core C to generate 4 scaffolds, named Fragment C 1 to C 4 , respectively. Thus, a total of 20 combinations were obtained. doi:10.1371/journal.pone.0048453.g002 Fragment C 1 has the strongest binding affinity with all the PPAR receptors, and hence it was singled out for further studies.
As shown in Figure 3, the favorable binding mode of the Cpd#1 was aligned with the guiding structure LY465608 for PPARa, PPARb and PPARc receptors. According to the crystal studies, a conversed hydrogen bonding network composed of 4 significant hydrogen bonds formed by the acidic head of both LY465608 and Cpd#1 to the active site residues (Ser280, Tyr314, Tyr464 and His440 in PPARa receptor; Ser289, His323, Tyr473 and His449 in PPARc receptor) of PPARa and PPARc receptors were observed in our docking results. However, in PPARb receptor 3 hydrogen bonds were detected between the acidic head of both LY465608 and Cpd#1 and the active site residues (His323, His449and Tyr473). Based on the previous experimental and theoretical studies, these hydrogen bonds were believed to play crucial roles in stabilizing the AF2 helix in the active conformation, which is essential for the ligand-binding and receptor activation.
Addition to the conserved hydrogen bonding network mentioned above, Cpd#1 formed an additional hydrogen bond with Thr288 in PPARb receptor, which could notably enhance the binding affinity compared to the guiding molecule LY465608. The hydrophobic tail of the Core C in both LY465608 and Cpd#1 were buried well in the hydrophobic pocket of the ligand-binding domains of PPARa, PPARb and PPARc receptors. However, due to the suitable size, Cpd#1 was more fitted to the hydrophobic arm II, resulting in the stronger binding affinities than LY465608 (also see Table 1).

Dynamics Behaviors of the Cpd#1-PPAR Receptors
Molecular dynamics can provide useful information for characterizing the internal motions of biomacromolecules with time [38,39,40]. To study the dynamics behaviors of the Cpd#1-PPAR receptors, 10-ns molecular dynamics simulations were performed on apo, LY465608-bound and Cpd#1-bound states of PPARa, PPARb and PPARc receptors. The root mean square (RMS) deviation from the initial structure is considered as an important criterion usually used to measure the convergence of the protein systems concerned. In the current case, the final RMS deviation values of the backbone structures for all the simulation systems were no more than 0.8 nm (Figure 4), giving an indication that the receptor structures had reached to the equilibrium states with little alterations during the entire simulations.  To further study the conversed hydrogen bonding network mentioned above, 200 snapshots for the LY465608-bound and Cpd#1-bound states of PPARa, PPARb and PPARc receptors were retrieved from the last 1-ns segment on the molecular dynamics simulation trajectories with an interval of 5 ps, and statistical analyses were further performed on the conversed hydrogen bonding network (for detailed information, please see Table 2 and Table 3). For the PPARa receptor, the hydrogen bonds formed by the key residues (Ser280, Tyr314, Tyr464 and His440) and both ligands (LY465608 and Cpd#1) were quite similar, no matter the average distances between the heavy atoms of the hydrogen donor and receptor and the hydrogen bond lifetime. However, huge differences were detected for PPARb and PPARc receptors. For the PPARb receptor, although the differences between the hydrogen bonds formed by LY465608 and Cpd#1 were quite small, Cpd#1 could form an additional hydrogen bond with Thr288 in the Arm III region. For the PPARc receptor, the hydrogen bonds formed by Ser289 with LY465608 and Cpd#1 were quite different. Compared with the LY465608-bound system, the hydrogen bond formed by Ser289 in the Cpd#1-bound system employed a much smaller average distance between the heavy atoms and a much larger life-time.
This observation indicated that the hydrogen bond formed by Ser289 and Cpd#1 was stronger than that of the LY465608bound system.

AMDE Predictions
Absorption, distribution, metabolism and excretion (ADME) describes the druggable dispositions of the lead compounds in pharmacokinetics and pharmacology. In our study, some pharmaceutically relevant properties of the new designed agonist derivatives as well as the original LY465608 compound, such as the ''molecular mass'' (Mol_MW), ''hydrogen bond donors'' (DonorHB), ''hydrogen bond acceptor'' (AccptHB), ''partition coefficient'' (logP o/w), ''van der Waals surface area of polar nitrogen and oxygen atoms'' (PSA), and ''aqueous solubility'' (logS), were predicted by means of the QikProp module packaged in the Schrodinger 2009. The predicted ADME results were listed in Table 4.   According to the Lipinski's Rule of Five, the orally active drugs should have no more than one violation of the following criteria: i) no more than 5 hydrogen bond donors (DonorHB #5); ii) no more than 10 hydrogen bond acceptors (AccptHB #10); iii) molecular mass less than 500 Daltons (Mol_MW #500); and iv) the octanolwater partition coefficient logP in 20.4 to 5.6 range (20.4$ logP o/w #5.6). All the criteria for the original LY465608 compound were within the acceptable ranges of the Lipinski's Rule of Five, except the logP o/w value which was comparatively higher and had high-risk to induce unfavorable distributions of the compound on fat and body fluid. This also might be the reason why the original LY465608 compound shows comparatively high toxicity. However, after rational design with the core hopping approach, the ADME properties of the 7 candidates are all within the acceptable range for human beings, indicating that these derivatives found in the current study can be utilized as potential candidates for the purpose of developing novel drugs.
In summary, we employed the core hopping approach to generate a series of novel compounds with an aim of finding new and powerful pan-agonists for PPARa, PPARb and PPARc receptors. After flexible docking procedures and molecular dynamics simulations, a set of 7 novel compounds were found using LY465608, a typical PPAR pan-agonist, as a guiding structure. Compared with the existing pan-agonist LY465608, the 7 candidates found in the core hopping processes not only had the similar function in activating PPARa, PPARb and PPARc receptors, but also assumed the conformation more favorable in binding to the PPAR receptors. It is anticipated that the 7 candidates, especially Cpd#1, may become potential lead compounds for designing novel pan-agonist against PPARa, PPARb and PPARc receptors with comparatively lower toxicities. Table S1 The chemical structures of the top 7 hits in the core hopping and glide docking. For comparison, the typical PPAR pan-agonists bezafibrate, Ly465608 and GW677954 are also involved. (DOC)

Supporting Information
Author Contributions