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

doi:10.21203/rs.3.rs-1046506/v1

PDF

Research Article

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

https://doi.org/10.21203/rs.3.rs-1046506/v1

This work is licensed under a CC BY 4.0 License

Version 1

posted 05 Jan, 2022

You are reading this latest preprint version

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 fit 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 fulfils 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.

Drug Discovery, Design, & Development

Bioinformatics

Isoliquiritigenin

PPARγ

Free Energy

Molecular Dynamics Simulation

Binding energy

Diabetes mellitus (DM) is a type of metabolic disorder that causes hyperglycemia. DM is classified in two main type, type 1 and type 2. Type 1 DM is the consequence of complete or near-total insulin deficiency this might be due to interactions of genetic, environmental, or immunologic factors that lead to the destruction of the pancreatic beta cells and insulin deficiency. Type 2 DM can be classified 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. Antidiabetic agents like, ɑ-glucosidase inhibitors, thiazolidinediones (PPAR gamma agonist; decrease Insulin resistance, increase glucose utilization), Dipeptidyl peptidase IV inhibitors, GLP-1 receptor agonists, insulin secretagogues, biguanides are in clinical use. PPARs (peroxisome proliferator activated receptors) are members of the nuclear receptor superfamily, act as ligand inducible transcription factors. PPAR𝛼, PPAR 𝛿, and PPAR𝛾, are homologous subtypes each encoded by different genes and with different tissue expression and ligand selectivity. All PPARs have roles in fat, carbohydrate metabolism, cell growth and differentiation. It shows critical involvement in adipogenesis and for adipocyte formation. The binding of ligand induces conformational changes in PPAR𝛾. PPAR𝛾 has 5 domains (A, B, C, D, E) 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 fingers are involved in the recognition of specific 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 finding 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 flexibility. 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 flavonoid), pterosupin (a dihydrochalcone), marsupin (a benzofuranone), pterostilbene, liquiritigenin (a stilbene), isoliquiritigenin, (2S)-7-hydroxyflavanone, 7, 4’-dihydroxyflavone, 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.

2.2 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 field ^[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 field 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 first 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, Induced-fit docking (IFD) workflow 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 flexible 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 non-polar 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 field [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 first, 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, finally 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]}.

2.3 MM/GBSA & MM/PBSA calculations

Computational methods which involve molecular mechanics energy and solvation models, i.e. Molecular Mechanics/Poisson-Boltzmann Surface Area (MM/PBSA) and Molecular Mechanics/Generalized Born Surface Area (MM/GBSA) used for Free energy calculation:

∆G_bind= ∆H - T∆S = ∆E_MM + ∆G_sol - T∆S

∆E_MM = ∆E_internal + ∆E_{electrostatic} + ∆E_vdw

∆G_sol = ∆G_PB/GB + ∆G_SA

G _non−polar = γSASA + b

where γ is a coefficient related to surface tension of the solvent and b is fitting 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]}.

2.4 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 (Absortion, distribution, Metabolism and Excretion) by QuickProp module of Schrodinger Suite & Molinspiration server. Molecular property of the Isoliquiritigenin was predicted for both physiochemical and pharmacological properties such as LogP, hydrogen bonding characteristics, molecular size, and rotatable bonds. The combined effect of physicochemical properties, pharmacokinetics and pharmacodynamics, and the drug-likeness model score were identified with Molsoft server. The scores were further compared with the standard anti-diabetic drug. Lipinski's rule of five was used to evaluate the acceptability of the Isoliquiritigenin ^{[38, 39]}.

PPARγ is very important drug target in type 2 DM. Thiazolidinediones (TZDs) have high affinity for PPAR𝛾 and used as potent insulin sensitizers. The first 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.

3.1 Molecular Docking

SiteMap study analysed five 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.0 and residues identified as 226, 227, 228, 247, 249, 255, 256, 258,259, 260, 280, 282, 288, 291, 292, 295, 296, 323, 326, 327, 329, 330, 333, 340, 341, 342, 343, 344, 345, 346, 348, 363, 364, 449, 469, 473 and 601 ^[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]}.

