Ligand and structure-based computational designing of multi-target molecules directing FFAR-1, FFAR-4 and PPAR-G as modulators of insulin receptor activity

Abstract Multi-agent therapies are an important treatment modality in many diseases based on the assumption that combining agents may result in increased therapeutic benefit by overcoming the mechanism of resistance and providing superior efficiency. Extensively validated 3D pharmacophore models for free fatty acid receptor-1 (FFAR-1), free fatty acid receptor-4 (FFAR-4), and peroxisome proliferator-activated receptor-G (PPAR-G) was developed. The pharmacophore model for FFAR-1 (r2 = 0.98, q 2 = 0.90) and PPAR-G (r2 = 0.89, q 2 = 0.88) suggested that one hydrogen bond acceptor, one hydrogen bond donor, three aromatic rings, and two hydrophobic groups arranged in 3D space are essential for the binding affinity of FFAR-1 and PPAR-G inhibitors. Similarly, the pharmacophore model for FFAR-4 (r2 = 0.92, q 2 = 0.87) suggested that the presence of a hydrogen bond acceptor, one negative atom, two aromatic rings, and three hydrophobic groups plays a vital role in the binding of an inhibitor of FFAR-4. These pharmacophore models allowed searches for novel FFAR-1, PPAR-G, and FFAR-4 triple inhibitors from multi-conformer 3D databases (Asinex). Finally, the twenty-five best hits were selected for molecular docking, to study the interaction of their complexes with all the proteins and final binding orientations of these molecules. After molecular docking, ten hits have been predicted to possess good binding affinity as per the Molecular Mechanics Generalized Born Surface Area (MM-GBSA) calculation for FFAR-1, FFAR-4, and PPAR-G which can be further investigated for its experimental in-vitro/in-vivo anti-diabetic activities. Communicated by Ramaswamy H. Sarma


Introduction
The peroxisome proliferator-activated receptors (PPARs) play a critical physiological role as lipid sensors and regulators of lipid metabolism (Berger & Moller, 2002). PPAR-G on distant tissues e.g. muscle and liver are likely to involve combined effects to enhance insulin-mediated adipose tissue uptake, storage of free fatty acids, induce the production of adiposederived factors with potential insulin-sensitizing activity and suppress the actions of insulin resistance causing adiposederived factors such as resistin (Berger & Moller, 2002). PPAR-G, a nuclear transcription factor, moderates the expression/ activity of G protein-coupled receptors (GPCRs) or free fatty acid receptors (FFARs), but its role in FFAR signaling is not clear (Yousefipour et al., 2009). FFARs have been proved to be good therapeutic targets for metabolic disorders (Bharate et al., 2009). The main reason for type 2 diabetes mellitus is the b-cell apoptosis (Butler et al., 2003) and decreased insulin secretion so, the best treatment of diabetes should be the b-cell restorative therapies and increased insulin secretion. FFAR-1 activation results in glucose-stimulated insulin secretion, b-cell growth, and further triggers the activation of PPAR-G (Bharate et al., 2009). FFAR-4 plays a key role in regulating the secretion of ghrelin (Dezaki et al., 2008), insulin, cholecystokinin (Lavine et al., 2010), incretin hormones like gastric inhibitory polypeptide (Iwasaki et al., 2015), and glucagon-like peptide-1 (Hirasawa et al., 2005).
One-third of almost currently marketed drugs assist as the target for the family of FFARs which provides the mechanism through which the cell receives the signals transmitted by extracellular factors (Bharate et al., 2009). In recent years, a number of pharmaceutical industries have been actively involved in the discovery of small molecule modulators of FFAR-1 and FFAR-4 receptors (Bharate et al., 2009). TAK-875, an FFAR-1 agonist is withdrawn from the market because of the severe hepatotoxicity induced by testing the drug taken for a longer duration (Kim et al., 2018). TUG-891, a potent agonist of FFAR-4 potential issues with ligand selectivity, response efficacy, and desensitization of the receptor (Hudson et al., 2013). KinDex Pharmaceuticals formed KDT-501, a double FFA4/PPARG agonist got from the hops plant for the treatment of T2D. However, there were no significant changes in oral glucose tolerance or insulin sensitivity measures (Kern et al., 2017). Thiazolidinediones (TZDs) comes under the class of insulin sensitizer that were used in the treatment of type 2 diabetes, but TZDs, the PPAR-G agonists have some severe side effects (Wright et al., 2014) including weight gain, fluid retention, congestive heart failure (American Diabetes Association, 2003) and bone fractures (Betteridge, 2011). Hence, there is a foremost need for safer and more effective PPAR-G directed therapeutics.
PPAR-G a widely used receptor in the treatment of T2D, however, the complete knowledge of signaling pathway of PPAR-G is still lacking (Ahmadian et al., 2013;Wang et al., 2015). The endogenous and exogenous agonists of PPAR-G bind and activate both PPAR-G and FFAR-1 which means it is an integrated two-receptor signaling pathway . PPAR-G requires FFAR-1 for the activation and formation of the complex which regulates the genes controlling metabolic activities. Potential advantages of FFAR-1 and FFAR-4 receptors as targets for the management of type 2 diabetes are; firstly, because FFAR-1 mediated insulin secretion is glucose-dependent, it is little or no risk of hyperglycemia, secondly, the glucose-dependent action might also result in better efficacy over other insulin secretagogues. Limited tissue distribution of FFAR-1 receptor suggests that there would be fewer chances of side effects and FFAR-4 mediated GLP-1 release might prevent a possible increase in body weight and b-cell failure which is often observed with insulin secretagogue therapy (Bharate et al., 2009) even it increases the b-cell replication (Gao et al., 2012;Leung & Zhang, 2014). Based on all the above points the FFARs and PPAR-G agonist should be an attractive approach for the management of type 2 diabetes. However, there is no single molecule that can bind on these targets and is based on a multi-target approach focusing on FFAR-1, FFAR-4, and PPAR-G modulating activities. Therefore, we aim to design these small molecules which can bind at PPAR-G, FFAR-1, and FFAR-4 providing all the benefits of these targets and could be a potential treatment in diabetes mellitus with minimal side-effects. In the present study, we have developed cross screened pharmacophore model which will have all the features for the FFAR-1, FFAR-4, and PPAR-G inhibitory activity. This generated pharmacophore model was validated and further utilized for the virtual screening to extract out new hits. These hits were further subjected to Lipinski's and ADMET filters followed by docking analysis in FFAR-1, FFAR-4, and PPAR-G crystal protein structures. The selected hits were subjected to MM-GBSA to access the binding affinity at the active sites of respective crystal structures (Scheme 1).

