Molecular Dynamics Free Energy Simulation study to Investigate Binding Pattern of Isoliquiritigenin as PPAR 𝛾 Agonist

Phytochemicals are rich source of bioactive constituents and can be used as another alternative to currently used drugs for diseases like Diabetes mellitus. The potential of Isoliquiritigenin (a constituent of Pterocarpus marsupium) as PPAR 𝛾 agonist was evaluated by in silico technique. Autodock results showed that Tyr327, and Tyr473 of the PPARγ forms H-bonds with Isoliquiritigenin (binding energy of -7.46 kcal/mol) and Troglitazone (known drug) showed H bond with Tyr327, Ser289, with binding energy of -11.01 kcal/mol. Isoliquiritigenin, binding energy in Extra precision (XP) was -6.74 kcal/mol while Troglitazone docking, gave binding energy in XP mode as -9.59 kcal/mol. The best Induced t docking (IFD) score of the optimised PPARγ- Isoliquiritigenin complexes was -9.39 Kcal/mol. The important residues in IFD forming H bond were Cys 285, Arg 288, Tyr 327 and Leu 340. The post docking MM/GBSA free energy for PPARγ with Isoliquiritigenin and Troglitazone was -49.29 and -71.48 Kcal/mol respectively. Binding interaction in MD simulation and Principal Component Analysis studies revealed stable binding throughout 100 ns simulation. Post Simulation MM/PBSA free energy was calculated. The results indicated that compound possessed a negative binding free energy with -114.37KJ/mol. It was observed that van der Waals, electrostatic interactions and non-polar solvation energy negatively contributed to the total interaction energy while only polar solvation energy positively contributed to total free binding energy. The Isoliquiritigenin fulls the criteria of drug-likeness property. Thus, study presents a systematic analysis on molecular mechanism of action of Isoliquiritigenin as PPARγ agonist in controlling Diabetes mellitus.


Introduction
Diabetes mellitus (DM) is a type of metabolic disorder that causes hyperglycemia. DM is classi ed in two main type, type 1 and type 2. Type 1 DM is the consequence of complete or near-total insulin de ciency this might be due to interactions of genetic, environmental, or immunologic factors that lead to the destruction of the pancreatic beta cells and insulin de ciency. Type 2 DM can be classi ed as group of disorders occurred due to insulin resistance, reduced insulin secretion, and increased glucose production and unfamiliar fat metabolism. Diabetes is a fast growing chronic metabolic disorder around the world. and A and B, domains involved in ligand-independent coregulator binding. The C domain is the DNA binding domain and is highly conserved. Two highly conserved zinc ngers are involved in the recognition of speci c DNA half-sites termed peroxisome proliferator response elements (PPRE). The E domain is ligand binding site. Ligand binding at E domain facilitates the interaction with coregulator molecules which leads to chromatin remodelling and activation of transcriptional machinery [1,2,3,4,5,6,7] .
Computational approach in nding out potential new drug and targets is less time consuming and cost effective. Molecular docking and simulation help in deducing detailed information of the behaviour of protein and ligand-protein complexes at an atomic level and is the most powerful tool to study protein exibility. The binding free energy, △Gbind, docking and scoring are generally used in drug design although it is not accurate methods to differentiate between binders and non-binders. MM/PBSA (molecular mechanics [MM] with Poisson--Boltzmann [PB] and surface area solvation) and GBSA, are group of methods which is inexpensive and more precise than the scoring functions [8,9,10,11] . The therapeutic property of plants is because of bioactive chemicals that regulates physiological property in the human body and also provide protection against various diseases. Chronic disorder such as Diabetes mellitus needs medication generally for longer periods and currently used drugs have side effects therefore, plant derived product may be a safer alternative. Pterocarpus marsupium use in lowering of blood glucose is an age long traditional practice in India. It also possesses beta cell protective and regenerative properties. The heartwood yields liquiritigenin, Isoliquiritigenin, alkaloid (0.017%), and resin (0.9%). Ethyl acetate extract of powdered dried heartwood of P. marsupium revealed the presence of following constituents: (-) epicatechin (a avonoid), pterosupin (a dihydrochalcone), marsupin (a benzofuranone), pterostilbene, liquiritigenin (a stilbene), isoliquiritigenin, (2S)-7-hydroxy avanone, 7, 4'dihydroxy avone, p-hydroxybenzaldehyde, (2R)-3 -(p-hydroxyphenyl)-lactic acid, and pm-33. Several studies have been conducted on various animal species viz., rats, dogs, and rabbits to study the hypoglycemic effect of P. marsupium. There are reports that P. marsupium restores the normal insulin secretion by reversing the damage to the beta cells and by repopulating the islets [12,13,14,15,16,17,18,19] . The aim of this study is to high light the molecular interaction properties of Isoliquiritigenin and PPARγ by molecular docking, dynamics simulation and binding energy calculation. The drug likeliness of Isoliquiritigenin was also performed.

