Introduction

G protein-coupled receptors (GPCRs) are critical regulators of cellular processes and thus represent prime drug targets (1–3). While traditional GPCR-targeted therapies focus on orthosteric sites, recent advances have revealed allosteric sites offering novel therapeutic avenues (4, 5)(5–7). Although exogenous synthetic allosteric modulators are known, endogenous counterparts remain poorly characterized (8–13). Developing high-affinity endogenous modulators requires integrating structure-based design, artificial intelligence (AI), and assays, yet traditional approaches like SAR analysis are hampered by limited GPCR allosteric modulator data (14, 15). While experimental techniques like FRET and BRET can validate allosteric compounds, their use in high-throughput screening for novel intracellular modulators is challenging (16–18). Identifying endogenous GPCR allosteric modulators is further complicated by factors like incomplete GPCR topology data (19) and vast chemical space (1, 20). This necessitates a hybrid computational approach combining allosteric site prediction, de novo ligand synthesis, and efficient screening, potentially enhanced by AI. Existing de novo drug design tools often lack practical applicability for this purpose due to computational limitations and high technical demands (21–27).

To overcome these challenges, we developed Gcoupler, a software suite (available as a Python package and Docker image) that integrates structural biology, statistical methods, and deep learning to identify GPCR allosteric modulators. We demonstrated the usability and applicability of Gcoupler in identifying novel endogenous modulators of GPCRs by exploiting the α-pheromone (α-factor)-induced mating or programmed cell death (PCD) pathway of S. cerevisiae. Notably, it is well documented that the highly elevated, non-physiological levels of α-factor trigger PCD in less than half of the MATa population (28). Moreover, an equivalent concentration of α-factor triggers distinct PCD kinetics across distinct laboratory strains; for example, BY4741 is more resistant than the W303 strain (29). We, therefore, hypothesized that a subset of pheromone-resistant cells might regulate the Ste2p-mediated PCD signaling via the endogenous intracellular metabolites by operating at the Ste2-Gα binding interface. Using Gcoupler, we identified a subset of intracellular metabolites that could potentially bind to Ste2p (GPCR) at the Gpa1 (Gα) binding interface and obstruct the downstream signaling. Our computational results further suggest that hydrophobic ligands such as sterols strengthen the Ste2p-Gpa1p binding and might trigger a cohesive response that potentially obstructs downstream signaling. Experimental evidence further supported these findings that the elevated intracellular levels of these metabolites rescue the pheromone-induced PCD. To evaluate the evolutionary conservation and possible clinically relevant translation of this mechanism, we tested these metabolites on human and rat isoproterenol-induced, GPCR-mediated cardiac hypertrophy model systems and observed attenuated response in the cardiomyocytes pretreated with GPCR-Gα-protein interface modulating metabolites.

Results

Development and validation of Gcoupler

Designing novel target molecules by integrating the topological, chemical, and physical attributes of protein cavities necessitates advanced neural networks. While existing approaches like Bicyclic Generative Adversarial Networks (BicycleGANs) (30) and Recurrent Neural Network (RNNs) (31) have demonstrated potential, end-to-end standalone tools for GPCR-specific ligand design remain scarce. To address this, we developed Gcoupler, a Python package and Docker image for structure-based, cavity-dependent de novo ligand design, employing advanced Graph Neural Networks (GNNs) and robust statistical methods. Gcoupler consists of four interconnected modules, i.e., Synthesizer, Authenticator, Generator, and BioRanker, all together offering an intuitive workflow for ligand design, screening, and prioritization.

The Synthesizer module identifies putative cavities in protein structures, providing users with flexibility to select cavities based on druggability scores or user-supplied critical residues. Leveraging LigBuilder V3 (32), it generates cavity-specific ligands influenced by topology and pharmacophores, outputting SMILES, cavity coordinates, and other requisite files to downstream modules for further steps (Figure 1a). By default, it synthesizes 500 molecules using a hybrid Growing-Linking mode optimized through genetic algorithms though it can also be user-defined. Notably, the fragment library of LigBuilder, comprising of 177 distinct molecular fragments in Mol2 format, allows the selection of multiple seed structures and extensions that best complement the cavity pharmacophores throughout multiple iterative runs. The Synthesizer module of Gcoupler enhances LigBuilder V3 for GPCR ligand design but lacks user-defined library screening, and proposes synthetically challenging molecules, and often requires post-processing to isolate high-affinity binders from a broad affinity range of synthetically designed compounds.

Figure 1:

To address these challenges, we developed Authenticator which validates Synthesizer outputs by segregating high-and low-affinity binders (HABs and LABs) via virtual screening using the AutoDock Vina (33) and statistically testing the differences in the binding energies distributions using the Kolmogorov–Smirnov, Epps-Singleton, and Anderson-Darling Test. This module refines synthetic ligand pools by analyzing binding energy distributions, facilitating informed segregation while allowing user-defined thresholds and hypothesis testing for flexibility. The Generator, the third module, employs state-of-the-art GNN models such as Graph Convolution Model (GCM), Graph Convolution Network (GCN), Attentive FP (AFP), and Graph Attention Network (GAT) to construct predictive classifiers using Authenticator-informed classes. With automated hyperparameter tuning and k-fold cross-validation, it ensures robust model training and minimizes overfitting. Notably, by default, Gcoupler employs three-fold cross-validation, but users can adjust this parameter.