Computational details
All computational operations were run on Lenovo workstation having a 64-bit operating system, x-64 based Intel i3 core processor-4005 U CPU @1.70 GHz with 4GB RAM.
(nm) concentration to produce uniform expression of biological activity data for each dataset and were consequently employed for the development of the pharmacophore model. The ligands were sketched in ChemDraw2d and the preparation of the ligands was done using the Ligprep, version 2.5, User Manual, Schr€ odinger, LLC, New York (2012) applying OPLS_3E forcefield (LigPrep, 2012). The ionization states were generated of each inhibitor molecule at pH 7.0 ± 2.0.

Pharmacophore generation
The selected datasets having optimized molecules were imported into the PHASE program and subsequently distributed into active, inactive, and moderate active categories on the basis of their biological activity values. For the generation of a series of common pharmacophore hypotheses (CPH) the generated conformational space of each molecule was employed that were subsequently ranked by their various automatically calculated scoring parameters in the PHASE program. The generated models were evaluated on the basis of fitness score, matched sites, survival score, or survival minus inactive (S-I) score. The model with a preferentially high ability to pick active molecules has a high value of the survival score. The maximum fitness score is 3 and a high value of the S-I score means it is perfectly matching with the pharmacophore model and able to discriminate clearly between actives and inactives, which should be selected. Each molecule was assigned with pharmacophoric sites that are default in PHASE for the generation of pharmacophore hypotheses. PHASE includes six different types of sites as hydrogen bond acceptor (A), hydrogen bond donor (D), hydrophobic (H), negatively charged (N), positively charged (P), and aromatic ring (R) features which define a pharmacophoric site. Further, for the development or identification of common pharmacophore hypotheses, an active analog approach was implemented, in which CPHs were selected from the generated conformational space of active ligand molecules and were ranked by the criteria mentioned above. The generated hypotheses had very similar active and inactive statistical scores. Thus, to avoid the selection of similar hypotheses for analysis, all generated hypotheses were clustered and the representative hypothesis of each cluster was selected based on the highest value of fitness score, Phase Hyposcore, and S-I scores.