Molecular Docking
PPARγ was obtained from RCSB protein data bank (PDB ID:4EMA), was used as initial structure. The structure of Isoliquiritigenin was obtained from pubchem followed energy minimisation was done in gas phase using with the OPLS force eld [19,20,21] . Docking calculations were done with Glide and AutoDock 4.0 suite. Protein preparation wizard module of Schrodinger suite was for minimazation of protein and water molecules more than 3 Å away from the ligand was removed OPLSe force eld used for restrained molecular minimization to make relaxed structure wherein water molecules, heteroatoms were deleted and added hydrogen atoms. Site map module was used to characterize feature of binding site PPARγ for Isoliquiritigenin. Hydrophilic and hydrophobic maps were generated. Glide-receptor grid generation module was used for docking studies. Co-crystallised ligand was rst selected and around which grid was generated. Isoliquiritigenin and Troglitazone were then docked with PPARγ. Extra Precision (XP) was used to evaluate the binding modes of ligands. Docking study was done with default parameters. IFD employed with their default values, 20 poses for ligand were generated on each iteration. Thus, Inducedt docking (IFD) work ow used to generate other conformations of the PPARγ. Only the most favourable docking pose for each ligand was selected for structural analysis. IFD uses reduced van der waals radii, increased Coulomb-vdw cutoff and removes temporarily highly exible side chain of PPARγ and docked with different poses of ligand. Residues and ligands were minimized ranked according to Glide score.
AutoDock 4.0 suite was used for docking study, hydrogen atoms was added to protein structure, all nonpolar hydrogen atoms were merged, Lamarckian genetic algorithm was used as search parameter.
Binding energy calculation uses autodock scoring function such as short range van der Waals and electrostatic interactions, hydrogen bonding, entropy losses. Flexible ligand and a rigid receptor were used for docking [22,23,24,25,26,27] .
2.2 Molecular Dynamics simulation GROMACS 4.6.7 package was used for MD simulation of the complex using the GROMOS96 43a1 force eld [28]. The initial conformation for MD simulation was the conformation with most negative binding energy in docking study. The protein topology parameters were created by using the Gromacs program.
The topology parameters of ligand were built by the Dundee PRODRG server. The simulation was performed by keeping complex in a cubic box of simple point charge (SPC) water molecules. Minimization was carried out by the steepest descent method of 10000 steps then by the conjugate gradient method for 10000 steps. Position-restrain simulation was carried out at rst, system was simulated with a time step of 2 fs every time. The position-restrained dynamics simulation (NVT and NPT) at 300 K for 300 ps was done to equilibrate the system, nally run for 100000 ps. Electrostatic interactions in the system was estimated by PME algorithm. MD simulation was done for 100 ns [29,30,31,32] where γ is a coe cient related to surface tension of the solvent and b is tting parameter. ∆G sol, ∆E MM and -T∆S, indicates the changes of the solvation free energy, gas phase MM energy and the conformational entropy upon binding, respectively. ∆E MM calculate ∆E vdw (van der Waals) energies and ∆E internal (bond, angle, and dihedral energies) and ∆E electrostatic (electrostatic). ∆G solv is the sum of the nonelectrostatic solvation component (nonpolar contribution), ∆G SA electrostatic solvation energy (polar contribution). ∆G PB/GB , GB or PB model were used for polar contributions. Solvent accessible surface area (SASA) is used for nonpolar energy contribution. Conformational snapshots of MD simulations used for calculation of -T∆S. Principal component analysis (PCA) used to study overall dynamics of MD trajectory and to evaluate the covariance matrix of protein atom displacement from dominant and collective modes of the protein [33,34,35,36,37] .

Prediction of physiochemical properties
Molecules with functional groups have certain properties that are similar to known drugs. Hence, calculating the molecular property would be helpful in producing a good oral drug, and it is considered as important feature in drug discovery and development. Isoliquiritigenin was studied for ADME properties

Results And Discussion
PPARγ is very important drug target in type 2 DM. Thiazolidinediones (TZDs) have high a nity for PPAR and used as potent insulin sensitizers. The rst TZD introduced in early1997 was Troglitazone and subsequently Pioglitazone and Rosiglitazone was introduced in 1999. Our age long traditional practices utilizes Ptercarpus marsupium (local name-Vijaysar) to lower down blood sugar and in vivo experimental studies by others also validated its antidiabetic potential. Isoliquiritigenin is an important constituents of P. Marsupium and chosen for protein ligand interaction study with PPARγ using molecular docking and dynamic simulation studies. PPARγ protein was taken as target since it plays an important role in glucose metabolism.

