Investigations on Binding Pattern of Kinase Inhibitors with PPARγ: Molecular Docking, Molecular Dynamic Simulations, and Free Energy Calculation Studies

Peroxisome proliferator-activated receptor gamma (PPARγ) is a potential target for the treatment of several disorders. In view of several FDA approved kinase inhibitors, in the current study, we have investigated the interaction of selected kinase inhibitors with PPARγ using computational modeling, docking, and molecular dynamics simulations (MDS). The docked conformations and MDS studies suggest that the selected KIs interact with PPARγ in the ligand binding domain (LBD) with high positive predictive values. Hence, we have for the first time shown the plausible binding of KIs in the PPARγ ligand binding site. The results obtained from these in silico investigations warrant further evaluation of kinase inhibitors as PPARγ ligands in vitro and in vivo.


Introduction
Peroxisome proliferator-activated receptors (PPARs) belong to the nuclear receptor super family and are ligand activated transcription factors, regulating the expression of a wide variety of genes. On activation by a ligand, they bind to the PPAR-responsive regulatory elements (PPRE) and/or PPAR associated conserved motif (PACM) as obligate heterodimers with retinoid X receptor (RXR) [1,2]. Similar to other nuclear receptor-family members, PPARs are multidomain proteins, consisting of an N-terminal transactivation domain (AF1), a highly conserved DNA-binding domain (DBD), and a C-terminal ligand binding domain (LBD) which has a ligand-dependent transactivation function (AF2) [3,4]. Three isoforms of PPARs (alpha, beta/delta, and gamma) have been identified so far in human, mouse, rats, xenopus, and hamsters [5][6][7] and among them, PPAR is the most intensively studied. PPAR has three alternatively spliced isoforms and all of them are expressed in adipose tissues [8,9].
It is primarily involved in the regulation of lipid metabolism and insulin sensitivity reactions and also plays an important role in carcinogenesis and cell physiology [10,11]. Also, PPARs have been shown to have ligand independent repression whereby they repress the transcription of direct target genes by recruitment of corepressor complexes which blocks the actions of coactivator complexes [12]. PPAR activation is involved in transcriptional regulation of genes involved in proliferation, angiogenesis, apoptosis, organogenesis, and energy metabolism and hence implicated in cell growth and viability [13][14][15][16]. PPAR signaling is modulated using different domains and various natural lipophilic agonists (ligands) such as unsaturated fatty acids, oxidized lipid species, eicosanoids, and prostaglandins [2,17,18]. Conformational changes caused by ligand binding lead to the modulation of PPAR activity by differential recruitment of cofactors [4,12]. PPAR exhibits high affinity towards thiazolidinediones (TZDs) [19]. TZDs including troglitazone, 2 PPAR Research rosiglitazone, and pioglitazone are FDA approved synthetic agonists of PPAR [20,21]. TZDs bind to the LBD of PPAR , activating the AF2 surface to accommodate the coactivators. The LBD of PPAR is largely a helical domain comprising 13 -helices (H1−H12, H2 ) and one parallel -sheet. The AF2 region (ligand-dependent transactivation function) is formed by helices 11 and 12 and consists of hydrophobic residues [22]. However, the use of TZDs has recently been restricted due to various side effects that include renal fluid retention, hemodilution, weight gain, edema, cardiomegaly and congestive heart failure, and loss of bone mineral density [23,24]. In view of the above reported undesirable side effects at therapeutic dosages, the focus has shifted to the discovery, development, and evaluation of "selective PPAR modulators" or SPPARMs as safer alternatives to PPAR full agonists. In contrast to the full agonists, the partial agonists show reduced transcriptional activity while having retained the insulin sensitization and hence show promising therapeutic potential with fewer side effects in animal models [25,26]. The acidic thiazolidinedione moiety of full agonists such as rosiglitazone forms strong hydrogen bonding network with the side chains of His323, His449, and Tyr473 from helices 5, 7, and 12, respectively, of PPAR and stabilizes AF2 to recruit coactivators [22]. However, partial agonists tend to stabilize the -sheet region by acidic substituents or form hydrophobic interactions and do not stabilize helix H12 via hydrogen bonding with Tyr473 [27,28].
On a separate note, PPAR undergoes several posttranslational modifications including phosphorylation of Ser273 by extracellular signal-regulated kinase ERK/cyclin-dependent kinase 5 (Cdk5) [29,30]. Moreover, the underlying mode of action for both full agonist and partial agonists to elicit antidiabetic property involves the inhibition of obesity-linked phosphorylation of Ser273 in PPAR . Hence rather than transactivation of the genes, the agonists cause conformational change in the LBD of PPAR preventing the kinase to phosphorylate the serine residue [29,30].
The associations of PPAR with signaling molecules including receptor and nonreceptor kinases corroborate the cross-talk function between the two signaling proteins [11,24,25,30,31]. Kumar et al. (2005) identified L-tyrosine derivatives as potential PPAR / inhibitors [32,33]. De Filippis et al. described the synthesis and the evaluation of PPAR activity of the new tyrosine derivatives, based on the combination of GW409544, a potent full agonist on both PPAR and PPAR and stilbene or phenyldiazene scaffolds [34]. Interestingly, some known ligands of PPAR (e.g., Honokiol, amorfrutin 1, amorphastilbol, and hydroxyhydroquinone) have tyrosine moiety like substructure and structural similarity to tyrosine kinase inhibitors. This led us to speculate on the possibility of TKI being ligands for PPAR . Structural superimposition of synthetic ligand rosiglitazone and selected TKI further confirmed our speculation. A question was posed whether these kinase inhibitors activate PPAR and could be potential PPAR ligands. The objectives of the present study were to (1) investigate the interaction of selected kinase inhibitors such as ibrutinib, sorafenib, sunitinib, erlotinib, gefitinib, and dabrafenib with PPAR in silico (Figure 1), (2) a comparative analysis of interaction of these KIs with rosiglitazone in the LBD of human PPAR using molecular docking studies, and (3) molecular dynamic simulation and MM/PB (GB) SA studies to evaluate the stability and conformational changes due to the interaction of the kinase inhibitors in the PPAR binding site.