Isoliquiritigenin and Troglitazone (known agonist) docked into a PPARγ target structure in a different positions, conformations and orientations and ligand-protein pairwise interaction energies were calculated. The best possible binding modes of the Isoliquiritigenin and Troglitazone, at targeted protein’s active sites are displayed in Fig. 1a& b. In our study, H-bonds were formed between Tyr327, and Tyr473 of the PPARγ and Isoliquiritigenin with binding energy of -7.46 kcal/mol. The van der waals interaction were seen with Phe282, Gln286, Ser289, His323, Leu333, Phe363, Met364, Lys367, His449, Glu453, Leu465 and Leu469. Similarily, the pi-Sigma interaction with Ile326, pi-Alkyl with Ala292, leu330 and Pi-sulphur interaction with Cys285. Molecular docking results with Troglitazone showed H bond with Tyr327, Ser289 with binding energy of -11.01 kcal/mol.

The van der waals interactions was seen with Phe282, Gly284, Gln286, His323, Leu333, Phe363, Lys367, His449, Glu453, Leu465, Leu469 and Tyr473. Similarly, the pi-Alkyl with Cys285, Arg288, Met364, Leu330, Leu333, Val339, Ile341 and pi-sulphur interaction with His323. Both the ligands Isoliquiritigenin and Troglitazone were further evaluated with another docking software Glide with Extra precision (XP) mode to find out accurate binding during PPAR interaction with different poses of ligand. In case of Isoliquiritigenin, the binding energy in XP was -6.74 kcal/mol. The important residues in H bond were Leu340 and Tyr473, in van der waals interactions was with Phe282, Gln286, Ser289, His323, Tyr327, Leu333, Phe363, Lys367, His449, Glu453, Leu465 and Leu469 (Figure 2a). The pi-Alkyl interaction with Arg288, Leu330, Val339 and pi-sulphur interaction with Met364 (Figure 2a). Troglitazone docking, gave binding energy in XP mode as -9.59 kcal/mol. The important residues in H bond were His323 and Ser289, while Phe282, Gly284, Gln286, Leu333, Phe363, Lys367, His449, Glu453, Leu465, Leu469, and Tyr473 showed van der waals interactions. The pi-Alkyl interaction with Cys285, Arg288, Leu330, Leu333, Val339, Ile341and pi-sulphur interaction with Met364 (Figure 2b).

Induced fit 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 flexible side chains were ignored while doing energy minimization to predict structure with different poses ^[45]. The best Induced fit 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 first 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 reflected in Figure 1a, 2a and 3 showed similar pattern of interaction between PPARγ and Isoiquiritigenin.

3.2 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⁻¹⁵) 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 significant 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 profiles 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 profile 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 figure 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 final 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 significant conformational changes in the PPARγ structure and structure of protein was compact.

Root mean square fluctuation (RMSF) of PPAR γ around its average conformations play an important indicator of protein activity in complex formations. Significant change in complex fluctuation occur around residues 240, 243, 253, 257,477 of N and C terminal respectively. Protein residue 357, 443, 447, 448, 457, 463 got stabilized on ligand binding. Residues critical in interaction with ligand were Arg 288, Tyr 327, Ser 289, His 323, Leu333, His449 and Tyr 473 (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 final 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).

3.3 Principal Component Analysis

Protein's function is due to its conformational change in its three-dimensional structure. It is difficult to identify particular movement responsible for protein actions. PCA identify combined motions of atoms and recognise the structures underlying the atomic fluctuations and it filter 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 first eigenvector. The second eigenvector was orthogonal to the first. 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 signifies the amount of variation that was not captured by PC1. PC1 and PC2 was denoted as the first 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 fluctuations in both simulations were well defined by the first 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 fluctuations from the first PC, and in orange line fluctuations from the second principal component. The most fluctuating 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 fluctuations 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 identifiable movement which suggest a stable conformation in the bound state.

Table 1

