key: cord-0023921-bmj3nf9b authors: Pal, Ajay; Curtin, James F.; Kinsella, Gemma K. title: Structure based prediction of a novel GPR120 antagonist based on pharmacophore screening and molecular dynamics simulations date: 2021-11-05 journal: Comput Struct Biotechnol J DOI: 10.1016/j.csbj.2021.11.005 sha: a4a075be387751787349bb201c7aab9f23c73033 doc_id: 23921 cord_uid: bmj3nf9b The G-protein coupled receptor, GPR120, has ubiquitous expression and multifaceted roles in modulating metabolic and anti-inflammatory processes. Recent implications of its role in cancer progression have presented GPR120 as an attractive oncogenic drug target. GPR120 gene knockdown in breast cancer studies revealed a role of GPR120-induced chemoresistance in epirubicin and cisplatin-induced DNA damage in tumour cells. Higher expression and activation levels of GPR120 is also reported to promote tumour angiogenesis and cell migration in colorectal cancer. Some agonists targeting GPR120 have been reported, such as TUG891 and Compound39, but to date development of small-molecule inhibitors of GPR120 is limited. Herein, following homology modelling of the receptor a pharmacophore hypothesis was derived from 300 ns all-atomic molecular dynamics (MD) simulations on apo, TUG891-bound and Compound39-bound GPR120S (short isoform) receptor models embedded in a water solvated lipid bilayer system. We performed comparative MD analysis on protein–ligand interactions between the two agonist and apo simulations on the stability of the “ionic lock” – a Class A GPCRs characteristic of receptor activation and inactivation. The detailed analysis predicted that ligand interactions with W277 and N313 are critical to conserve the “ionic-lock” conformation (R136 of Helix 3) and prevent GPR120S receptor activation. The results led to generation of a W277 and N313 focused pharmacophore hypothesis and the screening of the ZINC15 database using ZINCPharmer through the structure-based pharmacophore. 100 ns all-atomic molecular dynamics (MD) simulations were performed on 9 small molecules identified and Cpd 9, (2-hydroxy-N-{4-[(6-hydroxy-2-methylpyrimidin-4-yl) amino] phenyl} benzamide) was predicted to be a small-molecule GPR120S antagonist. The conformational results from the collective all-atomic MD analysis provided structural information for further identification and optimisation of novel druggable inhibitors of GPR120S using this rational design approach, which could have future potential for anti-cancer drug development studies. G protein coupled receptors (GPCRs) are potential drug targets, with more than 40% of approved drugs modulating the action of one or more GPCRs [1] . Through their interaction with a wide number of bioactive molecules, including ions, lipids, amino acids, peptides, proteins, and small organic molecules these receptors play a crucial role in many physiological processes. These processes include the transmission of light and of odorant signals, mediation of neurotransmission and hormonal action, cell growth and the immune defence [1] [2] [3] . Among GPCRs, GPR120, a member of the rhodopsin sub-family was deorphanized in 2005 [4] as a Free fatty acid receptor (FFAR), which has strong affinity for long-chain fatty acids (LCFAs). GPR120 (or free fatty acid 4 receptor, FFAR4) has recently been explored for its role in cancer cell proliferation and tumour angiogenesis through the binding of Omega-3 polyunsaturated fatty acids (x-3 PUFAs), mainly docosahexaenoic acid (DHA) and eicosapentaenoic acid (EPA) [5, 6] . In humans, GPR120 is present in two isoforms: a short isoform at 361 amino acids -(GPR120S) and a long isoform at 377 amino acids (GPR120L). GPR120L contains an additional 16 mammalian species [4, [7] [8] [9] [10] . In humans, GPR120S couples effectively to both Gq/11 and arrestin pathways whereas GPR120L is an arrestin-biased receptor [7, 10, 11] . The functional significance of the two isoforms in humans remains unclear as there is no tissue reported in which GPR120L is expressed selectively and both isoforms have identical endogenous substrate binding sites known as orthosteric binding pockets [8] . Both the short and long isoforms of GPR120 share the same endogenous substrates which bind at the orthosteric pocket and activate the downstream signalling through second messengers [8, 9] . Besides the orthosteric pocket, GPCRs also contain remote binding sites on the receptor surface known as allosteric sites. These allosteric sites can indirectly modulate the receptor signalling by changing the conformation of orthosteric binding pocket upon ligand binding at allosteric sites [9] . When GPR120S is activated by x-3 PUFAs (orthosteric endogenous ligands), it couples with Gq/11 and induces mobilisation of intracellular calcium [12] while activation by LCFAs and PUFAs also results in b-arrestin2 recruitment [13] . The basal phosphorylation of GPR120S is two-fold higher than that of GPR120L in the absence of endogenous ligands [4, 14] . GPR120 is a potent lipid mediator in metabolic diseases such as obesity and type-2 diabetes mellitus (T2DM) [15] . The antiinflammatory actions of GPR120 are induced through b-arrestin2/TAB1 interactions, independent of G protein -Gq/G11coupled pathway [15] . Recent studies have developed potent and selective GPR120 agonists such as GW9508 -a partial agonist of GPR120 [6] , TUG-891 -a selective potent agonist of GPR120 [12] and compound 39 (a benzofuran propanoic acid analogue) from Merck & Co., Inc [16] , but no orthosteric GPR120 antagonists are available to date (Fig. S1 ). However, AH7614 (Fig. S1 ) has been shown to act as a negative allosteric modulator of GPR120. [17, 18] . GPR120 has also been shown to regulate tumour growth and migration of various cancer types, including melanoma and prostate cancers [5] . GPR120 activation results in a signalling cascade that releases lysophosphatidic acid species via activation of cytosolic phospholipase A2 to induce resistance to cisplatin-induced DNA damage in tumour cells. Such inducedchemoresistance was limited in GPR120 'knockout' studies [19] . A similar study in breast cancer cells reported high expression of GPR120 in breast cancer cells upregulated the ABC transporters expressions which induce resistance to epirubicin-induced cell death. [20] In contrast, GPR120 knockdown in breast cancer cells reduced the epirubicin resistance induced by GPR120 agonist treatment. [20] The Human Protein Atlas (HPA) analysed RNA-seq data from 64 human cell lines and 37 human tissues samples to estimate protein expression of GPR120 [20] , with the highest TPM (transcripts per million) of GPR120 being in rectum and colon tissues [5, 6, 21] . Wu and co-workers screened a panel of malignant and nonmalignant human colorectal cell lines and reported higher expression of GPR120 in all studied colorectal cancer (CRC) cell lines compared to non-CRC cell lines [6] . GW9508 tested in CRC cell lines demonstrated activation of PI3K/Akt and nuclear factor kB [6] . GW9508 has been reported to be mitogenic for CRC cells [5, 6] , but not for prostate cancer cells [22] . Another report by Liu et al [22] showed that the GPR120 agonist EPA activates ribosomal protein S6 kinase b-1 (p70S6K1) in the CRC CaCo2 cell line. GW9508 (at 10 mM dose) stimulates migration of CRC cells whereas TUG-891 (a GPR120 selective agonist) inhibits migration of prostate cancer cells [6, 22] . These studies suggested that GPR120 activation (rather than antagonism as suggested by Wu et al 2013 [6] ) might be effective in anticancer therapy. These two results (on CRC cells and prostate cancer cells) stand in contrast and highlight the importance of exploring the mechanisms underlying the opposing effects of GPR120 agonists on proliferation in different types of cancer cells. The expansive expression pattern and diverse physiological roles of GPR120 has made it a potential target of multiple therapeutic interventions. For effective structure-based drug design better understanding of the role of residues in the orthosteric binding pocket of the receptor and of agonist-induced subtle conformational changes from ligand-protein interactions is needed. In the absence of an X-ray crystal structure of human GPR120, a 3D homology model of human GPR120S [23] was used to perform structure-based pharmacophore screening (SBPS) through ZINC-Pharmer with the ZINC database [24] as summarised in Fig. 1 . All-atomistic molecular dynamics (MD) simulations of 100 ns were performed on the GPR120S receptor models with SBPS selected hits docked at the orthosteric pocket. The MD studies explored the protein-ligand binding modes and their effect on receptor activation to design potential GPR120 antagonists. The Irish Centre for High-End Computing (ICHEC) Kay HPC cluster of 336 nodes (13,440 cores) 2.4 GHz Intel Xeon Gold 6148 (Skylake) processors (www.ichec.ie) was used to perform molecular docking and MD simulations. ICHEC provided 60,000 CPU core units for computational processing and 1000 gigabits (GB) of data storage space. The visual analysis and homology model building were carried out on an in-house 8 node (Intel Ò Core TM i7-4790 CPU @ 3.60 GHz  8) Linux cluster. For visualization and image rendering PyMol Open-source version 2.1.0 [25] and Biovia DS Client visualizer 2019 [26] were used. GPR120S homology model previously generated using template combination of Orexin OX2 receptor (PDB: 4S0V) and Opioid deltalike receptor (4N6H) (see supplementary appendix A) was used for docking the reference ligands. The generated homology model of GPR120S has previously been used to screen chemical libraries of small molecules in combination of with in vitro screening and structure-active relationship profiling [23] . The successful identification of GPR120-targeting small molecules using the generated GPR120S validated the use of GPR120S homology model in the present computer-aided drug discovery. TUG-891 [12] and compound 39 [16] were used as reference molecules for predicting binding poses of ligands in the orthosteric pocket by molecular docking. Discovery Studio Pipeline Pilot version 9.1 [26] was used to generate 3D conformations of TUG-891 and compound 39 [12, 16] with stereoisomers and tautomers, which were then energy minimized in the Avogadro suite with the MMFX96 force field [27] . The energy minimised homology model of GPR120S and ligands were pre-processed in AutoDock4 tools (www.autodock.scripps.edu) [28] to assign torsion angles and add polar hydrogens and charge. AutoDock SMINA [29] was used for rigid-flexible docking. AutoDock4 tools' Grid setting feature was employed to define the binding site grid box, based on the site-specific mutation study of the orthosteric binding pocket residues by Hudson et al 2014 [30] to include residues -R99 (TM2), W104 (ECL1), F115 (TM3), W207, F211 (TM5), W277 (TM6) and F304 (TM7), essential for ligand binding. The grid size dimensions used were 40  60  60 units with spacing of 0.375 Å per unit, and the centre point of the binding pocket set at the xyz coordinates of: x = 61.822, y = 59.75, z = 46.597. The coordinate files of a hydrated, equilibrated 128 molecule POPC lipid bilayer (1-palmitoyl-2oleoyl-sn- glycerophosphocholine) along with lipid parameters for the GRO-MOS 54a7 force field were obtained from the Automated Topology Force Field Builder (ATB) repository [31] . The coordinate file was then resized using Inflategro2 [32] compatible with GROMACS 5.1 [33] to produce a fully hydrated, 512 molecule lipid bilayer. InflateGro2 incorporated LAMBADA [32] to determine the protein's longitudinal configuration using a recursive optimization to test different protein orientations thus aligning the membrane and protein. The membrane protein was then automatically embedded into the lipid bilayer patches and clashing lipids removed. The model was stretched in the plane of the lipid bilayer and energy minimized again on contraction allowing optimization of the lipid/protein interface. The entire lipid bilayer was inflated and then slowly compressed around the protein. Each compression step was followed by a round of steepest descent energy minimization to relax the lipid molecules, keeping the protein restrained. The entire system was then solvated with a singlepoint charge (SPC) water model and neutralized with addition of counter-ions. MD simulations were performed using the GROMACS version 5.1 package with the GROMOS ffG54a7 force field, extended to improve the lipid components of the force field. The topology and other force field parameters for all ligands were obtained from the ATB repository server [31] and were examined manually for any discrepancies with charges. The final lipid bilayer generated consisted of a pre-equilibrated layer of 512 molecules of POPC molecules. The systems were then energy minimized for 10,000 steps using the steepest descent algorithm in the GROMACS package. A 500 ps position-restraining simulation was carried out to restrain the protein by 1,000 kJ mol-1 harmonic restraints to relieve the close contacts with POPC and water under NVT (constant Number of particles, Volume, and Temperature) ensemble conditions, with a Vrescale (modified Berendsen) temperature coupler [33] . This was followed by another 5 ns equilibration run under NPT (constant Number of particles, Pressure, and Temperature) ensemble conditions, before a final production run. The MD systems were run at 310 K, i.e., above the phase transition temperature of pure POPC, to ensure that the lipids maintained their proper density, and 1 bar pressure under isothermal-isobaric ensemble [32, 34] . Nosé-Hoover (which is used widely for membrane NPT simulations [32] ) temperature and Parrinello-Rahman pressure couplers were used to maintain the temperature and pressure values with the protein, ligands, lipids and water (plus ions) molecules coupled separately with a coupling constant of st = 0.1 ps. Semi-isotropic pressure coupling was set with sp = 2 ps, allowing the bilayer to deform in the x-y plane independently of the z-axis. A time-step of 2 fs was used throughout, with periodic boundary conditions. LINCS constraint algorithm was used to maintain the geometry of the molecules [35] . Long-range electrostatic interactions were calculated using the particle-mesh Ewald (PME) method. Van der Waal's interactions and Coulomb interactions were cut off at 12 Å with updates every five steps. Checkpoint files on production runs were saved every 50 ps. The trajectory analysis of protein systems was performed by inbuilt GROMACS tools and visualized in PyMol [25] and XMGRACE (http://plasma-gate.weizmann.ac.il/Grace/). The overall stability of the simulated systems was also checked with respect to temperature, pressure, and potential energy of the systems to check thermodynamic equilibrium during the production simulation runs, confirming the convergence of individual trajectories. The WAD-DAICA webserver was used to compute the protein-ligand binding free energies of protein-ligand snapshots extracted from the MD trajectory every 5 ns using the Binding Affinity by AI module of WADDAICA [36] . The protein-ligand interaction plots from the MD trajectory were generated using Molecular-dynamics-Interac tion-plot (available at https://github.com/tavolivos/Moleculardynamics-Interaction-plot). Multiple linear regression (MLR) was used to correlate the physicochemical descriptors (obtained from SwissADME [37] ) with the binding affinity predictions from WAD-DAICA [36] . The ChemMine tools webserver was used to develop a structural similarity clustering scatterplot [38] . ChemMine clustering tool converts the Tanimoto similarity matrix into distance matrix to cluster the molecules. A six-feature pharmacophore model was generated by manual curation of docked reference structures (TUG-891 and compound 39) using the ZINCPharmer web interface [24, 39] , which was then screened against the ZINC15 database containing $ 230 million purchasable lead-like 3D molecules [24] . The physicochemical and drug-likeness profiling of the screened hits was performed using the online webserver SwissADME [37] . The pharmacophore screening results from ZINCPharmer were downloaded as structure data files (sdf format). The 3D molecules were then prepared and docked against the GPR120S receptor model following the set protocol described in section 2.1. The selected compounds were prepared for 100 ns MD simulations following the protocol described in section 2.2 and 2.3 (see Fig. 1 ) In silico investigation of receptor activation commenced with docking of a known agonist to the receptor model in the inactive state (homology model). Previous molecular docking and sitespecific mutations studies of GPR120S have revealed that a single arginine residue in TM2 (R99) has a critical interaction between the receptor and the -COOH (carboxylate) of its ligands [7, 30] . Six other residues from the orthosteric binding pocket were defined essential for TUG-891 binding and interaction with GPR120S were: W104 (ECL1), F115 (TM3), W207, F211 (TM5), W277 (TM6) and F304 (TM7) following site-specific mutation studies performed by Hudson et al [30] (Fig. 2a, b) . These seven residues were selected as the main criterion for defining the orthosteric binding pocket as well as protein-ligand interactions for binding pose prediction of GPR120S agonists. Over the years, some GPR120 agonists have been developed but to date no orthosteric antagonist of GPR120 is available. For chemical scaffold diversification, besides TUG-891 (EC 50 43 nM) [12, 13, 30] , another known agonist Compound39 -benzofuran propanoic acid analogue (EC 50 97 nM) [16] was used in the molecular docking experiment (Fig. 2c, 2d) . The best docked pose for TUG-891 yielded a binding score À9.87591 kCal/mol (free energy of binding calculated by AutoDock SMINA). The carboxylate of TUG-891 forms a salt bridge with R99 and a strong T-type (perpendicular) pi-stacking interaction between F115 and the cyclic aromatic core structure of TUG-891 stabilize the ligand into the pocket formed between TM3, TM6 and TM7. Other equitable hydrophobic interactions were also observed to stabilize the docked TUG-891 in the orthosteric binding pocket (Table S1 ). A molecular docking study conducted by Hudson [12, 13, 30] used the HM of human GPR120S (based on beta-2 adrenoreceptor template; PDB id: 3P0G) to dock TUG-891. An extensive overlap is observed between the binding pocket resi-dues interacting with the agonist in both studies. Similarly, the best docked pose of Compound39 (binding score À9.82688 kCal/mol) also showed favorable hydrophobic interactions and displayed hydrogen bonds with (R99, W277) two out of the seven specific residues as well as hydrogen bonds (T125, N313) with other residues in the orthosteric binding pocket (Table S2) . The results from literature and the present methodologies provide valuable information concerning the optimal GPR120S model generation. To design new GPR120S selective ligands whether agonists or antagonists, optimal interactions with the non-conserved residues involved in the binding pocket need to be considered. These residues are specific to GPR120S and share proximity to the bound ligands. GPCR activation is initiated by binding of an agonist into the orthosteric binding pocket causing reorganization of interhelical interactions. These structural changes propagate towards the intracellular domain via activation switches as they relay the changes from the inactive to the active state [40] . Here the conformational investigation began with an energy minimized and pressure-volume equilibrated systems of human GPR120S (inactive state) embedded in POPC lipid bilayer and solvated SPC water molecules (Fig. 2a) . The 300 ns MD simulations of three different systems -an apo system and two agonist-bound systems, TUG-891 and Compound39 were performed to gain further insights on the topology of the binding site and agonist orientations over a time period. The protein backbone of starting structures in all the three systems reached stability with no significant changes in RMSD values after 200 ns timestamp (Figure S2 a) . The cysteine bridge between TM2 and ECL2 is highly conserved in most of the Class A GPCRs and is used as an anchor point in modelling the individual backbone of hydrophilic loops as well as forming a ligand trap [34, 40] . The continuous disulfide linkage between the two cysteine residues of GPR120S model (C111 and C194) in TM2 and ECL2 was observed throughout the 300 ns MD production run of all three protein systems (Figure S2 b) supporting the structural stability of the generated GPR120S model. The RMSF (Root Mean Square Fluctuation) analysis of the protein backbone (Fig. 3) illustrates stable helix domains and high rate of fluctuations in the loop regions only, especially the ECL2 (177-204) and ICL3 (236-252) domains of the ligand-bound protein models compared to the apo protein system. The Compound39bound system recorded a marked RMSF difference in the ICL1 (65-71) domain compared to TUG-891-bound system suggesting that both the agonists might be inducing conformational changes by two different mechanisms. The low range of fluctuations in the apo protein system suggest that the generated protein model was in a stable inactive state and remained in the inactive state as it might be in global minima due to absence of ligand induced conformational changes during MD simulation run. The comparative analysis of the putatively active state models (agonist-bound) with the inactive (apo-model) state recorded conformational changes in the residues involved in the ionic lock -R136 of TM3 and D259 of TM6 over the MD production run of 300 ns. The inactive state model retained the closed conformation of ionic lock as both the residues stayed within the range of 2-4 Åan appropriate distance to form and sustain the hydrogen bond or the salt-bridge, making the inactive state of the model stable. As expected, the active state model of TUG-891-bound human GPR120S recorded the two residues drifting apart over the MD trajectory of 300 ns with the ionic-lock broken within the first 10 ns of MD simulation. At the end of the production run the two residues were $ 9-10 Å apart disrupting the salt-bridge. (Fig. 4a) . The Compound39-bound protein model demonstrated unexpected behaviour of staying closer to the inactive conformation. The average distance between the two residues (R136 and D259) remained $ 4 Å throughout the production run (Fig. 4b) . Although the residues were not close enough for salt-bridge formation but with the residues being this close, the intracellular cavity of the receptor might not be able to accommodate the heteromeric G protein. Two possible inferences can be proposed from the Com-pound39-bound protein simulation: 1) The protein model might be in a local minimum of the energy landscape at that specific conformation [41] ; 2) The interaction network pattern of Compound39 carboxylate tail with T125, W277 and N313 (Fig. 4d, 5a) , not observed in TUG-891 docking interactions, led the MD trajectory to different conformational changes at the intracellular domain especially at TM3 -ER 136 M motif (D/ERY) involved in the ''ionic lock" formation. A 300 ns MD re-run from the same starting conformation of Compound39-bound protein system with random initial velocities verified that the variance was reproducible. The comparative analysis of distance between the residues (R136 and D259 -involved in formation of ''ionic lock") from the 300 ns MD rerun reported an increase (average $ 5 Å) with a higher range of fluctuations compared to the first 300 ns MD run (average $ 3 Å) but were significantly less than those of the TUG-891-bound system (average $ 9 Å). The average distance was recorded to be greater than 4 Å limiting the hydrogen bond interactions for the salt bridge / ''ionic lock" formation (Fig. 4c) . A recently published study investigated protein-ligand stability from a 200 ns MD simulation of SR13 (a chromane propionic acid analogue -derived from the same Merck & Co. patented series as Compound39) docked to homology model of GPR120S receptor (Uniprot ID: Q5NUL3-2) [42, 43] . This study highlighted the conformational changes the SR13 ligand adopts to enter the binding pocket [42] . However, to the best of our knowledge, these are the first simulations of GPR120S complexed with ligands (TUG-891 and Compound39) at a 300 ns timescale. As the human GPR120S receptor can bind to the flexible FFAs (PUFAs) as well as a diverse set of rigid compounds such as TUG-891 suggested the existence of different binding conformations in the orthosteric binding pocket of the receptor leading to a cascade effect inducing protein activation. The difference in binding pattern of TUG-891 and Compound39 highlighted important residues (T125, W277 and N313) in the orthosteric binding pocket. Residues T125 and N313 were not analyzed in the previous sitespecific mutagenic study while the W277A mutation resulted in loss of receptor activity [30] . The hydrogen bonding analysis of Compound39-bound protein simulations showed an average of $ 60% H-bond occupancy between Compound39 and the W277 sidechain during the 300 ns production run accompanied with $ 35% and $ 10% occupancy for the N313 and T125 sidechains respectively. The interaction network of W277 and / or N313 with Compound39 might be affecting the conformational changes as observed in TUG-891-bound protein simulations. Based on the inferred significance of the W277 and N313 interaction network with the carboxylate chain of Compound39, a single structure-based pharmacophore model was generated by enumerating the 3D conformation of functional features present in the receptor binding pocket. The docked conformations of TUG-891 and Compound39 were superimposed to generate the pharmacophore model using the ZINCPharmer package [39] which resulted in a six-featured hypothesis consisting of two HBA (Hydrogen bond acceptor), two Ar (Aromatic ring systems) and two Hb (hydrophobic) with preferred chemical features ( Fig. 5 ; Table S3 . The selected hypothesis was generated to focus screening on ligands interacting with the W277 and N313 (Fig. 5b) residues with scaffold features of TUG-891 and Compound39 being added to attain rigidity and better anchorage in the binding pocket by interacting with essential binding residues [13, 30] such as R99, W104, F115, W207, F211 and F304. The pharmacophore-based virtual screening (VS) resulted in 63 unique chemical hits identified from the ZINC15 commercial database. Further, selection of these 63 compounds was performed by analysis of physicochemical descriptors computed using Swis-sADME [37] to predict the druglike and / or lead-like nature of the compounds. In combination with analysis of predicted physicochemical descriptors, an in cerebro assessment was applied to select compounds with diverse and synthesizable scaffolds. Such a well-informed manual selection of ligands after virtual screenings have previously been reported to refine the screening performance [44, 45] . The final screening resulted in 9 best-hits (shown in Table 1 ) based on the SwissADME predictions over compounds' physicochemical properties, lipophilicity, bioavailability, and druglikeness (see supplementary material). These hit compounds were then prepared and docked into the receptor binding pocket following the set protocol (section 2.1) to prepare protein-ligand complexes. The protein-ligand interactions of these 9 compounds were analyzed (Fig. 6a) to confirm W277 / N313 H-bond interaction(s) before preparing the protein-ligand complexes for 100 ns MD production runs following the set protocol (section 2.2 and 2.3). Comparative MD simulations for the 100 ns duration of each of the co-complexes of the 9 hit compounds from ZINC database, TUG-891, Compound39 and Apo-model were carried out to evaluate the significance of the N313 interaction in terms of conformational changes to the protein. A MD time scale of 100 ns was chosen based on observations from other studies where100ns MD simulations were shown to be sufficient to illustrate the dynamics of protein-ligand fingerprinting and induced conformational changes in GPCR molecular studies [42, 47, 48] . The protein-ligand binding affinity of the selected 9 compounds over the Table 1 Assigned nomenclature and docked binding scores (kCal/mol) of the selected pharmacophore-based VS hits based on the SwissADME druglikeness profiling. Log Po/w -oil/water partition co-efficient; Leadlikeness -250 MW 350; XLOGP 3.5; Num. rotatable bonds 7 [37] . The full SwissADME parameters are provided in the supplementary material for Cpds 1- duration of 100 ns MD simulations predicts the strength of the binding interaction between the receptor and the compounds (Fig. 6b) . Cpd 2 and Cpd 8 were predicted to show a continuous decreasing gradient in the binding affinity while Cpd 9 was predicted to show a continuous increase in the binding affinity over time during the MD simulation. Other compounds did not show a consistent pattern of decreasing or increasing binding affinity predictions (Fig. 6b) . The binding affinity predictions from the 100 ns MD studies gave insights to the binding stability of docked hits and non-covalent interactions with residues in the binding pocket. Amongst all the compounds, Cpd 1, 7 and 9 were predicted to conserve W277 and N313 H-bond interactions during the 100 ns MD simulations (Fig. 6c) . As the binding affinities predicted for this study used only the protein-ligand complex (without solvent) snapshots extracted from the MD simulations, it should be cautioned that solvent (water) molecules also play a key role in noncovalent interactions -forming bridge interactions between ligand and protein [49] . Most of the interactions with Cpd9 (electrostatic or hydrophobic) were observed to be continuous during the simulation. Similarly, the protein-ligand binding affinity prediction over the 100 ns MD simulation suggested a stable binding of Cpd9 in the orthosteric binding pocket of GPR120S (Fig. 6b) . In contrast, Cpd 2 showed the lowest binding affinity in the binding pocket compared to the other ligands and Cpd 2 was the only ligand which did not form hydrophobic and /or H-bond interactions with R99 as well as N313 (Fig. 6a and 6b) at the starting conformation. Interestingly, Cpd 7, which had a N313 H-bond interaction (Fig. 6a) at the initial (0 ns) conformation was predicted to have a high binding affinity of À10.8 kCal/mol. However, the binding affinity of Cpd 7 reduced to À10.25 kCal/mol over the MD simulation as the number of N313 H-bond interactions of Cpd 7 reduced (Fig. 6c) . The binding affinity predictions (Fig. 6b) of Cpd 7 at 25 ns and 45-55 ns timeframes reported strong binding affinities (-11.48 kCal/mol and À10.85 to À10.91 kCal/mol respectively) compared to the rest of the MD production run. The protein-ligand interaction fingerprinting of the 100 ns MD production run showed that reformation of H-bonds between Cpd7 and N313 was accompanied by several hydrophobic interactions with other residues in the binding pocket such as F88, F115, W207, I280, I284 and F303. Similar interactions patterns were observed in other ligands (Cpd 1, 4, 5) where an improvement in binding affinities was predicted due to reformation of the N313 H-bond with ligands along with increased hydrophobic interactions. While in-depth analysis of extended MD simulations would further help to understand these changes the present correlative binding affinity predictions and N313 interactions could infer the key role of N313 residue in the ligand binding stability. As GPR120S was modelled from templates in the inactive form, the binding of ligands with antagonistic activity should keep the receptor in inactive state without causing major conformational changes at the intracellular domain of the receptor. The protein's structural stability evaluation by RMSF analysis of 100 ns trajectory showed (Fig. 7a) that Cpd 9 and Cpd 1 stabilized the protein backbone as well as reducing the fluctuations of ECL2 (177-203) and ICL3 (236-252) regions observed in Apo, TUG-891 and Com-pound39 bound protein systems while Cpd 7 bound protein system has recorded the highest range of fluctuations ($12 Å) in the ICL3 region. Binding of other ligand molecules have an overall similar effect on the protein backbone fluctuations, but higher than that of Cpd 9. Focusing on the ligands themselves Cpd 9 and Cpd 7 they were also found to be the most stable in the GPR120S orthosteric binding pocket with RMSD values below $ 1.5 Å and 2 Å respectively throughout the 100 ns MD trajectory (Fig. 7b, S3-S7) . Such low ligand RMSD values suggest that these two ligands are tightly bound to the orthosteric binding pocket without major changes in their initial docked orientations. It is important to mention that both ligands lead to contrasting effects on the protein-backbone RMSD values -Cpd 9 stabilized the protein-backbone in its initial inactive form whereas Cpd 7 deviated the protein-backbone away from its initial conformation. Cpd 9 stabilized the protein within the timespan of the first 20 ns, keeping the average RMSD of protein model below 4 Å. While the binding of Cpd 7 to GPR120S leads to the highest RMSD values ($8 Å). Ligand RMSD analysis of Cpd 1 also presented a range of fluctuations with protein backbone RMSD reaching above $ 5 Å (Fig. 7b) . Compared with other ligands, the Cpd 9 bound protein model predicted the least movement of the distance between R136 and D259 (involved in the ''Ionic lock") at the intracellular domain of GPR120S (Fig. 7c ) that is the site specific for G-protein coupling. The study published by Provasi et al. [50] used inactive and active crystal structures of GPCRs with ligands eliciting different pharma-cological actions performing 20 ns MD simulations. The study reported that depending on their pharmacological activity, the ligand bound to the receptor can shift the conformational equilibrium towards active or inactive state of the receptor. A similar conformational shift (Cpd7 -from inactive to active state) as well as stability (Cpd9 -stabilized inactive state) in receptor state was observed in our study. Here, during the 100 ns MD production run of Cpd9 the ''ionic lock" remained closed (Fig. 7c, S6) , inhibiting the coupling between receptor and G-protein, and hence keeping the receptor in the inactive state. For further analysis, MLR was applied to predict the contribution of physicochemical descriptors (independent variable) such as molecular weight, H-bond donors, H-bond acceptors, logP and topological polar surface area (TPSA) on binding affinity (dependent variable) [51, 52] (Fig. 8) . The average binding affinity (obtained from WADDAICA webserver) of the last 20 ns of MD snapshots (from 80 to 100 ns) was used for the MLR model. The performance of the MLR model is expressed in terms of R 2 , which was found to be 0.799 signifying that $ 80% of the data fit the regression model. Cpd 9 with the highest binding affinity value from the WADDAICA server [36] , was also predicted to have the highest binding affinity by the MLR model. Indeed, the MLR model agreed with the binding affinity predictions of the WADDAICA webserver but it should be noted that many of the variables used in the MLR model were obtained from other prediction algorithms -such as binding affinity, logP and TPSA. The selected compounds Cpd 1-9 were also clustered with the co-crystalized ligands of templates (4N6H-EJ4 and 4S0V-SUV) used for homology model generation and reference ligands used for generation of pharmacophore (TUG-891 and Compound39) to verify the unbiased nature of the screened compounds. The twodimensional scaling cluster ( Fig. 9 ; Table S4 ) confirms that screening of Cpd9 and Cpd7 was not biased towards the templates or reference ligands used in the study. Finally, with the design and discovery of novel scaffolds by in silico methods, the pharmacokinetic profile and synthesizability of these compounds can be a limiting factor. The proposed antagonist Cpd9 has a promising predicted ADME profile indicating leadlikeness but with a moderate predicted solubility (see supplementary material). Furthermore, a feasible retrosynthesis scheme obtained from CAS SciFinder (https://scifinder.cas.org) confirms the ease to synthesis of this compound ( Figure S8 ). Fig. 8 . Scatter plot of a multiple linear regression model of the contribution of physicochemical descriptors -molecular weight, number of H-bond acceptors and H-bond donors, TPSA and logP from SwissADME [37] over the WADDAICA [36] binding affinity for the 9 compounds. The full SwissADME parameters are provided in the supplementary material for Cpds 1-9. Fig. 9 . Clustering of template ligands (PDBs: 4S0V-SUV and 4N6H-EJ4); agonists (TUG-891 and Compound 39) and selected compounds 1-9 by Two Dimensional similarities with a similarity cutoff of 0.4 [38] . V1 and V2 denotes the variable 1 and 2 for the distance matrix obtained from the Tanimoto similarity matrix. In this study, the MD analysis of 300 ns production runs of agonist bound GPR120S models led to the generation of a pharmacophore hypothesis targeting W277 and N313 residues of GPR120S receptor to discover potential candidates as GPR120S antagonists. The structure-based pharmacophore hypothesis was validated by running MD simulations of pharmacophore identified synthesizable hits over a period of 100 ns suggesting that H-bond interactions of Cpd 9 (2-hydroxy-N-{4-[(6-hydroxy-2-methylpyri midin-4-yl)amino]phenyl}benzamide) with W277 and N313 stabilize the occupancy of the ligand in the orthosteric pocket and keeps the protein in the inactive form. While the interactions of a phenylimino-phenol analogue -Cpd 7 phases the protein from inactive to active form due to breakage of the ionic lock which can lead to G-protein coupling at the intracellular domain. Further site-specific mutation studies targeting the N313 residue can confirm the importance of N313 interactions in GPR120S antagonist design. Therefore, the insights from the present study can potentially be employed to enhance the selectivity of GPR120S ligands and target interactions with key residues (W277 and N313) to develop novel agonists and antagonists. Further development of a selective and potent GPR120S antagonist could be useful in management of various cancers where GPR120S has been implicated as a tumour-promoting receptor. The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Trends in GPCR drug discovery: new agents, targets and indications G protein-coupled receptor, an important target for drug design and screening The significance of G proteincoupled receptor crystallography for drug discovery Free fatty acids regulate gut incretin glucagon-like peptide-1 secretion through GPR120 The role of free-fatty acid receptor-4 (FFA4) in human cancers and cancer cell lines Identification of G-proteincoupled receptor 120 as a tumor-promoting receptor that induces angiogenesis and migration in human colorectal carcinoma Concomitant Action of Structural Elements and Receptor Phosphorylation Determines Arrestin-3 Interaction with the Free Fatty Acid Receptor FFA4 Cloning, expression, and pharmacological characterization of the GPR120 free fatty acid receptor from cynomolgus monkey: Comparison with human GPR120 splice variants New Binding Sites, New Opportunities for GPCR drug discovery Differential signaling by splice variants of the human free fatty acid receptor GPR120 International Union of Basic and Clinical Pharmacology. LXXXVIII. G Protein-Coupled Receptor List: Recommendations for New Pairings with Cognate Ligands The Pharmacology of TUG-891, a potent and selective agonist of the free fatty acid receptor 4 (FFA4/GPR120), demonstrates both potential opportunity and possible challenges to therapeutic agonism Discovery of a potent and selective GPR120 agonist Agonism with the omega-3 fatty acids a-linolenic acid and docosahexaenoic acid mediates phosphorylation of both the short and long isoforms of the human GPR120 receptor Free-fatty acid receptor-4 (GPR120): Cellular and molecular function and its role in metabolic disorders Discovery of benzofuran propanoic acid GPR120 agonists: From uHTS hit to mechanism-based pharmacodynamic effects Probe-Dependent negative allosteric modulators of the long-chain free fatty acid receptor FFA4 Development of free fatty acid receptor 4 (FFA4/ GPR120) agonists in health science Fatty acid 16:4(n-3) stimulates a GPR120-induced signaling cascade in splenic macrophages to promote chemotherapy resistance Fatty acid receptor GPR120 promotes breast cancer chemoresistance by upregulating ABC transporters expression and fatty acid synthesis A pathology atlas of the human cancer transcriptome Omega-3 fatty acids and other FFA4 agonists inhibit growth factor signaling in human prostate cancer cells In silico and in vitro screening for potential anticancer candidates targeting GPR120 ZINC 15 -ligand discovery for everyone The PyMol Molecular Graphics System, Schrödinger, LLC. Version 2.1 Discovery Studio Modeling Environment Avogadro: an advanced semantic chemical editor, visualization, and analysis platform AutoDock4 and AutoDockTools4: Automated docking with selective receptor flexibility Lessons learned in empirical scoring with smina from the CSAR 2011 benchmarking exercise The Molecular basis of ligand interaction at free fatty acid receptor 4 (FFA4/GPR120) Testing and validation of the Automated Topology Builder (ATB) version 2.0: prediction of hydration free enthalpies LAMBADA & InflateGRO2: Efficient membrane alignment and insertion of membrane proteins for molecular dynamics simulations Canonical sampling through velocity rescaling Real time monitoring of membrane GPCR reconstitution by plasmon waveguide resonance: on the role of lipids A parallel linear constraint solver for molecular simulation WADDAICA: A webserver for aiding protein drug design by artificial intelligence and classical algorithm SwissADME: a free web tool to evaluate pharmacokinetics, druglikeness and medicinal chemistry friendliness of small molecules ChemMine tools: an online service for analyzing and clustering small molecules ZINCPharmer: pharmacophore search of the ZINC database Diverse activation pathways in class A GPCRs converge near the G-proteincoupling region Minima hopping: An efficient search method for the global minimum of the potential energy surface of complex molecular systems A selectivity study of FFAR4/FFAR1 Agonists by molecular modeling Discovery of chromane propionic acid analogues as selective agonists of GPR120 with in vivo activity in rodents Combining in silico and in cerebro approaches for virtual screening and pose prediction in SAMPL4 Structure based discovery of small molecule suppressors targeting bacterial lysozyme inhibitors The principles of ligand specificity on beta-2-adrenergic receptor Virtual screening, molecular dynamics and structure-activity relationship studies to identify potent approved drugs for Covid-19 treatment Diverse GPCRs exhibit conserved water networks for stabilization and activation Ligandinduced modulation of the free-energy landscape of G protein-coupled receptors explored by adaptive biasing techniques Combined molecular docking and QSAR study of fused heterocyclic herbicide inhibitors of D1 protein in photosystem II of plants Virtual screening and statistical analysis in the design of new caffeine analogues molecules with potential epithelial anticancer activity The authors gratefully acknowledge to Irish Centre for High-End Computing (ICHEC) for the provision of computational facilities and support. The authors wish to acknowledge Technological University Dublin (TU Dublin) for providing financial assistance through Fiosraigh Dean of Graduate Research School Award doctoral fellowship to Ajay Pal (Grant No. PB04185). Supplementary data to this article can be found online at https://doi.org/10.1016/j.csbj.2021.11.005.