Protein and Ligand Preparation.
The starting structure for the simulations was taken from the X-ray structure of the ligand binding domain and coactivator assembly of PPAR (PDB code 2PRG, resolution: 2.3Å) [21]. The protein was prepared using Schrodinger's protein preparation wizard [35] by removal of crystallographic water molecules and addition of hydrogen atoms, followed by minimization and optimization using OPLS2005 force field of Schrodinger [36]. The shape and properties of the receptor were represented on a grid by several different sets of fields that provide progressively more accurate scoring of the ligand poses. The SDF files for the drugs rosiglitazone, ibrutinib, sorafenib, sunitinib, erlotinib, dabrafenib, and gefitinib were obtained from Pub-Chem database (https://pubchem.ncbi.nlm.nih.gov/). These molecules were then prepared in Schrodinger Ligprep wizard. During the ligand preparation all possible conformations were taken into account. The ligands were subjected to further predocking preparations where hydrogens were added followed by minimization and optimization in OPLS 2005 force field as implemented on Maestro software. Finally, 32 conformations of each ligand were generated and used for docking.

Docking of Ligands in PPAR Binding Domain Using
Glide. Molecular docking procedures were carried out after preparing the ligand library using Schrodinger Ligprep module and defining the grid corresponding to the ligand binding site of the protein. The grid was prepared using rosiglitazone at the center and the interacting residues as the ligand binding site residues along with a cubic space of 12 angstroms around the ligand rosiglitazone. The Glide module of Schrodinger uses Systematic and Simulation Method for searching flexible ligand poses. In a systematic method, it uses incremental construction for searching, and its output -Score is an empirical scoring function which is a combination of various parameters.

Molecular Dynamics
Simulations. The poses with highest Glide scores obtained from the docking simulations (proteinligand complexes) were further subjected to MD simulations. The purpose of performing the MD simulations was to determine the stability of the drug molecules and to understand the binding pattern and to identify the binding residues in the receptor to the drug molecule. Parameters for all the drug molecules (bound at the ligand binding site) were generated using antechamber module of AMBER suite [37,38]. The restrained electrostatic potential (RESP) was used to describe the partial atomic charges. Then the general AMBER force field (GAFF) [39] was used to describe the parameters of drug molecules. The standard AMBER force field for bioorganic systems (ff99SB) was employed to describe the protein, followed by the addition of hydrogen atoms and counterions to neutralize the system. The input files for energy minimization, dynamics, and analysis were prepared with xleap. Both systems were solvated using atomistic TIP3P water in a box with edges at least 12Å from the complex. All simulations were performed using AMBER molecular dynamics suite version [37]. Energy minimization was first conducted with the steepest descent method and then switched to conjugate gradient every 500 steps for a total of 5000 steps with 0.1 kcal/molÅ 2 restraints on all atoms of the complexes. Following this step, another two rounds of energy minimization were performed by only restraining the protein and further releasing all the restraints for 2000 steps of each round. Long-range Coulombic interactions were handled using the particle mesh Ewald (PME) summation. For the equilibration and subsequent production runs, the SHAKE algorithm was employed on all atoms covalently bonded to a hydrogen atom, allowing for an integration time step of 2 fs. The system was gently annealed from 0 K to 300 K over a period of 50 ps using a Langevin thermostat with a coupling coefficient of 1.0 ps and 50 ps of density equilibration with weak restraints. The system was again equilibrated for 500 ps without any restraints. The production phase of the simulations was run without any restraints for a total of 25 ns on each system. Coordinates and energy values were collected every 10 ps throughout the simulations.