Molecular Docking
SiteMap study analysed ve sites based on site score that considered various parameters such as amino acid exposure, hydrophobicity, hydrophilicity, donor/acceptor ratio, contact, size and volume. The amino acids in proposed binding active site having site score more than 1 [40] . PPARγ protein has ligand binding pocket in the centre of the ligand binding domain Pocket contains many polar residues such as Cys285, Ser289, His323, Tyr327, His449, and Tyr473. The central region of the ligand-binding pocket is surrounded mainly by nonpolar residues such as Leu330, Leu339, Leu353, and Met364. The Ω-pocket mainly contains hydrophobic residues Ile249, Met348, and Ile341, Leu255, Gly258, Ile262, Ile281 with few polar residues such as Glu259, Arg280, and Ser342 [41,42,43,44] .
Induced t docking was done to evaluate the conformational changes in PPARγ induced by Isoliquiritigenin binding. This method utilizes different poses, by using reduced van der waals radii and an increased coulomb-vdw cut off. Highly exible side chains were ignored while doing energy minimization to predict structure with different poses [45] . The best Induced t docking score of the optimised protein-ligand complexes was -9.39 Kcal/mol. The important residues in H bond were Cys285, Arg288, Tyr327 and Leu340. The van der waals interactions were observed with Gln286, Ser289, His323, Tyr327, Leu333, Val339, Ile341, Ser342, Phe363, His449 and pi-Alkyl interaction with Ile326 and leu330 ( Figure 3).
There were several classes of scaffolds developed for PPARγ agonists. The important ones are TZDs, Indoles and Benzimidazoles. It was observed that Tyr327, Arg288, and His323 plays critical role in activity of PPARγ. Some of the rst partial agonists developed for PPAR uses Indoles which interacts with residue Ser342 and van der waals interaction with Cys285 and Arg288. These compounds stabilize secondary structure by interacting with Ile341, also by hydrophobic contacts with Leu330 and Leu333.
Side chain of Ser289, Tyr327 and Tyr473 forms hydrogen bonds with Benzimidazoles. They also make hydrophobic contacts with Leu469 [46,47,48] . The present study as re ected in Figure 1a, 2a and 3 showed similar pattern of interaction between PPARγ and Isoiquiritigenin.

Molecular dynamics simulation (MD)
MD simulations utilizes theoretical models to study the spatial and energetic dynamics of PPARγ and Isoiquiritigenin at an atomic level. Binding free energy changes of a complex was examined to see the energetics and mechanisms of conformational change. The force between atoms is combination of bonded (Bond angle, bond length, Torsion) and non-bonded interactions (Coulomb & Leonard-jones).
During simulation, time is divided into discreet time steps e.g. 1 fs (10 −15 ) and the forces were calculated for each time steps while adjusting the position. Isoiquiritigenin collides with residues of PPARγ protein with preference for binding sites during 100 ns MD simulations and changes in structure and stability in water model was evaluated using GROMACS 4.6.7. This interaction in time dependent MD trajectories was recorded as RMSD, RMSF and radius of gyration. The conformation results obtained after simulation are more signi cant and stable than the docked conformation. Therefore, the binding orientation of Isoiquiritigenin with PPARγ predicted through MD simulation showed better correlation to their biological activity.
The structural variations in the PPARγ were analysed by root mean square deviation (RMSD) and the radius of gyration (Rg). RMSD of complex showed that after small rearrangement from the early conformation, the complex was fairly stable during complete MD simulation period. RMSD for the complexes of PPARγ with isoliquiritigenin showed that the structures of the systems equilibrated well after 6 ns of MD simulation (Fig. 4a). It showed that the RMSD pro les were always less than 0.4 nm for complex and PPARγ alone during the entire 100 ns simulation. Thus, the stability of PPARγ in Isoliquiritigenin bound state found as suitable candidate for post analysis. RMSD of Isoliquiritigenin showed that it stably bounded to PPARγ pocket. It showed stable pro le during the simulation (Fig. 4b).
The radius of gyration was calculated to analyse structural changes of PPARγ, when the Isoliquiritigenin was bound. The plot of radius of gyration in simulation time for PPARγ is given in gure 5. The average Rg values throughout the simulation time was 19.17 Å for the complex (Fig. 5) and 18.79 Å for PPARγ alone. Comparative analysis of nal pose of PPARγ -isoliquiritigenin complex after 10 ns molecular dynamics simulation with crystal structure of PPARγ revealed that binding of the Isoliquiritigenin did not bring signi cant conformational changes in the PPARγ structure and structure of protein was compact.  (Figure 6).
The H-bond between Isoliquiritigenin and PPARγ was analyzed and average of all H-bond candidates was calculated as 2.71 H-bonds as shown in Fig. 7.
The binding of the Isoliquiritigenin to PPARγ did not induce a large change in the solvent accessible surface area (SASA) of the complexes. SASA of PPARγ and PPARγ with Isoliquiritigenin was 139.53 and 143.12 respectively. The above results indicate the stability of the overall structure of protein when the inhibitor was bound.
Comparative analysis of different poses 100 ns simulation of PPARγ bound with Isoliquiritigenin was done. Two dimensional (2D) plots of nal pose were generated and compared with initial structure. The interaction of PPARγ with ligand Isoliquiritigenin was stable and remain in the same binding pocket. The important residues in H bond were Lys367, in van der waals interactions with Phe282, Gln 286, Ser 289, His 323, Tyr 327, Leu 333 and Phe 363 during entire simulation (Fig. 8a,b,c,d).