Validation of pharmacophore model
Validation was done to estimate the quality of the generated pharmacophore model by evaluating its ability to separate active compounds from inactives (Khedkar et al., 2007). Enrichment factor (EF), was calculated that measured the number of active hits obtained based on the recovery rate of actives against the decoy database. The goodness of hit score (GH) score was calculated to qualify model selectivity, the accuracy of hits, and the ability to segregate actives from a dataset consisting of known actives and inactives. The GH score generally ranges from 0 to 1 where 0 indicates it's a null model and 1 indicates it's an ideal model. EF and GH values were calculated using the below formulas (Khedkar et al., 2007): where Ht represented the number of hits obtained from the database, Ha indicated the number of active molecules in the hit list, D was a number of decoys and A was the number of actives. The best-generated pharmacophore models were also validated by using the calculation of specificity and sensitivity from the receiver operating curve (ROC). The pharmacophore models were subjected to testing the ability to correctly classify a list of compounds as actives or inactives was assessed by ROC analysis. The true positive rate (sensitivity) against the false positive rate (1-specificity) is plotted in a receiver operating characteristic (ROC) curve.

Virtual screening
The database of ligands was prepared by importing the molecules from commercial databases, namely Asinex 10452, Asinex synergy and elite, Asinex gold and platinum, and Asinex biodesign consisting of 10452, 94945, 203257, and 292926 compounds respectively. Database screening was done on the basis of validated pharmacophore models on which the rejection criteria are applied for those compounds which have an alignment score less than 1.2. All molecules which have a fitness value close to 3, matched sites more, Phase Hyposcore was high were selected. In the first step, individual models were screened through the database to obtain inhibitors against the respective targets. This is followed by screening the hits separately and alternatively with the other two models to achieve potential triple inhibitors. The triple inhibitors hence obtained were screened through the third pharmacophore model to achieve multi-targeted molecules. This led to the retrieval of hits that fit all of the features of the selected pharmacophore model of PPAR-G, FFAR-1, and FFAR-4 respectively. In order to find the novel molecules as potential anti-diabetic inhibitors, the database virtual screening was conducted summarized in Scheme 2.

Lipinski filter and QikProp properties
The screened hits were filtered on the basis of Lipinski rule of five and QikProp properties (Lipinski et al., 1997;QikProp, 2019). The duplicates of the same molecules were removed. After evaluating a large dataset (2,245 compounds) of marketed oral drugs, Lipinski et al. (1997) at Pfizer reached a series of conclusions on the optimum properties of drugs intended for oral administration (Lipinski et al., 1997). In order to achieve efficient oral absorption and cell permeability, drug candidates should follow the Lipinski rule.
QikProp predicts the widest variety of pharmaceutically relevant properties: octanol/water and water/gas log Ps, human oral absorption, blood/brain partition coefficient (log BB), aqueous solubility log S, overall CNS activity, Caco-2 and MDCK cell permeabilities, skin permeability (log Kp), number of likely metabolic reactions (metabol), log Khsa prediction for human serum albumin binding and log IC 50 for HERG K þchannel blockage, so that decisions about a molecule's suitability can be made based on a thorough analysis (QikProp, 2019).
After screening, molecules with good fitness scores, close to the scoring of active molecules in the set of pharmacophores, were subjected to structure-based drug designing techniques such as molecular docking.

