Probing the intermolecular interactions of PPARγ-LBD with polyunsaturated fatty acids and their anti-inflammatory metabolites to infer most potential binding moieties

Background PPARγ is an isoform of peroxisome proliferator-activated receptor (PPAR) belonging to a super family of nuclear receptors. PPARγ receptor is found to play a crucial role in the modulation of lipid and glucose homeostasis. Its commotion has been reported to play a significant role in a broad spectrum of diseases such as type 2 diabetes mellitus, inflammatory diseases, Alzheimer’s disease, and in some cancers. Hence, PPARγ is an important therapeutic target. Polyunsaturated fatty acids (PUFAs) and their metabolites (henceforth referred to as bioactive lipids) are known to function as agonists of PPARγ. However, agonistic binding modes and affinity of these ligands to PPARγ are yet to be deciphered. Methods In this study, we performed a comparative molecular docking, binding free energy calculation and molecular dynamics simulation to infer and rank bioactive lipids based on the binding affinities with the ligand binding domain (LBD) of PPARγ. Results The results inferred affinity in the order of resolvin E1 > neuroprotectin D1 > hydroxy-linoleic acid > docosahexaenoic acid > lipoxin A4 > gamma-linolenic acid, arachidonic acid > alpha-linolenic acid > eicosapentaenoic acid > linoleic acid. Of all the bioactive lipids studied, resolvin E1, neuroprotectin D1 and hydroxy-linoleic acid showed significant affinity comparable to proven PPARγ agonist namely, rosiglitazone, in terms of Glide XP docking score, H-bond formation with the key residues, binding free energy and stable complex formation with LBD favouring co-activator binding, as inferred through Molecular Dynamics trajectory analysis. Conclusion Hence, these three bioactive lipids (resolvin E1, neuroprotectin D1 and hydroxy-linoleic acid) may be favourably considered as ideal drug candidates in therapeutic modulation of clinical conditions such as type 2 DM, Alzheimer’s disease and other instances where PPARγ is a key player.