Finally, BioRanker, the last module prioritizes ligands through statistical and bioactivity-based tools. The first level ranking offered by BioRanker is composed of a statistical tool that encompasses two distinct algorithms, namely G-means and Youden’s J statistics, to assist users in identifying the optimal probability threshold, thereby refining the selection of high-confidence hit compounds (Supplementary Figure 1a). Additionally, bioactivity embeddings computed via Signaturizer (34) enable multi-activity-based ranking using a modified PageRank algorithm. Taken together, Gcoupler is a versatile platform supporting user-defined inputs, third-party tools for cavity selection, and customizable statistical analyses, enhancing its adaptability for diverse ligand design and screening tasks. This integrated framework streamlines cavity-specific ligand design, screening, and ranking, providing a comprehensive solution for GPCR-targeted drug discovery.

To evaluate Gcoupler’s performance, we tested its modules across five GPCRs (AA2AR, ADRB1, ADRB2, CXCR4, and DRD3) using experimentally validated ligands and matched decoys from the DUD-E dataset (35). Gcoupler accurately identified orthosteric ligand-binding sites and additional allosteric cavities across all targets, validating its de novo cavity detection algorithm (Supplementary Figure 1b). For orthosteric sites, the Synthesizer module generated ∼500 compounds per GPCR, and the Authenticator module classified these into high- and low-affinity binders using thresholds optimized for each GPCR (-7 kcal/mol for AA2AR, CXCR4, DRD3; -8 kcal/mol for ADRB1, ADRB2) (Supplementary Figure 1c). Statistical validation confirmed significant separation between these groups (p < 0.0001), enabling the Generator module to construct graph-based classification models with high values of AUC-ROC (>0.95), sensitivity, and specificity (Figure 1b-c, Supplementary Figure d-e). These models reliably distinguished ligands from decoys, demonstrating Gcoupler’s accuracy in identifying high-affinity ligands. We extended this validation to allosteric sites using three GPCR-ligand complexes (β2AR-Cmpd-15PA, CCR2-CCR2-RA-[R], CCR9-Vercirnon) from the PDB (36). Gcoupler successfully identified allosteric binding sites and generated classification models for synthetic compounds with consistently high AUC-ROC values (>0.95) (Figure 1d, Supplementary Figure 2a-c). Projection of experimentally validated ligands onto these models further confirmed their predictive accuracy (Figure 1e), underscoring Gcoupler’s robustness and versatility for orthosteric and allosteric ligand discovery. We next benchmarked Gcoupler’s efficiency against AutoDock using the alpha-1A adrenergic receptor (993 agonists from ChEMBL31). Gcoupler was 13.5 times faster, leveraging its deep learning-based Generator module and AutoDock Vina’s efficiency. Both methods provided comparable predictions for active compounds, demonstrating Gcoupler’s speed and accuracy, making it ideal for large-scale ligand design and drug discovery (Supplementary Figure 3a-h).

Finally we used Gcoupler to evaluate the ligand space conservation (functional conservation) of the GPCR-Gα interface. We analyzed multiple human GPCR-Gα complexes from the PDB (Figure 1f, Supplementary Table 1), identifyed conserved motifs (DRY, CWxL, NPxxY) and binding pockets through sequence and structural analyses (Figure 1g). Structural conservation was quantified by normalized RMSD values (mean: 1.47 Å; median: 0.86 Å), confirmed robust topological similarity at the GPCR-Gα interfaces (Figure 1h-j, Supplementary Figure 3i-k, Supplementary Table 2). Using Gcoupler, ∼50 ligands per GPCR were synthesized, with physicochemical property analysis (via Mordred) and revealed high cosine similarity, which further supports the functional conservation of GPCR-Gα interface (Figure 1h, k, Supplementary Figure 3l, Supplementary Table 3). In summary, we used Gcoupler to systematically evaluate and analyze the ligand profiles of the GPCR-Gα-protein interface and observed a higher degree of sequence, topological, and functional conservation.

Gcoupler reveals endogenous, intracellular Ste2p allosteric modulators

We next utilized Gcoupler to test the hypothesis that the intracellular metabolites could potentially and directly regulate the GPCR signaling by directly interacting with the GPCR-Gα-protein interaction interface (Figure 1l). To test this hypothesis, we utilized a well-characterized yeast mating pathway mediated via the Ste2p-Gpa1p interface (Supplementary Figure 4a). Gcoupler was employed to screen metabolites from the Yeast Metabolome Database (YMDB) against the yeast phospholipid bilayer-embedded molecular dynamics simulated cryo-EM structure of Ste2p (Supplementary Figure 4b-i). This led to the identification of 17 potential surface cavities on Ste2p, with two intracellular regions, IC4 and IC5, accounting for over 95% of the Ste2p-Gpa1p interface (Figure 2a-b, Supplementary Figure 4j-n). Conservation analysis of these cavities across 14 yeast species further confirmed their structural significance (Supplementary Figure 4o). We independently ran Gcoupler for IC4 and IC5, compared their chemical diversity to ligands for the extracellular cavity (EC1), revealing distinct pharmacophore properties, and confirmed cavity-specific chemical spaces through principal component analysis (PCA) (Supplementary Figure 5a-e). We further assessed the reproducibility of the Gcoupler by synthesizing 100 in silico compounds per run across five runs for the Ste2p IC4 cavity. Analysis of chemical diversity via 2D/3D PCA and pairwise Tanimoto Similarity (atom-pair fingerprints) revealed heterogeneous yet overlapping chemical compositions among ligands across runs (Supplementary Figure 5f-h).