Binding Free Energy
Calculations. The binding free energies of PPAR for all the KI molecules were analyzed by the MM/PB (GB) SA scripts, integrated in the AMBER 12 software package. In this procedure, snapshots were first extracted from the obtained trajectories. For each snapshot, free energy is calculated for the protein, ligand, and complex using single trajectory approach. The binding free energy was computed as the difference:

Per Residue Interaction Decomposition.
To determine the contribution of each residue to the binding energy, the MM-GBSA method was used. MM-GBSA method decomposes the interaction energies for each residue by considering molecular mechanics and solvation energies without consideration of the contribution of entropies. Each residue contribution includes three terms: van der Waals contribution (Δ vdw ), electrostatic contribution (Δ ele ) in a vacuum, and solvation contribution (Δ solvation ).
All energy components in the above equation were calculated using 1000 snapshots from the last 10 ns of the MD simulation. The calculations and the computational methods used in this paper are well documented in the literature and have been used previously for small molecules [40].

Molecular Docking Studies.
A ligand is stabilized energetically at the binding site in the protein structure by weak intermolecular interactions such as hydrogen bonds and hydrophobic interactions [31]. In the present study, the binding mode for KIs in the ligand binding site of human PPAR was investigated using molecular docking studies. The docking results and the docked conformations of KIs are illustrated in Table 1 (Table 1) forms H-bond interaction between the oxygen atom of diphenyl ether side chain and the phenolic hydroxyl group of Tyr327 (helix 5) of protein (Figure 2(b)). The interaction of ibrutinib is stabilized by hydrophobic interactions with -sheet and turn residues Ile341, Ser342, and Met348 (Table 1, Figure 2(b)). Molecular overlay of binding pose for ibrutinib/PPAR in the surface volume of rosiglitazone/PPAR indicates that the conformation of ibrutinib in the PPAR binding pocket is similar to that of rosiglitazone (Figure 3(h)). Gefitinib/PPAR docked complex with Glide XP score −9.10 kcalmol −1 and Glide energy −52.45 kcalmol −1 shows H-bond interactions between fluorine substituent of 3-chloro-4-fluorophenyl group of the ligand and main chain nitrogen atom of Glu343 ( -turn) of the protein. Gefitinib forms another H-bond interaction between NH linker connecting the quinazoline ring and 3-chloro-4-fluorophenyl ring and the main chain carbonyl oxygen atom of Leu340 of PPAR (Figure 2(c)).
The docking conformation of gefitinib in the PPAR binding site is stabilized by hydrophobic interactions (Table 1, Figure 2(c)). The orientation of gefitinib in the PPAR binding pocket is more towards the -sheet region (Figure 3(c)). Figure 2(d) and Table 1 211  215  219  223  227  231  235  239  243  247  251  255  259  263  267  271  275  279  283  287  291  295  299  303  307  311  315  319  323  327  331  335  339  343  347  351  355  359  363  367  371  375  379  387  391  395  399  407  411  415  419  423  427  431  435  439  443  447  451  455  459  463  467  471  475  different physiological conditions. In order to evaluate the possible deviation in the structure during simulation, RMSD were calculated based on the initial backbone coordinates of the protein-ligand complexes. The plot for RMSD indicates that stability in the interactions for all the protein-ligand complexes is attained after 15 ns of simulation. The most stable interaction was observed for PPAR -gefitinib complex as it reached the equilibrium faster with lowest RMSD level. The curves for PPAR -sorafenib complex were quite similar to the standard PPAR -rosiglitazone, showing stable RMSD after 10 ns. PPAR -ibrutinib complex shows a high RMSD value initially attains stable conformation after 10 ns at a slightly higher level than PPAR -rosiglitazone complex. This is followed by PPAR -dabrafenib and PPAR -erlotinib complexes that have highest initial RMSD values and attain stability after 15 ns. Moreover, after 15 ns the average RMSD for all the complexes was approximately 1.5Å to 2.5Å and remains stable throughout the 25 ns simulation. However, it is interesting to observe that there is a sudden decrease in the stability of PPAR -rosiglitazone complex around 20 ns with high value of RMSD (2.5Å) for the complex. The results from the RMSD plot indicate that the PPAR -KI and the standard PPAR -rosiglitazone systems could be satisfactorily explored in the allocated nanosecond simulation studies. The solvent accessible surface area (SASA) plot for PPAR -KI and the standard PPAR -rosiglitazone complexes show no significant changes in the equilibration curve throughout 20 ns MD simulation run (Figure 4(b)). This indicates that surface of all the protein-ligand complexes has maintained a similar accessibility to the solvents.