Background
Peroxisome proliferator-activated receptor (PPAR) comprises of three isoforms, α, β and γ belonging to a super family of nuclear receptors [1]. PPARs are ligand-activated transcription factors that regulate genes playing a vital role over a broad spectrum of physiological and pathological conditions [2]. PPAR receptors are expressed by various tissues including muscles, hepatocytes, adipocytes and endothelial cells. Though the three isoforms of PPAR (α, β and γ) share high level of sequence and structural similarity they are distinct in terms of expression and tissue distribution [3]. Molecular 3D structure of PPAR constitutes DNA binding domain at the N-terminus and ligand binding domain (LBD) at the C-terminus. PPAR's interaction with its agonist leads to heterodimer formation with retinoid X receptor (RXR) [4]. PPAR-RXR heterodimer gets bounded at peroxisome proliferator response elements (PPREs) occupying the promoter region of target specific genes. Further, this process leads to recruitment of various transcriptional cofactors involved in the initiation of transcription process, thereby, triggering expression of several genes involved in diverse physiological and pathological processes [5][6][7]. Each PPAR subtype plays a unique physiological role in different tissues; however, all the three isoforms are well known to be involved in lipid and glucose homeostasis [8]. Of all the three isoforms, PPAR α and γ are most extensively studied when compared to PPAR β.
PPARα is predominantly expressed in tissues involved in metabolic activities of various tissues including muscles, heart, liver, intestine and brown adipose cells. PPARα activation leads to a decrease in lipid levels. PPARα receptor functions as a lipid sensor and helps in controlling energy combustion [9][10][11][12]. PPARγ is widely expressed in adipocytes, thereby, playing a crucial role in adipogenesis, lipid synthesis and in maintaining energy balance. PPARγ activation improves insulin sensitivity. In addition, PPARγ is expressed in spleen, large intestine, white and brown adipose tissues that may account for the involvement of these tissues in the pathobiology of type 2 Diabetes Mellitus and metabolic syndrome [10] [13,14]. PPARβ is abundantly expressed in liver and abdominal adipose tissues by which it regulates blood cholesterol, glucose levels and influences fatty acid oxidation in cardiac and skeletal muscles [15,16].
In the present study, we focussed on PPARγ receptor, as it plays a critical extensive role over broad spectrum of diseases such as type 2 diabetes mellitus, inflammatory diseases, Alzheimer's disease, and in some cancers [17][18][19][20][21][22][23]. PPARγ comprises a Y-shaped ligand binding domain (LBD), which is segmented into three arms, arm I, arm II and arm III. Arm I is extended towards helix12 (H12), known to be polar and widely conserved across the PPAR isoform [24][25][26]. It also harbors transcription activation function-2 (AF-2) at C-terminal region and is held in its active conformation by the hydrogen bonding network with Arm I favouring ligand binding [27,28]. Arm II and Arm III are found to be less conserved compared to Arm I and are hydrophobic in nature [29]. It has been proposed that diverse ligands bind to PPARγ with different binding modes to LBD with most of the ligand binding scenarios displaying a hydrophilic interaction with Arm I region and hydrophobic interactions with either arm II or arm III regions [30].
It has been proposed that full agonists of PPARγ reside over a large area of LBD of PPARγ with a U-shaped conformation and ideally comprising a polar head and a hydrophobic tail. The polar head of the full agonist forms a network of hydrogen bonds with the ARM-I residues (His449, His323, Ser289, Tyr327 and Tyr473) of PPARγ side chains. The hydrogen bonds formed with these residues are responsible for the conformational change of H12 and activation of PPARγ activity [31]. In contrast, partial agonists activate PPARγ by an H12 independent mechanism, wherein, the key residues in LBD are completely different to that of the full agonists in that it leads to decrease in H12 stability, thereby affecting the coactivators binding, which, in turn, reduces the transcriptional activity of PPARγ [32,33]. Most of the previous studies suggested that partial agonists form hydrogen bond with Ser342 of LBD [31,33].
PPARγ agonists reduce lipid levels, enhance insulin sensitivity and thus, show anti-diabetic and antiinflammatory actions. PPARγ, is well documented to be activated by a wide range of fatty acid molecules and their metabolites, of which PUFAs and its metabolites play a major role in exerting beneficial effects [34]. PUFAs are also reported to play significant roles in inflammatory and immune responses [35], lowering the levels of total cholesterol, Triglycerides and Low-density lipoproteins (LDL) [36,37], in mediating apoptosis in colon cancer cells [38] and also shown to reduce the risk of early atrial fibrillation during post cardiac surgery [39]. Molecular modelling and structural bioinformatics studies serve as powerful and efficient tools for studying intermolecular interactions and dynamic behaviour of molecules in an in silico simulated conditions [40][41][42]. In the present study, we performed molecular docking of PPARγ against various PUFAs and its metabolites (henceforth called as bioactive lipids) [43], keeping cocrystal bound Rosiglitazone (BRL) as a reference structure (PDB ID: 4OF8). Further, the bound complexes were subjected to molecular dynamics simulation to compare the binding efficacies and to infer the agonistic binding modes of fatty acids, which are potential therapeutic molecules [44,45], towards identification of the most potent and efficient bioactive lipid agonist targeting PPARγ (Fig. 1).

Protein Preparation
As a preliminary step, the crystal structure of PPARγ in complex with Rosiglitazone (PDB ID: 4O8F) at a resolution of 1.9 Å was pre-processed using Protein preparation wizard of Schrödinger suite towards optimizing the stereochemistry by assigning proper bond order, removing steric clashes, adding hydrogen atoms, fixing the disulphide bonds, missing residues, atoms and loops. Further refinement was performed by adjusting the terminal chi rotation of Asparagine, Glutamine and Histidine residues. Optimal protonation states for Histidine residues were also assigned followed by the removal of unwanted hetero groups. Finally, energy minimization was performed with OPLS2005 to obtain the optimal geometry favourable for the commencement of plausible docking studies.

Receptor grid generation
The receptor grid was generated using the grid generation option of the Schrodinger suite across the PPARγ LBD domain to perform a targeted docking study. The grid was fixed across residues that play a critical role at the three different arms of the LBD by setting the vanderwaals radii for the receptor with a scaling factor of 1.0 Ǻ and a partial cut off of 0.25 Ǻ.

Ligand Preparation
To expedite the protein-ligand docking workflow, all the ligand molecules were optimized using Schrodinger Ligprep module by fixing the vanderwaals radii with the scaling factor of 0.80 Ǻ and a partial charge cut-off of 0.15 Ǻ. Finally, the optimized compounds were energy minimized with OPLS 2005 as force field.