Figure 2:

Post these performance/reproducibility checks of Gcoupler, we finally generated high-affinity (HABs) and low-affinity (LABs) binders, screened yeast metabolites using a hyperparameter-tuned Graph Neural Network (Attentive FP), and predicted potential binders to IC4 and IC5 with a >0.5 binding probability. AutoDock-based molecular docking of YMDB metabolites against the cavities further refined our list, revealing metabolites with consensus binding energies ≤ -7 kcal/mol and binding probabilities >0.5 (Figure 2c-f, Supplementary Figure 6a-f, Supplementary Table 4). Furthermore, HAMs (High-Affinity Metabolites) and LAMs (Low-Affinity Metabolites) displayed distinct atomic fingerprints, with enriched functional groups, including R2NH, R3N, ROPO3, ROH, and ROR, observed in HAMs (Supplementary Figure 6g-h). Pathway enrichment analysis identified metabolites involved in steroid, sesquiterpenoid, and triterpenoid biosynthesis, along with folate pathways (Figure 2g). BioRanker module further pinpointed sterols, including zymosterol (ZST), ubiquinone 6 (CoQ6), and lanosterol (LST), as top candidates, exhibiting high prediction probabilities (>0.99) and structural similarity to HABs (Figure 2h, Supplementary Figure 6i-j). To validate Gcoupler-identified allosteric modulators, we performed control analyses of the Authenticator and Generator modules and blind docking of YMDB metabolites with Ste2p. Removing HAB and LAB labels from IC4 cavity compounds and applying random data splitting yielded poor model performance compared to Authenticator-guided splitting, confirming its robustness (Supplementary Figure 6k). Increasing training data size (25%, 50%, 75%, 100%) improved model accuracy for YMDB metabolite projections (Supplementary Figure 6l-m). Blind docking using AutoDock failed to distinguish HAM and LAM, unlike cavity-specific analysis via AutoDock and Gcoupler, which showed clear segregation (binding energy cutoff: -7 kcal/mol, Gcoupler probability: 0.5; Supplementary Figure 6n). All these rigorous control analyses and blind docking validated the reliability of Gcoupler’s predictions, confirming its robustness in identifying cavity-specific modulators.

To experimentally validate the role of Gcoupler-predicted metabolites in Ste2p signaling, we performed a genetic screen of metabolic mutants. Mapping the predicted metabolites to KEGG and MetaCyc pathways revealed key enzymes involved in their metabolism. Screening 53 single metabolic deletion mutants showed that most deletion mutants resisted α-factor-induced PCD, with some displaying accelerated growth, suggesting a metabolic modulation of Ste2 signaling (Figure 2i-k, Supplementary Figure 7a-c, Supplementary Table 5). Notably, the STE2 knockout mutant (ste2Δ) displayed no significant PCD induction, validating the importance of Ste2p in this process. To further investigate the metabolic pathways associated with PCD resistance, we conducted high-resolution metabolomics on yeast cells exposed to varying α-factor concentrations. Survivors of PCD induction exhibited differential metabolite profiles, with CoQ6 and ZST showing consistent enrichment across conditions (Figure 2l-n, Supplementary Figure 7d-j, Supplementary Table 6). Pathway analysis of these metabolites indicated involvement in glyoxylate, dicarboxylate, purine, and vitamin B6 metabolism (Supplementary Figure 7k). These findings, in conjunction with the genetic screen, highlight the interplay between metabolism and Ste2 signaling, with computationally predicted metabolites like ZST and CoQ6 potentially conferring resistance to α-factor-induced PCD.

Elevated Endogenous Metabolic levels selectively inhibit GPCR signaling

To evaluate the stability of the interactions of zymosterol (ZST), lanosterol (LST), and ubiquinone 6 (CoQ6) at the Ste2p-Gpa1p interface, we performed molecular dynamics (MD) simulations (100 ns) for Ste2p-metabolite complexes. Results revealed thermodynamically stable interactions across cavities IC4 and IC5, except for CoQ6, which exhibited fluctuating RMSD in IC5 (Figure 3a-c, Supplementary Figure 8a-b). Extended 550 ns simulations confirmed stability for ZST and LST in IC4, while CoQ6 showed RMSD fluctuations post 100 ns (Supplementary Figure 8c). To gain further insight into the contributing residues from the MD simulations, we performed a residue-wise decomposition analysis that provides information about the energy contributions from the different residues to total binding free energies (Supplementary Figure 9a). These results suggest that IC4 and IC5 specific residues predominantly contribute to the total binding free energies. Notably, the binding free energies are obtained as an average of over 500 configurations corresponding to the last 50 ns of the MD simulations.

Figure 3:

To gain deeper insight into the mode of action of these metabolites in inhibiting Ste2p signaling, we first analyzed their impact at the orthosteric site, i.e., α-factor binding. We performed protein-peptide docking of Ste2p and α-factor and observed that metabolite-binding at the Ste2p-Gpa1p interface favors the α-factor interaction, as inferred from the binding free energies ΔG (kcal/mol) (Supplementary Figure 9b-c). We also analyzed protein-protein interaction between the Ste2p (GPCR) and miniGpa1p (55) by selecting GPCR configurations with and without metabolite-induced altered cavity topologies (IC4 and IC5), respectively (frames output from the aforementioned Ste2p-metabolite complex simulations) and computed the dissociation constant (Kd), binding affinity (ΔG) and the structural changes in the overall Ste2p topology (Figure 3d-e, Supplementary Figure 9b-d). These computational analyses revealed that in contrast to metabolite-free Ste2p-Gpa1p interaction, referred to as wild-type (WT) condition, the Kd value is multi-fold lower in the presence of metabolite, indicating a cohesive response induced by these metabolites. A multi-fold lower Kd value further indicates and potentially explains that the metabolites binding favors Ste2p (GPCR) and miniGpa1-protein interaction and enables the establishment of a stable complex that might influence the shielding of the effector-regulating domains of the Gpa1p or influence its binding with the Ste4p (Gβ)-Ste18p (Gγ) complex.

We next asked whether the observed resistance toward α-factor-induced PCD in single metabolic mutants might be the direct consequence of the identified metabolites or a pleiotropic response due to an altered genome. To test this, we next exogenously elevated ZST, CoQ6, and LST levels (prediction probability > 0.99 and post-filtration using BioRanker) in wild-type yeast and observed selective modulation of Ste2p signaling. Pre-loaded cells exhibited rescue from α-factor-induced PCD across growth kinetics, PI-based viability, and FUN^TM1 assays while maintaining α-factor specificity over acetic acid-induced PCD (Figure 3f-g, Supplementary Figure 10a-c). Mating assays further supported reduced Ste2p signaling in metabolite-preloaded cells, with CoQ6 showing no significant effects, consistent with its instability observed in molecular dynamics simulations (Figure 3h, Supplementary Figure 10d). Further, we evaluated the deactivation of the pathway at the MAPK signaling level, where the α-factor-induced p-Fus3 levels were suppressed by ZST and LST (Figure 3i) but not CoQ6 (data not shown). Moreover, by using the P_FUS1-eGFP-CYC1 strain-based reporter assays, metabolite pre-treatment confirmed the inhibition of α-factor signaling (Figure 3j). These findings suggest endogenous metabolites as modulators of Ste2p signaling by stabilizing Ste2p-Gpa1p interactions.

Site-Directed Ste2 Mutants Abrogate Metabolite-Mediated Rescue Phenotype

To gain a deeper understanding, we investigated the role of metabolite binding in modulating Ste2p-Gpa1p interaction dynamics and employed both computational and experimental approaches. In silico site-directed mutagenesis of Ste2p revealed key metabolite-binding residues: S75 and R233 for CoQ6, L289 for ZST, and T155, V152, and I153 for LST. Prioritizing strong hydrogen bonds within an ideal 2.7–3.3 Å range, we selected mutants S75A, L289K, and T155D. These mutations significantly increased the dissociation constant (Kd) of the Ste2p-Gpa1p complex, indicating weakened interactions compared to wild-type Ste2p (Figure 4a-b, Supplementary Figure 11a). This computational evidence supports the hypothesis that metabolite binding stabilizes the Ste2p-Gpa1p complex, facilitating the rescue response to α-factor-induced programmed cell death (PCD).

Figure 4:

Next, we performed experimental validation by generating site-directed mutants and confirmed these computational predictions. Mutants (S75A, L289K, T155D) expressed in a ste2Δ background (Supplementary Figure 11b) exhibited no rescue from α-factor-induced PCD upon pre-treatment with their respective metabolites, as shown by fluorometry-based cell death assays, FUN^TM1 staining, and p-Fus3 signaling analysis by mapping MAPK pathway activation (Figure 4c-e). The lack of rescue highlights the direct role of metabolite binding at the Ste2p-Gpa1p interface in regulating downstream signaling. Interestingly, S75A mutants showed no α-factor-induced effects, likely due to significant structural disruptions. Shmoo formation assays further corroborated these findings, with no rescue effects observed in the mutants despite metabolite pre-loading (Supplementary Figure 11c-e).

To gain mechanistic insights, we performed RNA sequencing on LST-preloaded and untreated (control) cells with and without α-factor exposure. Transcriptomic analysis revealed differential expression of genes implicated in PCD and mating responses, including GSY1, whose downregulation was linked to α-factor resistance (Figure 4f-g, Supplementary Figure 12a-c, Supplementary Tables 7-10). Loss-of-function studies using knockout strains validated the involvement of key transcripts in metabolite-mediated rescue, particularly YCR095W-A, which exhibited a pronounced phenotype loss (Supplementary Figure 12d-f). These results collectively suggest that metabolite binding at the Ste2p-Gpa1p interface directly drives rescue responses, with secondary contributions from differentially expressed genes in attenuating α-factor-induced cell death.