PPAR Research
structural stability is observed for all the PPAR -ligand complexes around the ligand binding site residues including Cys285, Gln286, Ser289, His323, Leu330, Met348, and Tyr473.  [41]. The main objective of this method is to find the difference in the free energies of bound and unbound state of protein-ligand complexes. All the thermochemical properties were calculated by MM-PBSA and MM-GBSA approach using AMBER program, for each coordinate at every 10 ps sampling frequency throughout the MD trajectory for all the protein-ligand complexes. The complexes with lowest binding energy are considered to be most stable ( Table 2). The total free energies (Δ TOT ) obtained from MM-GBSA for the protein-ligand complexes show comparable values ranging from −49.82 kcalmol −1 to −33.10 kcalmol −1 (−49.82 kcalmol −1 for PPAR -erlotinib complex and −33.10 kcalmol −1 for PPAR -sorafenib complex) ( Table 2). Van der Waal's energy (Δ VDW ), electrostatic energy (Δ ELE ), nonpolar contribution (Δ NPOL ), and total energy of solute (Δ GAS ) have negative values and show favorable contribution to the total free energy. However, total energy of solvation (Δ SOLV ) and polar solvation contribution (Δ POL(GB) ) possess positive values and therefore contribute unfavorably towards the total free energy. Table 3 shows total free energies (Δ TOT ) obtained from MM-PBSA approach for the protein-ligand complexes.  PPAR -erlotinib complex appears to be the most stable complex having a binding energy value of −72.13 kcalmol −1 , followed by PPAR -ibrutinib complex (−62.62 kcalmol −1 ). PPAR -gefitinib complex possesses the lowest free energy value (−39.54 kcalmol −1 ). Further, the entropy contribution (Δ ) was calculated using quasi-harmonic approximation ( Table 4). The values for entropy contribution were very similar, indicating no clear involvement of this term in the determination of free energy. Here it is important to mention that it has been suggested that "entropy contributions can be neglected if only a comparison of states of similar entropy is desired such as two ligands binding to the same protein" [42]. Hence the entropy calculations can be ignored.
purpose, Generalized Born (GB) model was used to efficiently compute solvation free energy. The individual residue contributions to the binding free energy vary in the range of +1.0 to −10.0 kcal/mol. Among these 270 residues in the crystal structure of PPAR (PDB ID: 2PRG), seven residues (Cys285, Arg288, Glu295, Ile326, Leu330, Ile341, and Ser342) have binding free energy contribution (Δ res ) more than −2.0 kcal/mol ( Figure 6). These residues are thus predicted as energetically important for ligand binding inside the ligand binding site via hydrophobic or hydrogen bond interactions in almost all complexes.

Conclusions
In the present study, we for the first time investigate the role of kinase inhibitors (KIs) as PPAR ligands. The molecular docking studies show that KIs could bind effectively in the PPAR pocket by forming important H-bond contacts and hydrophobic interactions, similar to known agonists of PPAR . The stability of interaction of KIs and PPAR ligand binding site residues predicted by docking experiments was further assessed by molecular simulation studies. Further, the free energy and entropy calculations provide insight into the binding affinity of KIs towards PPAR . Overall, we propose KI as promising lead compounds that could be PPAR ligands. We expect that the structural insights obtained in this study will facilitate the design of novel KI based PPAR gamma ligands. More elaborate studies are needed for validation of these predictions and better understanding of the mechanism of the interaction of KIs with PPAR and the future clinical applications of KIs.