Protein selection and preparation
For the determination of ligand-receptor interactions and filtering of virtual screened molecules, a glide docking program of Schr€ odinger software was employed. This docking program includes four docking levels, i.e. simple precision (SP), extra precision (XP), high throughput virtual screening (HTVS), and SP-peptide and for the present analysis, XP level was employed.
The crystal structures were selected on the basis of good resolution that is less than 2.0 Å. The crystal structures of each protein structure were downloaded and prepared using the 'protein preparation' program of Schr€ odinger software. The optimization of crystal structures includes deletion of solvent molecules (water), completion of the bond order, the addition of polar hydrogen atoms, correction of ionization and tautomeric states of amino acid residues, 3D framing of missing side chains and loops. Finally, all crystal complexes were energy minimized in 'impref' program using an OPLS_3E force field to RMSD value of 0.30 Å for the improvement of existing (any) steric clashes in the protein complexes due to the addition of external hydrogen atoms. The receptor grid was generated using the centroid of the co-Scheme 2. The filters incorporated and the molecules obtained after every filter in the protocol followed during the process of virtual screening of chemical databases.
crystallized ligands or mentioning any specific residue in the amino acid chain which has to be present in the docking (Jasuja et al., 2014).
Protein selection is done after cross-docking operation. Cross-docking is the process of extraction of all ligands from their crystal structures and then re-dock those ligands in each individual protein crystal (Jasuja et al., 2014). The extracted original crystal ligands were aligned over the redocked ligands and the deviation between two ligands was assessed as root mean square deviation (RMSD) (Jasuja et al., 2014). In order to signify the quality of crystal protein, the average RMSD of all the re-docked ligands in each protein was determined. A lower RMSD value indicated the ability of the crystal protein to dock the molecules more precisely (Jasuja et al., 2014). On the basis of protein selection, the crystal structures of PDB ID: 3P0G; resolution: Å, PDB ID: 4PHU; resolution: 2.33 Å and PDB ID: 1I7I, resolution: 2.35 Å were obtained from the protein databank (PDB). The selected crystal structures of respective proteins were then used for the docking analysis of study molecules as well as for the virtual screened molecules to explore the essential molecular interactions required for the binding and biological activity (Jasuja et al., 2014).

Docking analysis
To obtain information about the binding interactions between the macromolecules and the potential ligands, docking studies were performed. All the obtained hits were then docked in respective selected and prepared proteins for PPAR-G, FFAR-1, and FFAR-4. The results were analyzed for the presence of hydrogen bonding, hydrophobic, and p-p interactions between molecules and the active site of the proteins (Isogaya et al., 1999;Lu et al., 2010). The common interactions in all the complexes were scrutinized. Finally, the top-ranked poses were retained for Molecular Mechanics Generalized Born Surface Area (MM-GBSA) calculation (Scheme 3).

Binding affinity and MM-GBSA
The Prime MM-GBSA panel can be used to calculate ligandbinding energies and ligand strain energies for a set of ligands and a single receptor. The MM-GBSA score was used as the scoring function to evaluate the free energy of the binding of hits to the macromolecules. This methodology calculates the binding energy of ligand with the receptor as the sum of gas-phase internal energy, estimated using a molecular mechanics force field and the solvation free energy, calculated using an implicit solvent model (Generalized Born Surface Area). Additionally, in some cases, a change in entropy upon binding is also considered in the calculations to improve the accuracy of the binding affinity predictions. MM-GBSA method is significantly faster than FEP (free energy perturbation) calculations, however, it suffers from larger uncertainties than FEP, but in most cases is able to predict relative binding affinities in reasonable agreement with experimental data (Parenti & Rastelli, 2012).
MM-GBSA can be calculated by the formula: where MMGBSA dG Bind means the binding energy of Complex À Receptor -Ligand and NS means, no strain.

Pharmacophore generation
The generation of pharmacophore took place by dividing the molecules into active, inactive, and moderately active molecules. For FFAR-4 and FFAR-1 inhibitors, the threshold value of pIC 50 for the actives was set as more than or equal to 7 and inactive as less than or equal to 6, and for PPAR-G inhibitors threshold value of pIC 50 for the actives were set as more than or equal to 7 and inactive as less than or equal to 6. A total of 9 molecules were considered as actives and 6 molecules were considered as inactive and the rest of the molecules were considered as moderately active molecules for forming the pharmacophore model based on structures of active and inactive for FFAR-4 inhibitors. A total of 13 molecules were considered as actives and 7 molecules were considered as inactive and the rest of the molecules were considered as moderately active molecules for forming the pharmacophore model based on structures of active and inactive for FFAR-1 inhibitors. A total of 22 molecules were considered as actives and 9 molecules were considered as inactive and the rest of the molecules were considered as moderately active molecules for forming the pharmacophore model based on structures of active and inactive.
Based on common features among all the molecules and the best conformational pose the molecules were aligned. Under the hypothesis setting, the hypotheses generated from the pharmacophore model should match at least 50-70% of an existing dataset. The number of features was varied from four to seven for hypotheses generation. To get the maximum possible alignment and all maximum common features we need to generate the conformers.
On the basis of common features in total 731 hypotheses were generated GPR-40, 596 for GPR-120, and 1222 for PPAR-G inhibitors and simultaneously ordered according to their survival score and survival minus inactive score (S-I) that correspond to score active and score inactive respectively. The active and inactive scores so obtained in the hypotheses were very similar therefore, in order to avoid selection of similar kind of hypothesis, all the generated hypotheses were clustered and a representative of each cluster was selected on the basis of the highest S-I score. Therefore, the top five hypotheses each for GPR-40, GPR-120, and PPAR-G representing varied clusters were selected whose statistical parameters of these models are reported in Table 1.