Finally, we explored whether intracellular allosteric modulators such as ubiquinone 6 (CoQ6), zymosterol (ZST), lanosterol (LST), fucosterol (FST), and ubiquinone 10 (CoQ10) modulate GPCR signaling in higher vertebrates such as human and rat beta 1/2-adrenergic receptors signaling. Notably, computational docking at the GPCR-Gα interface revealed strong binding affinities for these metabolites, similar to Ste2p-metabolite interactions (Supplementary Figure 13a-b). Sequence conservation analysis of the GPCR-Gα interface across yeast Ste2p and adrenergic receptors in humans and rats further confirmed a high degree of evolutionary conservation at the metabolite-binding residues (Supplementary Figure 13c). Finally, to test this functional relevance, we evaluated the effect of these metabolites on isoproterenol-induced adrenergic receptor-mediated cardiac hypertrophy in human AC16 cardiomyocytes and neonatal rat cardiomyocytes. Preloading cells with these metabolites significantly attenuated hypertrophic responses, as evidenced by reduced single-cell surface area in quantitative assessments (Figure 4h-j). Notably, to further evaluate the evolutionary conservation of this phenomenon, we also analyzed 75 unique GPCR-Gα complex structures from six species, selected from the PDB database. Dynamic docking was performed using five metabolites (CoQ6, ZST, LST, FST, and CoQ10) identified by Gcoupler as potential allosteric modulators and five negative controls predicted as poor binders. Results revealed significantly lower docking scores (<-7 kcal/mol) for Gcoupler-recommended metabolites compared to negative controls, irrespective of GPCR type or species (Supplementary Figure 13d-h). These findings demonstrate that intracellular metabolite modulation of GPCR activity is a conserved mechanism extending beyond yeast to higher vertebrates.

Discussion

Over the last few decades, extensive research has focused on identifying allosteric modulators of GPCRs due to their relevance in drug discovery (7, 20). Most known modulators are exogenous and target extracellular sites, while intracellular allosteric sites, identified recently through structural biology, offer novel avenues for regulation (3, 4, 6–8, 20). These sites, overlapping with G-protein and β-arrestin coupling regions, highlight the potential for intracellular allosteric modulation (1, 5, 8). Intracellular modulators, including chemically diverse agents like auto-antibodies and sodium ions, remain poorly understood, emphasizing the need for systematic exploration of these sites (4–6, 8, 9, 11–13, 37). However, the lack of data on intracellular modulators limits the feasibility of conventional computational approaches (14, 19, 38)(39–42).

To address this gap, we developed Gcoupler, a computational framework integrating de novo cavity identification, ligand synthesis, statistical analysis, graph neural networks, and Bioactivity-based ligand prioritization. Unlike existing tools, Gcoupler does not require cavity-specific experimentally validated compounds for model training. Gcoupler’s precision in cavity mapping, flexibility for user-defined queries, and ability to screen large chemical libraries make it a versatile and efficient tool. Additionally, Gcoupler’s generic design allows application beyond GPCRs, contrasting with existing platforms that often have limitations in modularity, precision, or open-source availability. Noteworthy, in contrast to other known allosteric sites identification tools for GPCRs, such as Allosite (43), AllositePro (44), AlloReverse (45), that largely leverage the ML-based models or require orthosteric ligand-bound structure as input, the cavity detection feature of Gcoupler (LigBuilderV3), a critical step for the entire workflow, is not limited to only the allosteric sites; instead, it identifies all possible cavity-like regions on the protein surface, which then gets classified into druggable, undruggable, or amphibious based on their individual scoring and ligandability, thus making it unbiased and more specific towards query protein. Notably, the rationale for opting for LigBuilder V3 for cavity identification over similar tools such as Fpocket (46) is that the former uses a hydrogen atom probe, moving along the protein surface grid of 0.5 Å for cavity detection, being much more precise in detecting cavity boundaries, in both breadth and depth mapping; in contrast, the latter considers clusters of alpha spheres (Supplementary Table 11).

To date, only a few methods leverage generative AI models for cavity/pocket-based drug design. Gcoupler is an open-source, end-to-end platform integrating Ligand-Based Drug Design (LBDD) and Structure-Based Drug Design (SBDD) for drug design and large-scale screening. Unlike Pocket Crafter (47), which requires proprietary tools (e.g., MOE QuickPrep) and lacks predictive model-building modules, Gcoupler offers comprehensive functionality. Similarly, DeepLigBuilder (23) and Schrodinger’s AutoDesigner are either closed source or limited in features compared to Gcoupler. Comparative analysis highlights Gcoupler’s unique advantages in precision, flexibility, and functionality (Supplementary Table 11).

Using Gcoupler, we investigated the molecular basis of innate resistance to α-factor-induced programmed cell death (PCD) in yeast. Unlike humans, yeast possess only two GPCR systems, making their pheromone-sensing pathway ideal for focused study. Previous research predominantly identified downstream regulatory mechanisms (29, 48, 49); however, our findings suggest an upstream, receptor-level regulation via endogenous intracellular metabolites. Computational and experimental evidence pinpointed specific metabolites binding to the Ste2p-miniGpa1 interface, modulating signaling. Site-directed mutagenesis confirmed the functional relevance of these metabolite-interacting residues. Of note, previous mutagenesis experiments also revealed multiple critical amino acid residues that overlap with IC4 and IC5, suggesting their functional relevance in Ste2p downstream signaling (Supplementary Table 12). Despite these advances, certain limitations remain, including the potential pleiotropic effects in metabolic gene knockouts and challenges in replicating natural metabolite concentrations. Mechanistic insights into sterol biosynthesis mutants (ergΔ) revealed impaired mating responses due to heterogeneous defects, such as reduced sterol accumulation and shmoo formation, impaired membrane fusion, and decreased FUS1 expression (50–54). This highlights a novel role for sterols in GPCR regulation and their broader implications for yeast microbial factories and stress tolerance (55, 56). Sterols and other endogenous metabolites were shown to modulate GPCR activity by targeting conserved Gα-binding sites, reinforcing the evolutionary conservation of this mechanism across species, as demonstrated in human and rat hypertrophy models in this study.