MMPBSA Free Energy Analysis of PPARγ – Isoliquiritigenin Complex After 100ns MD Simulation
	Vander waal (KJ/mol)	Electrostatic (KJ/mol)	Polar Solvation (KJ/mol)	SASA (KJ/mol)	Binding energy (KJ/mol)
Complex	-179.48±35.20	-38.97±11.65	107.92±28.11	-16.14±3.17	-126.67±24.92

Figure 12 displayed energy decompositions of all residues in the complex. It was calculated to qualitatively find out the key residue that contributed significantly in the Ligand binding. The major contributing residues were found to be 69, 106, 109, 142, 146, 248.

3.4 Drug-likeness of Isoliquiritigenin

Isoliquiritigenin was studied to predict its drug-likeness properties by considering following molecular descriptors: logP (partition coefficient), molecular weight (MW), topological polar surface area (TPSA), hydrogen bond acceptors and donors count in a molecule. Isoliquiritigenin was found to have logP as 2.77, TPSA as 77.75, total number of atoms (n Atoms) : 19, Mol .Wt : 256.26, number of hydrogen bond acceptors was 4, number of hydrogen bond donors: 3, and 0 violation of the Lipinsky’s rule. Lipinski, Ghose and Veber rules states that membrane permeability must possess logP ≤ 5, number of hydrogen bond acceptors ≤10, number of hydrogen bond donors ≤ 5, molecular weight ≤500, topological polar surface area (TPSA) < 140 Ǻ and number of rotatable bonds (n rotb) < 10 (measures molecular flexibility) and also, the total number of atoms between (n Atoms) 20 and 70. This was used as filter for drug-like properties ^{[38,39,40,41,42,43,44,45,46,47,48, 49]}.

PPRs acts as ligand inducible transcription factors belongs to members of the nuclear receptor superfamily. They are mainly involved in energy, carbohydrate and lipid metabolism. Lee et al., (2017) showed binding modes of Pioglitazones and Lobeglitazones with PPAR γ. The head groups of TZDs forms H bonds with Tyr 473 of PPAR γ in active conformations, which correlates with full agonism of the drugs. Similar binding mode was observed in molecular docking results ^[42]. The pi interactions with Phe 264, Phe 363 and h bond with Leu 340, Tyr 473. in silico identification of PPAR γ agonist from Chinese medicine showed similar interactions as H bond with residue Tyr 327, Lys 367 etc ^[50]. Quantitative parameters such as RMSD showed structural stability of the protein ligand complex. Profile 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 fluctuations more in N and C terminal ends whereas very low fluctuations in the area where amino acid and ligand were interacting. Lower RMSD and small fluctuations 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 affinity 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].

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 find 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 affinity for Isoliquiritigenin than the Troglitazone(marketed drug). Drug likeliness study also favours the drug like property of Isoliquiritigenin.

Authors hereby declare that there is no conflict of interest.