3D QSAR studies with validation of the QSAR models
The 3D-QSAR models for each of the datasets were generated and divided the dataset molecules into test and training set molecules considering the uniform variation of biological activity of the molecules ( Table 2). The best hypothesis of FFAR-1 was ADHHRRR_1, indicating that FFAR-1 inhibitors have one hydrogen bond acceptor HBA (A), one hydrogen bond donor HBD (D), two hydrophobic groups (H), and three aromatic rings (R) features. This hypothesis was selected on the basis of the highest Q2 value i.e. 0.904 and also showed a high R2 value i.e. 0.985 with Fischer test value (F-value) 772.0. The large value of F indicated that the model is a statistically significant regression model. The Pearson-r value can range from À1 (the perfect negative linear relationship between variables) to 1 (the perfect positive linear relationship between variables). Similarly, best hypotheses of FFAR-4 (AHHHNRR_1) and PPAR-G (ADHHRRR_1) were selected which showed high Q2 values aS-I survival minus inactive score corresponds to score as inactives and high value of S-I represents the hypothesis that is more likely to pick active molecules than inactives. A: hydrogen bond acceptor; D: hydrogen bond donor; H: hydrophobic; R: aromatic ring; N: negative group. i.e. 0.873 (FFAR-4) and 0.886 (PPAR-G), high R2 0.919 (FFAR-4), 0.892 (PPAR-G) and F value 899.2 (FFAR-4) and 560.3 (PPAR-G). The AHHHNRR_1 comprised six features including one HBA (A), three hydrophobic rings (H), one negative atom (N), and two aromatic rings (R). The ADHHRRR_1 comprised seven features including one HBA (A), one HBD (D), two hydrophobic rings (H), and three aromatic rings (R). The selected models rendered good predictive power over other models of each of the targets respectively. The statistical results of generated QSAR models are mentioned in Table 3. The spatial arrangement of features along with their distance present in pharmacophore models of FFAR-1, FFAR-4, and PPAR-G and its mapping over their corresponding highest active molecules is shown in Figure 1 The correlation graphs obtained between the experimental and the predicted activity of training set molecules and test set molecules obtained from the best models are displayed in Figure 2   sensitive to omissions from the training set. The predicted activity of the training set and test set molecules for FFAR-1, FFAR-4, and PPAR-G is provided in Tables S1-S3 (Online Resource). In Y-randomization, the selected 3D QSAR models of FFAR-1, FFAR-4, and PPAR-G exhibited lower values of R 2 scramble i.e. 0.805, 0.777, and 0.693, respectively, as compared to corresponding original R 2 training values, confirming the trueness of the selected models. Also, in cross-validation of the selected 3D QSAR models of FFAR-1, FFAR-4 and PPAR-G exhibited lower values of R 2 CV i.e. 0.873, 0.536, and 0.832, respectively, as compared to corresponding original R 2 training values, confirming the trueness of the selected models.