In summary, our work uncovers a novel regulatory mechanism for GPCRs mediated by intracellular metabolites and presents a computational framework, Gcoupler, to explore unexplored allosteric sites. The proposed model suggests that selective metabolites binding to GPCR-Gα interfaces induce local conformational changes, stabilizing GPCR-G-protein complexes and potentially obstructing downstream signaling. Alternative mechanisms, such as orthosteric site modulation, kinase/arrestin interaction interference, or alterations in membrane dynamics, remain to be explored, warranting further investigation. A critical limitation of our study is the absence of direct binding assays to validate the interaction between the metabolites and Ste2p. While our results from genetic interventions, molecular dynamics simulations, and docking studies strongly suggest that the metabolites interact with the Ste2p-Gpa1 interface, these findings remain indirect. Direct binding confirmation through techniques such as surface plasmon resonance, isothermal titration calorimetry, or co-crystallization would provide definitive evidence of this interaction. Another critical limitation of our findings is its reliance on tools like AutoDock and PRODIGY for preliminary binding affinity estimates, which lack the thermodynamic precision of advanced methods. To address this, we employed MD simulations with MM/GBSA, incorporating factors like protein flexibility and solvation effects for more accurate ΔG calculations. While computationally intensive approaches were beyond this study’s scope, we ensured reported ΔG values reflected system conformational flexibility by basing them on pre-simulated docked structures from molecular dynamics simulations. Further, our results suggest that the metabolite binds to the Ste2p-Gpa1 interface and modulates receptor activity upon pheromone stimulation, as supported by various assays. However, the precise sequence of interactions between Ste2p, the metabolite, and Gpa1 remains unexplored, as it requires sequential experiments beyond this study’s scope. Taken together, addressing these limitations in future work would significantly strengthen our conclusions and provide deeper insights into the precise molecular mechanisms underlying the observed phenotypic effects.

Material and methods

Backend code for the Gcoupler

The back-end code for Gcoupler is implemented entirely in Python (3.8) and comprises four modules: Synthesizer, Authenticator, Generator, and BioRanker. Synthesizer employs LigBuilder V3.0 (32) for de novo in silico ligand synthesis, identifying protein cavities likely to be active or allosteric sites using a hybrid GROW-LINK approach with a Genetic Algorithm. The module autonomously selects one cavity for ligand synthesis based on user-defined residue positions. Using the CAVITY function of LigBuilder (27), it classifies 3D grid points around the protein into occupied, vacant, and surface points and integrates geometric and physicochemical properties to identify binding sites (57). Synthesizer outputs ligand structures in SMILES and PDBQT formats, alongside cavity grid coordinates for downstream modules. Authenticator validates the synthesized ligands using AutoDock Vina (1.2.3) for virtual screening (33). Binding energy calculations classify ligands into High-Affinity Binders (HABs) and Low-Affinity Binders (LABs), preserving balance for subsequent deep learning analysis. The default energy threshold is set to -7 kcal/mol, but users can explore alternative cutoffs, visualize distributions, and perform statistical comparisons within the workflow. These steps ensure precise identification and prioritization of ligand candidates for further analysis. The Authenticator uses the Kolmogorov-Smirnov test (58), Anderson-Darling test (59), and Epps-Singleton test (60) for hypothesis testing for the comparison of the distributions. The Authenticator module visualizes ligand distributions using overlapping density plots and Empirical Cumulative Distribution Function (ECDF) curves. If the default threshold fails to produce statistically meaningful separation, users can supply alternative negative datasets, such as decoys generated via Gcoupler’s inbuilt RDKit Chem module or custom datasets (61). The Generator module builds deep learning-based classification models using the DeepChem (2.6.1) library (62). It accepts High-Affinity Binders (HABs) and Low-Affinity Binders (LABs)/decoys from the Authenticator module to train four graph-based models: Graph Convolution Model (GCM), Graph Convolution Network (GCN) (63), Attentive FP (AFP) (64), and Graph Attention Network (GAT) (65). Class imbalance is addressed through upsampling techniques. Generator tests all models using default hyperparameters, returning performance metrics for user selection. Hyperparameters can be tuned via manual settings or using default values followed by k-fold cross-validation. The final optimized model, trained on the complete synthetic dataset (HAB + LAB/decoys), enables large-scale screening of user-supplied compounds based on their SMILES representations. The BioRanker module performs post-prediction analysis for functional activity-based compound screening. Positively predicted compounds are selected using a stringent probability threshold or adaptive methods such as G-means and Youden’s J statistic, which optimize sensitivity and specificity. The selected compounds are projected into biological activity spaces (Chemistry, Targets, Networks, Cells, Clinics) by comparing their biological activity descriptor vectors with those of HABs using cosine similarity (34). A modified PageRank algorithm ranks compounds based on activity-specific scores, with support for multi-activity ranking to refine results based on user-defined biological properties, ensuring precise and context-relevant compound prioritization.

Additional information about the backend code for Gcoupler along with methodolody for Runtime analysis, Sequence-structural-functional level analysis, Molecular Dynamics Simulation, Molecular docking (AutoDock), Functional enrichment analysis, and Protein-protein docking can be accessed in the Supplementary Information.

Gcoupler Benchmarking

To assess batch effects across Gcoupler runs for a specific cavity, we utilized the standard Gcoupler Docker image. Intracellular cavity 4 (IC4) of the Ste2 protein of yeast was used for benchmarking. A total of 100 molecules were in silico synthesized by the Synthesizer module of Gcoupler iteratively. Post-generation, atom pair fingerprints (ChemmineR; R package) were calculated for the synthesized molecules from each run, and the data was visualized using Principal Component Analysis and pairwise comparison using Tanimoto Similarity (ChemmineR, R package).