Famm, K., Litt, B., Tracey, K. J., Boyden, E. S. & Slaoui, M. Drug discovery: A jump-start for electroceuticals. Nature496, 159–161. https://doi.org/10.1038/496159a (2013).
Tontonoz, P. , Spiegelman, B.M. Fat and beyond: the diverse biology of PPAR𝛾, Annual Review of Biochemistry, 77, 289–312 (2008).
Zoete, V., Grosdidier, A., Michielin, O. Peroxisome proliferator-activated receptor structures: ligand specificity, molecular switch and interactions with regulators, Biochimica et Biophysica Acta—Molecular and Cell Biology of Lipids, 1771,8, 915–925 (2007).
Chandra, V., Huang, P., Hamuro, Y. et al. Structure of the intact PPAR-𝛾-RXR-𝛼 nuclear receptor complex on DNA,” Nature, 456,220, 350–356 (2008).
Hummasti, S., and Tontonoz, P. The peroxisome proliferator activated receptor N-terminal domain controls isotype-selective gene expression and adipogenesis,” Molecular Endocrinology, 20 (6) 1261–1275 (2006).
Bugge, A. and Mandrup, S. Molecular mechanisms and genome-wide aspects of PPAR subtype specific transactivation, PPAR Research, ID169506,12 (2010).
Helsen, C. and Claessens, F. Looking at nuclear receptors from a new angle, Molecular and Cellular Endocrinology, 382, 1, 97–106 (2014).
Jpenberg, A.I., Jeannin, E., Wahli, W. and Desvergne, B. Polarity and specific sequence requirements of peroxisome proliferator activated receptor (PPAR)/retinoid X receptor heterodimer binding to DNA. A functional analysis of the malic enzyme gene PPAR response element, The Journal of Biological Chemistry , 272 (32) 20108–20117(2010).
Gohlke, H. and Klebe, G. Approaches to the description and prediction of the binding affinity of small-molecule ligands to macromolecular receptors. Ang Chem Int Ed,. 41, 2644-76 (2002).
Shirts, M.R., Mobley, D.L., Chodera, J.D. Alchemical free energy calculations: ready for prime time. Ann Rep Comput Chem, 3, 41-59 (2007).
Homeyer, N., Stoll, F., Hillisch, A. et al. Binding free energy calculations for lead optimization: assessment of their accuracy in an industrial drug design context. J Chem Theory Comput. 10, 3331-44 (2014).
Kollman, P. A., Massova, I., Reyes, C., Kuhn, B., Lee, S. H. L., Duan, T., Y., Wang, W., Donini, O., Cieplak, P., Srinivasan, J., Case, D. A., Cheatham, T. E. Calculating structures and free energies of complex molecules: Combining molecular mechanics and continuum models. Acc Chem Res, 33 (12), 889–897 (2003).
Aggarwal, B.B., Prasad, S., Reuter, S., Kannappan, R. , Yadev, V.R., Park, B. , Kim, J.H., Gupta, S.C., Phromnoi, K., Sundaram, C., Prasad, S., Chaturvedi, M. M., Sung, B. Identification of novel anti-inflammatory agents from ayurvedic medicine for prevention of chronic diseases. Curr Drug Targets. 1, 12(11), 1595–1653 (2011).
Kameswara, B., Rao, G., Kesavulu, R., M.M. and Apparao, C. Effect of oral administration of bark extracts of Pterocarpus santalinus L. on blood glucose level in experimental animals. Journal of Ethnopharmacology, 74,1, 69–74 (2001).
Chakravarthy, B.K., Gupta, S., Gambhir, S.S. and Gode, K.D. The prophylactic action of (-)-epicatechin against alloxan induced diabetes in rats, Life Sciences. 29 (20) 2043–2047(1981).
Maurya, R., Ray, A.B. Constituents of Pterocarpus marsupium. J Nat Prod. 47,(1), 179-81 (1984).
Tripathi, J. , Joshi, T. Flavonoids from Pterocarpus marsupium. Planta Med. 54 (4), 371-372 (1988).
Halagappa, K. , Girish, H.N., Srinivasan, B.P. The study of aqueous extract of Pterocarpus marsupium Roxb. on cytokine TNF-α in type 2 diabetic rats. Indian J Pharmacol. 42(6), 392–396 (2010).
Yadav, V.K. & Mishra, A. In vitro & in silico study of hypoglycemic potential of Pterocarpus marsupium heartwood extract, Natural Product Research, 33(22), (2019) 3298-3302.
Glide, version 6.7, Schrödinger, LLC, New York, NY, (2015).
LigPrep, version 3.4, Schrödinger, LLC, New York, NY, (2015).
MacroModel, version 10.8, Schrödinger, LLC, New York, NY, (2015).
Singh, A. & Mishra, A. Molecular dynamics simulation and free energy calculation studies of Coagulin L as dipeptidyl peptidase-4 inhibitor, Journal of Biomolecular Structure and Dynamics, DOI: 10.1080/07391102.2020.1822917 (2020).
Sastry, G.M., Adzhigirey, M., Day, T. et al. Protein and ligand preparation: parameters, protocols, and influence on virtual screening enrichments. J Comput Aided Mol Des. 27, 221–234. doi: 10.1007/s10822-013-9644-8 (2013).
Banks, J.L. , Beard, H.S., Cao, Y., Cho, A.E., Damm, W., Farid, R. , Felts, A.K., Halgren, T.A., Mainz, D.T., Maple, J.R. R. Murphy, D.M. Philipp, M.P. Repasky, L.Y. Zhang, Berne, B.J., Friesner, R.A., Gallicchio, E., Levy, R.M. Integrated Modeling Program, Applied Chemical Theory (IMPACT). J Comp Chem. 26, 1752–1780. doi: 10.1002/jcc.20292 (2005).
Morris, G.M., Huey, R., Lindstrom, W., Sanner, M.F., Belew, R.K., Goodsell, D.S. et al. AutoDock4 and AutoDockTools4: Automated docking with selective receptor flexibility. J Comput Chem, 30, 2785–91 (2009).
M.A. Lill, M.L. Danielson, Computer-aided drug design platform using PyMOL. J Comput Aided Mol Des, 25, (2011) 13–9.
Wallace, A.C., Laskowski, R.A., Thornton, J.M. LIGPLOT: A program to generate schematic diagrams of protein-ligand interactions. Protein Eng. 8, 127–34 (1995).
Berendsen, H.J.C., van der Spoel, D., Drunen, R.van. GROMACS: A message-passing parallel molecular dynamics implementation. Computer Physics Communications, 91, 1–3, 2, 43-56 (1995).
van Gunsteren, W.F., Billeter, S.R., Eising, A.A., Hunenberger, P.H., Kruger, P.K.H.C., Mark, A.E., Scott, W.R.P., Tironi, I.G. Biomolecular Simulation: The GROMOS96 Manual and User Guide, VdfHochschulverlag AG, Zürich. (1996).
Pronk, S. et al. GROMACS 4.5: a high-throughput and highly parallel open source molecular simulation toolkit. Bioinformatics, (2013).
Schuttelkopf, A.W., van Aalten, D.M.F. PRODRG: a tool for high-throughput crystallography of protein-ligand complexes. Acta Crystallogr. 60, 1355–1363 (2004).
Verma,S., Singh. A., Mishra, A. Molecular Dynamics Investigation on the Inhibition of MDM2-p53 Interaction by Polyphenols, Molecular informatics, 32, 2, 203-212 (2013).
van Alten, D.M., Findlay, J.B., Amadei, A., Berendsen, H.J. Essential dynamics of the cellular retinol-binding protein—evidence for ligand-induced conformational changes. Protein Eng, 8, 1129–1131 (1995). PMID:8819978.
Wang, J. M., Hou, T. J., Xu, X. J. Recent advances in free energy calculations with a combination of molecular mechanics and continuum models. Curr. Comput.Aided Drug Des, 2 (3), 287–306 (2006).
Amadei, A., Linssen, A.B., Berendsen, H.J. Essential dynamics of proteins, Proteins. 17, 412–425 (1993). PMID:8108382 33.
Kumari, R., Kumar, Open Source Drug Discovery Consortium ,LynnA.g_mmpbsa:A GROMACS Tool for High-Throughput MMPBSA Calculations. J Chem Inf Model. 54, 1951–1962 (2014). doi:10.1021/ ci500020mPMID:24850022.
Munnaluri, R., Sivan, S.K. and Manga, V. Molecular docking and MM/GBSA integr.ated protocol for designing small molecule inhibitors against HIV-1 gp41 Med.Chem. Res. 24, 829 (2015).
Lipinski, C.A., Lombardo, F., Dominy, B.W., Feeney, P.J. Experimental and computational approaches to estimate solubility and permeability in drug discovery and development settings. Adv Drug Deliv Rev. 46, 3–26 (2001).
Jarrahpour, A., Fathi, J., Mimouni, M., Hadda, T.B., Sheikh, J., Chohan, Z. et al.,. Petra, Osiris and Molinspiration (POM) together as a successful support in drug design: Antibacterial activity and biopharmaceutical characterization of some azo schiff bases. Med Chem Res. 21, 1–7 (2012).
Halgren, T. Identifying and Characterizing Binding Sites and Assessing Druggability, J. Chem. Inf. Model., 49, 377–389 (2009).
Kroker, A.J. and John Bruning, B. Review of the Structural and Dynamic Mechanisms of PPARγ Partial Agonism. PPAR Research, ID 816856, doi.org/10.1155/2015/816856 (2015).
Lee, M.A., Tan, L., Yang, H., Im, Y.G. , Im Y.J. Structure of PPRgamma complexed with lobeglitazone and pioglitazone reveal key determinants for the recognition of antidiabetic drugs. Scientific reports, Nature. 7, 16837, (2017) DOI:10.1038/s41598-017-17082.
Amano, Y., Yamaguchi, T., Ohno, K., Niimi, T., Orita, M., Sakashita, M. Hypertens. Res. 35, 715-719 (2012). Lu, I.L., Huang, C.F., Peng, Y.H., Lin, Y.T., Hsieh, H.P., Chen, C.T. et al. Structure-based drug design of a novel family of PPAR partial agonists: virtual screening, X-ray crystallography, and in vitro/in vivo biological activities, J. Med. Chem. 49, 2703-2712 (2006).
Prime, version 4.0, Schrödinger, LLC, New York, NY, (2015).
Liberato, M.V., Nascimento, A. S., Ayers, S.D. et al. Medium chain fatty acids are selective peroxisome proliferator activated receptor (PPAR) 𝛾 activators and Pan-PPAR partial agonists, PLOS ONE, 7,(5) IDe36297 (2012).
Einstein, M., Akiyama, T.E., Castriota, G.A. et al., The differential interactions of peroxisome proliferator-activated receptor𝛾 ligands withTyr473 is a physical basis for their unique biological activities. Molecular Pharmacology, 73 (1) 62–74(2008).
M Sime, Allan, A.C., Chapman, P. et al. Discovery of GSK1997132B a novel centrally penetrant benzimidazole PPARγ amma partial agonist. Bioorganic and Medicinal Chemistry Letters , 21(18), 5568–5572 (2011).
Liu, X., Jing, Z., Jia, W.Q., Wang, S.Q., Ma, Y., Xu, W.R., Liu, J.W., Cheng, X.C. Identification of novel PPAR alpha gamma dual agonists by virtual screening, ADMET prediction and MD simulations. JBSD. DOI: 10.1080/07391102.2017.1373706 (2017).
Chen, K.C., Chen, C.Y.C. in silico identification of potent PPAr agonist from traditional medicine a bioactivity prediction virtual screening and MD simulation study; PLO S, 192452, DOI; 10.1155/2014/192452(2014).
Katari, S. K., Natarajan, P., Swargam, S., Kanipakam, H., Pasala, C., Umamaheswari, A. Inhibitor design against JNK1 through e-pharmacophore modeling docking and molecular dynamics simulations, https://doi.org/10.3109/10799893.2016.1141955 (2016).
Kira, A., Armacost, S. Riniker, and Cournia, Z. Exploring Novel Directions in Free Energy Calculations, Journal of Chemical Information and Modeling 60 (11), 5283-5286 DOI: 10.1021/acs.jcim.0c01266 (2020).

No competing interests reported.

PDF

Version 1

posted 05 Jan, 2022

You are reading this latest preprint version

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

Status:

Version 1

Abstract

Figures

1. Introduction

2. Methods

2.2 Molecular Docking

2.2 Molecular Dynamics simulation

2.3 MM/GBSA & MM/PBSA calculations

2.4 Prediction of physiochemical properties

3. Results And Discussion

3.1 Molecular Docking

3.2 Molecular dynamics simulation (MD)

3.3 Principal Component Analysis

3.4 Drug-likeness of Isoliquiritigenin

4. Conclusion

Declarations

References

Additional Declarations

Status:

Version 1