Principal Component Analysis
Protein's function is due to its conformational change in its three-dimensional structure. It is di cult to identify particular movement responsible for protein actions. PCA identify combined motions of atoms and recognise the structures underlying the atomic uctuations and it lter overall protein motions from local, fast domain motions. PCA offers matrix of eigenvectors and a diagonal matrix of eigenvalues. Each of the eigenvectors describes a direction and eigenvalue describe magnitude of motion. The axis which had maximum movement was the rst eigenvector. The second eigenvector was orthogonal to the rst. Plotting the projections against each other showed the path complex travelled. Principal components explain the direction in which movement is more spread out. Therefore, PC2 signi es the amount of variation that was not captured by PC1. PC1 and PC2 was denoted as the rst and second important conformation changes of the complex during the binding. The Fig. 9 and Fig. 10, showed the comparison of the PCA results for the PPAR-Isoliquiritigenin complex and the PPAR alone. It was evident from both the graph that dynamics of the protein from its initial conformation, remain unchanged with either ligand bound or independent. The overall uctuations in both simulations were well de ned by the rst ten eigenvectors, which account for PPAR alone and in complex as 72.46% and 79.60% of the total variance respectively. The PC1 & PC2 in Isoliquiritigenin bound and unbound form was different from each other. The graph was plotted between PC1 and PC2 of the covariance matrix of PPARγ alone showing more distribution in y axis than PPARγ bound with Isoliquiritigenin in x axis. There was difference in RMSF plot for PC1 and PC2 in Isoliquiritigenin bound and unbound form of protein( Fig. 11a  & b). The blue line indicates uctuations from the rst PC, and in orange line uctuations from the second principal component. The most uctuating atoms were His 217, Tyr 219, Phe 226, Arg 288, Ser 289, His 327 in PPARγ unbound state and bound state showed more at Asp 110, Ala 159. High uctuations of residues during the simulations, signifying a movement of these segments. However, the molecular dynamics of Isoliquiritigenin bound form of the protein did not have any identi able movement which suggest a stable conformation in the bound state.  [50] . Quantitative parameters such as RMSD showed structural stability of the protein ligand complex. Pro le of the protein and complex was found to be relatively stable about 0.34nm and 0.44nm similar which in acceptable range [51] . RMSF of given protein in MD trajectories showed uctuations more in N and C terminal ends whereas very low uctuations in the area where amino acid and ligand were interacting. Lower RMSD and small uctuations in RMSF and Rg of docking complex are good indication of system stability. MMPBSA was used to calculate free energies of Isoliquiritigenin, PPAR γ and complex. All the ensemble averages the change energy of bonded, nonbonded, polar and non-polar interactions for a binding a nity calculation, led to very stable binding pattern between protein and ligand under study as shown in Table 1 [11,34] . Drug likeliness property showed Isoliquiritigenin as an ideal candidate as it did not violate criteria given by Lipinsky. Similar reports on virtual screening and ADMET analysis showed by Liu et al., 2017 [49] .

Conclusion
The development of new drugs is a very complex and demanding interdisciplinary process. Virtual screening can be an inexpensive method to explore our plant resources to nd out new drug lead. Phytochemicals are rich reservoir of compounds which have protective property to our system. Pterocarpus marsupium (Vijaysar) have been used to lower down blood glucose in Indian traditional system, but there is no clear-cut understanding on its mechanism of action at molecular level. In the present study we performed Molecular Docking, Dynamics Simulation and Binding energy calculation on interaction of Isoliquiritigenin (an important constituent of Pterocarpus marsupium) with PPAR . Results indicated that PPAR had strong binding a nity for Isoliquiritigenin than the Troglitazone(marketed drug). Drug likeliness study also favours the drug like property of Isoliquiritigenin.

Declarations
Authors hereby declare that there is no con ict of interest.

Figure 1
Autodock results to show protein ligand interaction of PPARγ with Isoliquiritigenin and Troglitazone (a, b).

Figure 2
Glide docking studies on PPARγ with Isoliquiritigenin and Troglitazone conformation analysis to show protein ligand interaction   Free energy decomposition showing contribution of residue in terms of binding energy