For model benchmarking, Gcoupler was validated on GPCRs from the DUD-E dataset, alongside the information about the active ligands and their randomly selected number-matched decoys (35). Additionally, Gcoupler’s performance in identifying experimentally elucidated allosteric sites and modulators was tested using PDB complexes obtained from the RCSB PDB database.

Metabolomics

Wildtype (BY4741) and ste2Δ yeast strains were grown in YPD medium at 30°C, 150 rpm, through primary and secondary cultures (16 hours each). Equal cell numbers (1.5 mL) were aliquoted into a 96-well deep well plate. α-factor (Sigma-Aldrich) was added at final concentrations of 10, 20, 30, 40, and 50 μM (eight replicates each). DMSO served as the solvent control, while untreated WT and ste2Δ conditions received no treatment. Plates were incubated (30°C, 150 rpm, 4 hours) under a breathable membrane. A 50 μL aliquot was taken for Propidium Iodide (PI; 11195, SRL) assay as described in Supplementary Information. Following the PI assay, four pooled replicates were pelleted (6000 rpm, 5 min, RT), treated with zymolyase (40 U/mL, 1X PBS, 30°C, 1 hour), washed with PBS, and metabolomics analysis was performed. Data analysis included peak normalization, omission of metabolites with constant or >50% missing values, and kNN-based imputation (MetaboAnalyst). Data were IQR-filtered, and differentially enriched metabolites (DEMs) were identified by calculating log2 fold change (|log2FC| ≥ 1, p < 0.05 via Student’s t-test). Pathway Over-Representation Analysis (ORA) was performed using MetaboAnalyst with hypergeometric or Fisher’s exact tests to assess pathway enrichment against background metabolite distributions. Further details about methodology is available in Supplementary Information.

Genetic Screening

Fifty-three knockout strains from the Yeast Deletion Collection, along with WT and ste2Δ, were treated with α-factor (30 μM), while DMSO served as the solvent control. Plates were incubated for 4 hours. A 50 μL aliquot was used to measure Propidium iodide-based cell viability assessment assay as described in Supplementary Information. Fluorescence data were normalized to blank-adjusted OD600, followed by two additional rounds of normalization with unstained and HK controls. Percentage fold change for the treated group was calculated relative to the untreated group, and statistical significance was determined using a one-sample Student’s t-test. Further details about the methodology used is available in Supplementary Information.

Pre-loading of Yeast cells with metabolite

Yeast cells were cultured in YPD medium at 30°C, 200 rpm for 16 hours in primary and secondary cultures. Equal cell densities (5 μL) from secondary cultures were inoculated into 96-well plates containing 145 μL YPD with metabolites Coenzyme Q6 (CoQ6, 900150O, Avanti^® Polar Lipids), Zymosterol (ZST, 700068P, Avanti^® Polar Lipids), and Lanosterol (LST, L5768, Sigma-Aldrich) at 0.1 μM, 1 μM, and 10 μM concentrations. Plates were incubated for 24 hours at 30°C, 200 rpm, with multiple biological replicates. Ethanol-treated wells served as solvent controls. For site-directed STE2 mutants, the mutants were grown in YPD for primary and secondary cultures, but the metabolite pre-loading was performed in YPGR instead of YPD to induce Ste2 expression. After pre-loading, the following assays were performed: growth kinetics, propidium iodide-based assay, FUN^TM1 staining, Mating assay, Phospho-MAPK activity based-western blot and Transgenic reporter assay. The detailed protocol for each of these assays is available in Supplementary Information.

RNA-Sequencing

Yeast cells (BY4741) were cultured in YPD medium at 30 °C, 200 rpm, with or without lanosterol (LST, 1 μM) in biological duplicates. Cells were subsequently treated or untreated with α-factor (30 μM) for 2 hours. RNA was isolated following Mittal et al. (66). Sequencing quality was assessed using MultiFastQ, and paired-end reads were trimmed and aligned to the S. cerevisiae reference genome (ENSEMBL R64-1-1; GCA_000146045.2) using the Rsubread package (v2.6.4). Gene expression counts were generated via featureCounts and normalized using TMM (Supplementary Table 7). Differential expression analysis was performed with NOISeq (v2.38.0) (Supplementary Tables 8–10). Functional enrichment was assessed using the Gene Ontology Term Finder (v0.86) from the Saccharomyces Genome Database. Raw FASTQ files and normalized expression data are available on Zenodo. Additional details are available in the Supplementary Information.

Site-Directed Mutagenesis

The gene encoding wild-type STE2 was PCR amplified from the genome of Saccharomyces cerevisiae and further cloned into plasmid pRS304 under galactose-inducible GAL1 promoter to generate plasmid pRS304-P_GAL1-STE2-CYC1 using Gibson assembly (67, 68). The mutants were generated by PCR amplifying the gene with primers consisting of respective mutations and cloned into plasmid pRS304 under GAL1 promoter to generate pRS304-P_GAL1-STE2 S75A-CYC1, pRS304-P_GAL1-STE2 T155D-CYC1, and pRS304-P_GAL1-STE2 L289K-CYC1. Wild-type (rtWT) and mutant STE2 were integrated by digesting the pRS304 vector with restriction enzyme BstXI to generate linearized plasmid and transformed into Saccharomyces cerevisiae BY4741 ste2Δ strain. Additional details about the methodology is available in Supplementary Information