Pharmacophore validation
The selected pharmacophore models were further validated using various techniques such as the calculation of EF, GH, and ROC. EF and GH scores are the matrices for qualifying the precision of hits and recall of actives from databases consisting of known actives and decoys. A significant model should have a GH score higher than 0.5. EF of all the models ranges from 2-3, indicating that these models have 2-3 times the capacity to pick actives than inactive. Details about the GH and EF scores are given in Table 4.
Further sensitivity and specificity parameters of the pharmacophore models were also calculated to determine the ROC of FFAR-1, FFAR-4, and PPAR-G. The ROC score and area under accumulation curve (AUAC) provides a practical way to measure the overall performance of the model as it helps us to understand the relation between model sensitivity (ability to discover true positives) and specificity (ability to avoid false positives). BEDROC score measures the early recognition enrichment of actives among the ranked compounds from the internal database. The ROC curve analysis of FFAR-1 yielded a ROC score of 0.92 ( Figure 3A), which means that in 9 out of 10 cases, a randomly selected active FFAR-1 inhibitor is ranked higher than an inactive one, the ROC curve analysis of FFAR-4 yielded ROC score as 1.00 ( Figure 3B), which means that in 10 out of 10 cases, a randomly selected active FFAR-4 inhibitors is ranked higher than an inactive one and the ROC curve analysis of PPAR-G yielded the ROC score 0.82 ( Figure 3C), which means that in 8 out of 10 cases, a randomly selected active PPAR-G inhibitors is ranked higher than an inactive one.

Virtual screening
Finally, the validated pharmacophore models were used for the database screening for FFAR-1, FFAR-4, and PPAR-G based on some filters like alignment score less than 1.2, fitness value high or close to 3, high Phase Hyposcore, more matched features with each pharmacophore models, and Lipinski rule with QikProp properties.
Nearly 40% of drug candidates fail in clinical trials due to poor absorption, distribution, metabolism, and excretion (ADME) properties. These late-stage failures significantly lead to the rapidly increasing cost of new drug development. Early detection of problematic candidates can histrionically lessen the extent of futile time and resources and simplify the overall development process (QikProp, 2019).
Accurate prediction of ADME properties prior to expensive experimental procedures, such as HTS, can eliminate unnecessary testing on compounds that will ultimately fail; ADME prediction was also used to emphasize lead optimization efforts for enhancing the desired properties of a given compound. Finally, incorporating ADME predictions as a part of the development process can generate lead compounds that are more likely to exhibit satisfactory ADME performances during clinical trials (QikProp, 2019).

Molecular docking
The hits obtained from the FFAR-1 were subsequently cross screened for FFAR-4 and PPAR-G and vice versa to obtain the triple inhibitors. To avoid the selection of false-negative molecules that did not interact properly with proteins and false-positive molecules that show the interactions but not favorable, docking analysis was carried out. The best molecules mapped over the FFAR-1, FFAR-4, and PPAR-G pharmacophore models are shown in Figure S1 (Online Resource).
On the basis of results obtained from cross-docking experiment crystal structure PDB ID: 4PHU (crystal structure of human FFAR-1 bound to allosteric agonist TAK-875), PDB ID: 3P0G (structure of a nanobody-stabilized active state of the beta2 adrenoceptor) and PDB ID: 1I7I (crystal structure of the ligand-binding domain of human PPAR-G in complex with the agonist AZ 242) for FFAR-1, FFAR-4, and PPAR-G were selected for further docking analysis of hits retrieved after the virtual screening. All the proteins were prepared as discussed in the experimental part and the whole process yielded a total of 10 potential triple inhibitors. Binding interactions of these molecules docked within the active sites of FFAR-1, FFAR-4, and PPAR-G followed important features identified in the three of the pharmacophore models.
After performing the cross-docking experiment, the active site of PDB 4PHU was determined using the co-crystallized ligand 2YB (allosteric agonist TAK-875) and found that the Arg183, Arg258, Tyr91, and Tyr240 are important sites for   interaction. The openness of the FFAR-1 binding pocket allowed the incorporation of polar groups such as spiro structure into the tail region. In contrast, the large and wide binding pockets of PPARs are usually relatively flat and enclosed inside the protein; thus, they can also accommodate polar groups. Major hit molecules show two hydrophilic residues, Arg183 and Arg258 that are electrostatically bonded and act as an anchor for the hydroxyl group of the head region of the ligand. Apart from the known interactions, in the hit molecules with high dock score, a strong hydrogen bond is also seen with the Ala 83, Leu 138, and Leu 135 in the hydrophilic pocket of the active site ( Figure 4). The protein structure 3P0G fits well with the configuration of the G-protein activation state due to the co-crystallized nanobody in the structure. High docking scores were achieved by some hit molecules when docked at the active site of 3P0G. These favorable high docked scores were received because of the reduction of the distance between TM3 and TM5 which arises upon activation of the active site proved from the incidence that the distance between the Ca carbons of Ser207 and Asp113 goes from 12.2 Å in the inactive structure to 11.4 Å in the activated 3P0G structure. The binding pocket between TM3 and TM6 is generally larger which allows a particularly prominent selectivity for agonists having bulky substituents to fit well in the active state of protein structure. The hit molecules occupied the same binding region as that of the co-crystallized nanobody found between TM III, TM V, TM VI, and TM VII regions. Through docking procedure, the important residues identified (Ser203, Ser207, Asp113, Lys305, Asn312, Tyr308, and Asp192) were indicated as very important for the agonistreceptor interactions. Most of the hit molecules show a common hydrogen bonding of the amido group with Asn 312, and Phe 193 and Phe 290 forms p-p interaction with the aromatic head regions. Generally, the linker of the hit molecules like triazole ring forms the hydrogen bond with Lys 305 which plays an important role in binding ( Figure 5).
Re-docking of the compound AZ 242 in the 1I7I crystal structure using the optimized docking parameters gave RMSD of 0.4209 Å respectively with a good dock score. Tyr473 is chiefly responsible for activation of PPAR-G, whereas His323, Ser289, and His449 are also important for the binding of the ligand in the ligand-binding domain. In 1I7I, the key residues His323, Tyr473, and Ser289 showed hydrogen bond with most compounds and His449 showed p-p interaction with the ligands (Figure 6).