Receptor-Ligand Docking
Docking studies establish interactions between protein and ligand molecules, thereby, aids in identifying the most favourable binding pose forming a stable complex with significant binding affinity, as scored by the docking score [46]. In this study, protein-ligand docking was performed using Schrodinger Glide version [7.0] across the grid set over the LBD region. Here, a flexible docking was performed against 10 ligand molecules, in which the best binding pose among 100 generated poses for each molecule was determined based on the Glide XP docking score. A rigid Re-docking was also performed for the rosiglitazone (BRL) against its bound crystal structure by setting the grid across the LBD so as to infer the predictive efficacy.

MMGBSA Scoring
Molecular Mechanics Generalized Born Surface Area (MMGBSA) scoring was also performed for all the docked complexes to calculate the binding free energies by implementing the equation through Schrödinger Prime Module: Where electrostatic solvation energy (ΔG GB ) was calculated by using the GB Models, ΔE MM is the difference between the minimized energies of ligand-protein complex and the total energies of protein and ligand in free form. ΔG solv is the difference in the GBSA solvation energies of the ligand-receptor complex and the sum of the solvation energies of receptor and ligand in the unbound state. ΔG SA is the difference in the surface area energies for the free receptor and the ligand. It is used to identify plausible binding conformations among the docked complexes in terms of binding free energy towards further stringent ranking of the complexes [47].

Molecular Dynamics simulation
Molecular dynamics simulation was performed for apo, co-crystal structure (PDB id: 4O8F) and for all the docked complexes using Desmond 3.6. The system was built using a cubical box solvated with Simple Point Charge (SPC) water model. Subsequently, the system was neutralized by adding 4 Na + ions for apo and 5Na + ions for complexes at a concentration of~6.22 mM. Further, this system was energy minimized with OPLS2005 as force field [48]. SHAKE algorithm was applied to restrain the geometry of water molecules, bond lengths and angles of heavy atoms and to constrain covalent bonds during MD simulation [49]. Periodic Boundary Conditions (PBC) were applied to stimulate a continuous system [50] and Particle Mesh Ewald method (PME) was applied for long range electrostatics [51]. Further, the system was equilibrated with NPT ensemble by setting temperature and pressure parameter to 300 K and 1.0 Bar, respectively. Nose-Hoover chain and Martyna-Tobias-Klein was chosen as a coupling algorithm for temperature and pressure, respectively [52,53]. Further, the equilibrated system with a total of 48,294 atoms was simulated for 20 ns (nanosecond) with a time step of 2 fs (femtosecond) and trajectories were recorded after every 1.0 ps. The Root Mean Square Deviation (RMSD) was calculated for the backbone atoms and were graphically analysed at a time point scale [54,55]. Similarly, root mean square fluctuation (RMSF) for each residue was also calculated to compare the major conformational changes in the residues between apo form and docked complex forms by keeping the rosiglitazone (BRL) bound crystal structure as a reference [56]. The radius of gyration was also calculated to infer the compactness of protein-ligand complexes for the comparison with apo form [57]. The 2D inter-molecular interaction plots depicting the complex stability during the MD run was also generated to infer the stability of all the complex structures.

Results
Comparative Molecular Docking studies of PPARγ with bioactive lipids Molecular Docking was performed for all the 10 bioactive lipid compounds against the PPARγ-LBD. Glide XP dock score and MMGBSA binding free energy score were calculated for all receptor-ligand docked complexes, which revealed the scores to be in the range of −4.8 to −9.9 kcal/mol and −71.147 to −106.046 kcal/mol (Table 1), respectively. Further, 2D interaction maps with a cut-off of 4 Å for each docked complex was generated to visualize the intermolecular interactions. This inferred that all the compounds to be majorly stabilized by hydrogen bonds during complex formation (Fig. 2). The re-docking of Rosiglitazone to PPARγ-LBD also showed an agreeable deviation of 0.4 Å, inferring the predictive accuracy. The 2D maps were further scrutinized for ligand contacts with the hotspot residues spanning Arm-I, Arm-II, Arm-III and AF-2 domain towards classifying potential full and partial agonists.