Cardiomyocytes Hypertrophy Models

Human AC16 cardiomyocytes were cultured in DMEM-F12 (Thermo Scientific) with 12.5% fetal bovine serum (FBS) at 37°C and 5% CO2. Cells were seeded in a 24-well plate for size measurements, treated after 24 hours with metabolites (CoQ6, ZST, LST, FST, CoQ10) at 2.5 μM, and incubated overnight with 1% FBS. The medium was refreshed with fresh metabolites and isoproterenol (25 μM) for 48 hours. Cells were washed with PBS, fixed with 4% paraformaldehyde, and stained with wheat germ agglutinin (Thermo Scientific) and DAPI. Images were captured using a Leica DMI 6000 B microscope at 20X magnification, and cell area was measured using ImageJ. Neonatal rat cardiomyocytes were isolated from 1-3-day-old SD rat pups using Collagenase Type II. After heart explantation and digestion, the cells were centrifuged and pre-plated for 90 minutes to remove fibroblasts. The cardiomyocytes were seeded in a gelatin-coated 24-well plate, incubated overnight with 2.5 μM metabolites and 1% FBS, and then treated with metabolites (2.5 μM) and isoproterenol (10 μM) for 72 hours. Cells were fixed, stained with alpha-sarcomeric actinin and DAPI, and images were captured using a Leica DMI 6000 B at 20X magnification. Cell area was quantified using ImageJ. Additional details about methodology is available in Supplementary Information.

Statistical Analysis

Statistical analyses were performed using Past4 software or R-Programming. The Mann-Whitney U test was applied to compare medians between two distributions (non-parametric), while Student’s t-test was used for pairwise comparisons of means. P-value correction was performed using the Bonferroni method when necessary. A significance threshold of 0.05 was set, with *, **, ***, and **** indicating p-values <0.05, <0.01, <0.001, and <0.0001, respectively.

Data Availability

The processed untargeted metabolomics data is provided as Supplementary Information. The raw RNA sequencing files are available at ArrayExpress under accession E-MTAB-12992.

Code Availability

A Python package for Gcoupler is provided via pip https://test.pypi.org/project/Gcoupler/. A docker container pre-compiled with Gcoupler and all of its dependencies can be found at https://hub.docker.com/r/sanjayk741/gcoupler. The source code of Gcoupler is available on the project GitHub page: https://github.com/the-ahuja-lab/Gcoupler and also at Zenodo with DOI: 10.5281/zenodo.7835335, whereas the raw sequencing files can be accessed using DOI: 10.5281/zenodo.7834294.

Ethical statement

The local IAEC (Institutional Animal Ethics Committee) committee at CSIR-Central Drug Research Institute approved all the animal experiments (IAEC/2020/38) following the guidelines of the Committee for the Purpose of Control and Supervision of Experiments on Animals (CPCSEA), New Delhi, Government of India.

Acknowledgements

The authors thank the IT-HelpDesk team of IIIT-Delhi for assisting with the computational resources. We thank all the members of the Ahuja lab for their intellectual contributions at various stages of this project. We thank Prof. G.P.S Raghava for providing critical comments on our manuscript. We thank Dr. Arjun Ray for providing intellectual support. We also thank Dr. Martin Graef and Dr. Kaushik Chakraborty for sharing yeast strains and the NIPER Guwahati central facility for helping us with high-resolution metabolomics. The Ahuja lab is supported by the Ramalingaswami Re-entry Fellowship (BT/HRD/35/02/2006), a re-entry scheme of the Department of Biotechnology, Ministry of Science & Technology, Government of India, Start-Up Research Grant (SRG/2020/000232) from the Science and Engineering Research Board, and a research grant from IHUB Anubhuti (Project Grant/23) and an intramural Start-up grant from Indraprastha Institute of Information Technology-Delhi. The INSPIRE faculty grant from the Department of Science & Technology, India, funds the Sengupta lab. Gupta Lab is funded by the Ramalingaswami Re-entry Fellowship (BT/RLF/re-entry/14/2019) from the Department of Biotechnology, Government of India.

Additional information

Author Contributions

The study was conceived by G.A. Computational analysis workflows were designed by G.A. Statistical guidance was provided by D.S. Docking, and Molecular Simulations were supervised by N.A.M and conducted by N.A.M, and S.K.M. Data analyses was performed by S.K.M., A.M., S.S., S.K., S.D., A.K.S., S.A., V.G., T.S.G., N.K.D and K.S.. S.K.M built the Gcoupler Python package. G.A. and A.M. designed yeast experimental workflows and performed by A.M. Site-directed yeast mutants were created by N. and D.S. Human cell culture-based experiments were performed and analyzed by A.G., and S.K.G. Illustrations were drafted by A.M., S.K.M., and G.A. G.A. wrote the paper. All authors have read and approved the manuscript.

Additional files

Supplementary Information. The supplementary Information provides details of the nucleotide and amino acid Multiple Sequence Alignments (MSAs) of the wild-type STE2 and the site-directed missense mutants of STE2. Each row starts with the sequence identifier followed by the aligned sequence (in chunks) with the ending position of the aligned sequence provided at the end, separated by tabs. The mutation sites are highlighted in red (S75A), yellow (T155D) and blue (L289K). Supplementary Information also contains detailed materials and methods.