MM-GBSA binding energy and glide score
A total of ten hits were retrieved that were observed to contain the above said essential interactions. The docked poses of these molecules were then subjected to MM-GBSA score calculation for calculating the binding energies of the respective complexes (Table 5).
All the hits were found to have docked and binding energy scores in the acceptable limits, i.e. the range of Glide XP G-Score for hits in complex with FFAR-1 was À7.820 to À9.293 and MM-GBSA binding energy was À51.9953 to À96.6919 kcal/mol. For the hits complexed with FFAR-4, the range of Glide XP G-Score was À5.814 to À7.551 and the MM-GBSA binding energy range was À wit2802 to À97.9551 kcal/mol. For the hits complexed with PPAR-G, the range of Glide XP G-Score was À10.066 to À12.092 and MM-GBSA binding energy range was À61.5427 to À83.2305 kcal/ mol. The overall Glide score and MM-GBSA binding energy of each hit molecule of the individual target are good and therefore it should proceed further to designing, synthesis, and biological evaluation.
The hits were also checked whether they have been earlier evaluated for FFAR-1, FFAR-4, and PPAR-G inhibitory activity using the SCI-Finder tool. The search results confined that the hits have not been previously tested experimentally for the inhibition of FFAR-1, FFAR-4, and PPAR-G. But as their synthesis scheme is of multiple steps and almost every hit molecule has a chiral center it would be tough to synthesis the required stereo-isomeric form rather than getting a racemic mixture which might not prove to be potent inhibitors. These results concluded that these hits can be used as lead molecules for designing the potential inhibitors of FFAR-1, FFAR-4, and PPAR-G.
The fatty-acid receptors FFAR-1 and FFAR-4 are implicated in the regulation of insulin secretion and insulin sensitivity, respectively. Although the properties of FFAR-4 and FFAR-1 converge on glucose homeostasis, few studies have addressed possible interactions between these two receptors in the control of b-cell function. Croze et al. (2020) found that FFAR-4 and FFAR-1 cooperate in an additive manner to regulate glucose homeostasis and insulin secretion (Croze et al., 2020). Furthermore, Sankoda et al. (2019) found that secretion of the gastric inhibitory polypeptide (a modulator of insulin resistance) is triggered by FFAR-4 and FFAR-1, without changing K-cell number or K-cell characteristics (Sankoda et al., 2019). Hence, a correlation amongst FFAR-1 and FFAR-4 has been established by these evidences.
The binding site of both FFAR1 and PPAR-G can be divided into three regions from top to bottom; polar hydrogen bonding region followed by two aromatic/hydrophobic regions. Acidic heterocylic groups such as indole, pyridine, quinalinone, thiazolidinedione, oxazolidinedione, tetrazole, and carboxylate can represent possible alternatives for binding with the critical polar residues within FFAR1 or PPAR-G binding sites (Helal et al., 2014). In addition, the distance between the two aromatic rings of the ligand should be carefully tailored to match the two aromatic/hydrophobic regions of the receptors (Helal et al., 2014). Based on the molecular binding interaction results, we assumed that a 4atom spacer is the most suitable. Also, branching at the terminal atom of the spacer, as in the case of rosiglitazone, is advantageous. Introduction of a mid-sized alkyl group on this atom can provide a better hydrophobic interaction with  the hydrophobic pocket which has a Y-shape in both FFAR1 and PPAR-G (Helal et al., 2014). These observations led us to propose a common pharmacophore for dual PPAR-G and FFAR1 agonists (Helal et al., 2014). This represents a valuable opportunity for the design of dual-acting agonists targeting both PPAR-G and FFAR1 for the treatment of diabetes (Figure 7). Some studies have suggested that, in some clinical cases, modulation of PPAR-G using TZDs is not sufficient for a full control of diabetes. For example, TZDs increased the insulin sensitivity in only two-thirds of patients at high risk of diabetes in a double-blind study (Troglitazone in Prevention of Diabetes \(TRIPOD\) study, 2004), a combination therapy containing a TZD (rosiglitazone) and a sulfonylurea (glimepiride) is now commercially available (Pf€ utzner et al., 2007). These findings suggest the need for co-administering insulin sensitizer and insulin secretagogue for the management of some diabetic cases. Design of dual modulators of PPAR-G and FFAR1 should lead to the discovery of a potent hypoglycemic agent that could act, for the first time, as insulin sensitizer as well as insulin secretagogue (Helal et al., 2014).
After working on thiazolidinedione hybrids, benzhydroland indole-based hybrids were synthesized and high intrinsic activities were noticed, with significantly larger maximum fold-induction compared to that of the PPAR-H full agonist reference, rosiglitazone (Darwish et al., 2018). The discovery of novel hybrid molecules possessing a balanced dual agonism on PPAR-G and FFAR1 receptors could provide less complex pharmacokinetic/pharmacodynamic relationships. Such profiles would lead to more predictable variability between patients and less extensive clinical studies (Morphy et al., 2004). Additionally, the afforded moderate activity of these hybrids, down to one digit lM, could be beneficial to limit the potential adverse effects raised with potent PPAR-G agonists. Nevertheless, such hybrids could keep an acceptable overall anti-diabetic efficacy as a result of their collective synergism on both receptors (Darwish et al., 2018).
The design of drugs with a dual mode of action is a valid approach. The concept of increasing insulin sensitivity and insulin secretion simultaneously has been successfully implemented using various marketed drug combinations as discussed above (Darwish et al., 2016(Darwish et al., , 2018Helal et al., 2014;Morphy et al., 2004;Pf€ utzner et al., 2007;Snitker et al., 2004). It is worth noting that the endogenous ligands of both PPAR-G and FFAR1 are fatty acids and considering a correlation of FFAR-1 with FFAR-4 and PPAR-G, targeting PPAR-G FFAR1 and FFAR-4 simultaneously should be a suitable approach towards the management of diabetes.

Conclusion
In recent years, many new targets for anti-diabetic drug development have been studied including PPAR-G, PTP-1B, FFAR-1, FFAR-4, DPP-4, and many others. In the current study, a structure-based pharmacophore modeling approach combined with other computational techniques like virtual screening, docking, and MM-GBSA studies were employed to identify novel FFAR-1, FFAR-4, and PPAR-G inhibitors. Molecular docking and binding affinity studies were carried out to study the binding modes and affinity of hits within the active pocket of FFAR-1, FFAR-4, and PPAR-G. Combined therapies in type 2 diabetes are of paramount importance which can preserve b-cell function, avoid adverse effects like weight gain, cardiovascular disorders, hepatotoxicity, hypoglycemia, decrease insulin resistance or increase insulin sensitivity depending on the glucose concentration and can provide a long duration of glycemic control, all these goals can be achieved from these hit molecules targeting FFAR-1, FFAR-4 and PPAR-G receptors. From the datasets of lakhs, we have narrowed down to ten compounds having triple receptor modulating activity, and based on the computational results these compounds can be used as leads for designing hit molecules which should be further synthesized and evaluated for biological activities. The top designed hits identified using combined pharmaco-informatic techniques may serve as potential agents for the treatment of diabetes.