Molecular dynamics simulation of the docked complexes
An unrestrained molecular dynamics simulation study of apo form and all docked complexes were performed to infer the backbone stability, residue fluctuations, structural compactness measured in terms of RMSD, RMSF, and Rg, respectively. The RMSD trajectory revealed that all complexes to be stable during the entire production run with the system convergence at~15 ns (Fig. 4). The RMSF trajectory inferred maximum fluctuations at L1 (238-251) and L2 (260-276) regions across all the apo and docked complexes (Fig. 5a). Rg trajectory also revealed structural compactness in apo and protein-ligand complexes (deviation within 1 Å) (Fig. 5b). A trajectory of Intra-molecular hydrogen bond counts (mean average of 210 H-bonds) was plotted for the entire production run. As expected, there was a gain in h-bonds in docked complexes in comparison with holo forms (Fig. 5c). For the top three ranking hits, Protein-ligands contact bar charts were plotted, which inferred the contribution of hotspot residues in establishment of intermolecular contacts, and was majorly found to be h-bond mediated (Fig. 6).

Receptor-Ligand Docking analysis
Molecular Docking studies were subsequently performed for all the bioactive lipids against the assigned grid surface on the protein. The docked complexes were analysed for the Glide XP score, MMGBSA score and Hbond interactions, to collectively infer the ligand binding affinity [58,59]. The Glide receptor-ligand docking results are tabulated in (  (Table 1). RsvE1, NPD1 and H-LA formed hydrogen bonds with the arm-I residues Hie449 (H11), Ser289(H3) and Tyr327, which are reported to be important for producing full activity of the compound by direct stabilization of H12 helix and are responsible for the transactivation activity of PPARγ [24,31,60,61]. HLA also formed a hydrogen bond with Tyr473 (H12) located on the arm I that harbours transcription activation function-2(AF-2), which is obligatory for ligand binding and as well as in complimenting PPARγ function [27]. Tyr473 is also hypothesized to play a crucial role in AF-2 stabilization, as it occupies helix12 which closes the ligand binding site upon ligand binding. This activity of helix12 favours in the reduction of conformational fluctuations that sets an optimal LBD structure for co-activator binding. It has also been reported that mutation at Tyr473 leading to the loss of agonistic activation by the ligands.
Alpha-linoleic acid (ALA) showed a Glide XP dock score of −5.174 kcal/mol by forming a single H-bond with the side chain of Hie449 (H11). Arachidonic acid (AA) and docosahexaenoic acid (DHA) formed a single H-bond with the side chain of Hie323 with an XP dock score of −7.040 kcal/mol and −7.925 kcal/mol, respectively. Eicosapentaenoic acid (EPA) and gammalinoleic acid (GLA) displayed a XP dock score of −7.126 kcal/mol and −5.068 kcal/mol by forming two hydrogen bonds with the side chain residues of Hie323 and Ser289, respectively. Lipoxin A4 (LXA4) with a docking XP score of −5.796 kcal/mol, displayed two hydrogen bond interactions with the active site residues Hie 449 and Tyr327. Linoleic acid (LA) showed the least XP dock score of −4.820 kcal/mol among all the docked bioactive lipids compounds by forming an H-bond with the back bone of Glu342 (Table 1).
The ligand interactions with ARM-I side chains of Hie323 (H3), Hie449 (H11), Ser289(H3) and Tyr327 of LBD is reported to play a significant role in enhancing agonistic activity as it is found to indirectly influence the H12 stabilization [24,31,60,61]. Of the 10 docked complexes, 9 bioactive lipids are predicted to act as full agonists for PPARγ by forming an H-bond interaction with the ARM-I hydrophilic cavity residues (Hie 323, Hie 449, Tyr473, Ser289 and Tyr327). Linoleicacid has been found to be a partial agonist for PPARγ as it formed H-bond interactions with Glu343 of ARM-III [24]. The superimposition of all ten ligand molecules (bioactive lipids studied in the present study) along with the rosiglitazone as a reference ligand, inferred all the ligands to be well bound within the three arms of agonist binding domain by interacting with key residues (Figs. 2 and 3) [32,62].

Molecular dynamics studies of the protein-ligand complexes
The RMSD trajectory of the backbone atoms of apo form and all the complex structures showed a convergence after a time frame ranging from~13,000 ps to1 5,000 ps with a maximum mean value of~3.5Ǻ and S.D of 0.675Ǻ (Fig. 4). The RMSF fluctuation has reached the maximum peak limit of 6-7Ǻ at the residues residing over the lengthy loop region (Loop1 with 14 residues and Loop2 with 17 residues) and remaining   residues showed fluctuations within 1Ǻ. The H12 (Tyr473) residue showed a slightly higher fluctuation in the apo form compared to all the other protein-ligand complexes. This corroborates with the hypothesis that upon ligand binding to H12 hotspot to confer closure of LBD, thereby stabilizes the LBD conformation in the ligand bound state (Fig. 5a) [27]. The apo form and all the 10 (bioactive lipids) complexes were found to be compact; by displaying a radius of gyration less than~1.3 (Fig. 5b). The intra-molecular H-bond graph also inferred no major loss in the protein secondary structure conformations, thereby, reinforcing that the proteinligand complex structures to be well stabilized (Fig. 5c). However, the protein complexed with the ligands Rsv E1, NPD1 and H-LA were found to establish strong protein-ligand contacts by constituting plausible intermolecular H-bonds with the active site residues of the LBD domain throughout the period of simulation. RsvE1, NPD1 and H-LA displayed a consistent and plausible H-bond interactions with the critical residues of ARM-1 (Hie323, Tyr327, Hie449, Lys367 and Tyr473) for~50% during the production run (Fig. 6a-c). The three bioactive lipids (RsvE1, NDP1 and H-LA) interacting with two conserved side chain amino acids Tyr327 and Lys367, are reported for favouring the formation of H-bonds at the Keto-group of the ligands. These amino acid side chains are known to play a vital role in determining the specificity of ligands that can efficiently couple with the receptor [62] (Fig. 6a and b). It is noteworthy that these residues of arm-1 (His323, His449, Tyr473, Tyr327 and Lys367) are hypothesized for indirect and direct stabilisation of H12 (Tyr473), thereby favouring for co-activator binding and for contributing towards the transactivity of PPARγ.

Conclusion
In the present study, we performed a molecular docking analysis with bioactive lipid compounds which are documented to be agonists of PPARγ [30]. All the 10 docked complexes were found to be well bound within the three different arms of LDB. Thereby in-order to perform a further comparative study on the protein-ligand interaction stability and the conformational stability of the LBD, molecular dynamics studies were carried out for apo form, LBD co-crystallized with rosiglitazone and for all the docked complexes with bioactive lipids.
The cumulative analysis based on the Glide score, Prime-MMGBSA free binding energy score and MD trajectories analysis infer that RsvE1, NDP1 and H-LA to be the best docked compounds as these three showed relatively a significant score to that of the re-docking score of the reference ligand Rosiglitazone in terms of glide XP score, MMGBSA score and H-bond forming pattern with the key residues ( Table 1). The MD trajectory also reinforce that RsvE1, NDP1 and H-LA to be the best compounds based upon the RMSF graph (Fig. 5a), which clearly depicts that upon RsvE1, NDP1 and H-LA ligands binding to the LBD, it displays a much lesser fluctuation at H12 region in comparison to that of the reference ligand Rosiglitazone (BRL) and thereby maintaining a stable LBD structure for coactivator binding.
The protein-ligand contacts of the RsvE1, NDP1 and H-LA with PPARγ throughout the production run of 20 ns also corroborates well by displaying a strong intermolecular interaction with the key active site residues of the LBD domain over~50% of the simulation time, which confirms the stability upon the ligand binding and its intactness with the receptor.
Moreover, RsvE1 is found to mimic the insulin sensitizing and anti-steatotic activities of omega-3-PUFAs, and also found to induce adiponectin expression similar to rosiglitazone [63]. NPD1 interacting with the PPARγ is hypothesized to be play a significant role in suppressing neuroinflammation and thus, is considered to be of benefit in some neurodegenerative diseases [64]. Hydroxy-linoleic acid (H-LA) is reported to have anti-proliferative action on tumor cells [65]. In the present study, we have ranked all these PUFA compounds collectively based on the Glide docking score, H-bonding formation with the key residues, Prime MMGBSA score and MD simulation trajectory analyses and observed that their relative PPARγ agonist activity to be as follows: RsvE1 > NPD1 > H-LA > DHA > LXA4 > GLA > AA > ALA > EPA > LA. However, Further experimental studies need to be performed to verify the results obtained in the present bioinformatics analysis. Such studies could include the assessment of ability of these bioactive lipids in the pathobiology of obesity, insulin resistance, type 2 DM, metabolic syndromes and cancer.