Exploring TAS2R46 biomechanics through molecular dynamics and network analysis

Understanding the intricate interplay between structural features and signal-processing events is crucial for unravelling the mechanisms of biomolecular systems. G protein-coupled receptors (GPCRs), a pervasive protein family in humans, serve a wide spectrum of vital functions. TAS2Rs, a subfamily of GPCRs, play a primary role in recognizing bitter molecules and triggering events leading to the perception of bitterness, a crucial defence mechanism against spoiled or poisonous food. Beyond taste, TAS2Rs function is associated with many diseases as they are expressed in several extra-oral tissues. Given that the precise functioning mechanisms of TAS2R remain poorly understood, this study employed molecular dynamics simulations combined with network-based analysis to investigate local conformational changes and global structural correlations in different states of the receptor. The focus was on the human TAS2R46 bitter taste receptor, recently resolved experimentally, both in the presence and absence of strychnine, a known bitter agonist. The results showed that the ligand-bound state of the receptor exhibited more correlated dynamics compared to the apo state, and the presence of the agonist mediated the allosteric network between two helices (TM3 and TM6) which mainly convey the signal transferring from the extracellular to the intracellular region. By elucidating the hallmarks of the conformational changes and allosteric network of TAS2R46 under varying conditions, this study has enabled the identification of the unique structural and dynamics features of this receptor, thereby establishing a foundation for a more profound characterisation of this intriguing class of receptors.

1 Introduction

The sense of taste, also known as gustation, is crucial for mammals in evaluating the taste and composition of foods (Roper, 2017). Bitter, along with sweet, sour, salty, and umami, is one of the five basic taste modalities and allows to distinguish toxic molecules, providing the last checkpoint before the ingestion of potentially harmful substances. Bitter taste arises from the interaction of organic bitter molecules with type 2 taste receptors (TAS2Rs), which is a subfamily of G protein-coupled receptors (GPCRs) (Chandrashekar et al., 2000). These taste receptors are mainly expressed on the functional gustatory transduction units, i.e., the taste buds of the tongue, which are contained in gustatory papillae. Bitter compounds, binding to TAS2Rs, induce receptor conformational changes and initiate a downstream cascade of events inside the cell typical of GPCR signalling pathways, which ultimately leads to bitter taste perception. It is worth noticing that TAS2Rs are not only located in the taste buds of the tongue, but also several extra-oral tissues express them, such as heart, skeletal and smooth muscle, upper and lower airways, gut, adipose tissue, brain, and immune cells (Behrens and Meyerhof, 2013; Lee et al., 2019). Therefore, their function is not limited to taste evaluation, but it is associated with immune response and various diseases, such as obesity, diabetes, asthma, and cancer (Dotson et al., 2008; Liggett, 2014; Workman et al., 2015; Lu et al., 2017; Pan et al., 2017; Lee et al., 2019; Behrens and Lang, 2022). Therefore, extra-oral TAS2Rs could also represent a promising target for pharmacological intervention for specific diseases or health conditions. In this scenario, the understanding of the molecular mechanisms driving TAS2R functions is not limited to the taste perception field but can also improve our knowledge of pathologies and relative treatments.

From a structural point of view, TAS2Rs include a short extracellular N-terminus domain, an intracellular C-terminus domain, and seven transmembrane α-helixes (7TMs), which are connected by three extracellular loops (ECLs) and three intracellular loops (ICLs) (Zhang et al., 2017). These receptors present an orthosteric binding pocket which is in the EC part of the 7TMs bundle, involving the extracellular region of TMs II, III, V, VI, VII (Behrens and Meyerhof, 2009; Brockhoff et al., 2010; Behrens and Ziegler, 2020; Pallante et al., 2021). The secondary structure of these receptors is composed mainly of alpha-helix associated with the transmembrane bundle (about 70%–75%), whereas there are about 20% of bend, coil, and turn and, in some cases, a minor part of beta-sheets (1%–2%), composing the EC and IC domains. The structural similarity between TAS2R receptors and class A GPCRs has in the past led to their classification within the same family. Nevertheless, several recent investigations have pointed out specific characteristics of the bitter receptors, indicating that TAS2Rs can form a distinct family within the GPCRs (Di Pizio et al., 2016; Tokmakova et al., 2023). Indeed, the TM similarity between TAS2Rs and GPCRs is lower than 30% and most of the important conserved motifs of class A GPCRs (DRY motif in TM3, CWxP in TM6, and NPxxY motif in TM7) are missing in TAS2Rs. In TAS2Rs, the DRY motif observed in TM3 of class A GPCRs appears to be substituted by a highly preserved FYxxK motif, while the NPxxY motif in TM7 is replaced by HSxxL (Di Pizio et al., 2016). Moreover, recent literature indicated that TAS2Rs have unique conserved motifs, such as the TM1-2-7 interaction (E1.42-R2.50-S7.47 for TAS2R14), differing from those found in other class A GPCRs, such as the highly conserved N1.50-D2.50-N7.49, with the latter in the NPxxY motif. TM residues are denoted throughout the text using a superscript numbering system based on the Ballesterose-Weinstein (BW) method (Ballesteros and Weinstein, 1995). In this method, the residue corresponding to the most conserved residue in TM X of class A GPCRs is designated as X.50, and subsequent residues are numbered in relation to this position. Moreover, the typical “ionic lock”, which stabilises the inactive state and involves a salt bridge interaction in class A GPCRs, is replaced by a weaker hydrogen bond (HB) between Y3.50 and R6.36 in TAS2R14. Another important difference is the variation of an important residue related to activation which is W6.48 and Y6.48 for class A GPCRs and TAS2Rs, respectively (Tokmakova et al., 2023). The comparison between class A GPCRs and TAS2Rs is particularly challenging since no experimental structures of the bitter taste receptors have been determined until recently when cryo-electron microscopy structures of human TAS2R46 in both strychnine-bound and apo states were revealed (Xu et al., 2022). In particular, strychnine is a toxic bitter alkaloid known to be one of the main agonists that activate the TAS2R46-G-protein pathway (Brockhoff et al., 2007).

The understanding of the activation mechanisms of GPCRs is also still unclear. This process is an allosteric activation since it transduces various extracellular (EC) stimuli into the intracellular (IC) phenomena: the activation is induced by agonist binding and the subsequent G-protein recruitment. In the last decades, regarding class A GPCRs, a plethora of scientific publications have analysed and discussed GPCR activation mechanisms and allosteric networks (Venkatakrishnan et al., 2013; Venkatakrishnan et al., 2016; Bhattacharya and Vaidehi, 2014; Nivedha et al., 2018; Zhou et al., 2019; Bock and Bermudez, 2021). It is commonly recognised that the outward movement of the transmembrane helix 6 (TM6), which allows binding of the C-terminal part of the G-protein α subunit, is a common feature of GPCR activation triggered by ligand binding in the orthosteric site (Venkatakrishnan et al., 2016; Zhou et al., 2019; Bock and Bermudez, 2021). More in detail, the common activation mechanisms of class A GPCRs from an inactive to an active state involve the elimination of TM3-TM6 contacts, the formation of TM3-TM7 contacts, and the rearrangement of TM5-TM6. Globally, upon binding various agonists, this process triggers the outward movement of the cytoplasmic end of TM6 and the inward movement of TM7 toward TM3 (Nygaard et al., 2009; Rasmussen et al., 2011; Venkatakrishnan et al., 2016; Zhou et al., 2019). TAS2Rs seem to have peculiar features characterising their conformation and allosteric activation, presenting remarkable differences from class A GPCRs. The typical conformational changes related to class A GPCRs were not observed in the recently solved structures, thus leading to the classification of TAS2Rs as class T GPCRs (Xu et al., 2022). Furthermore, the only large conformational change detected was the different localization of the ECL2, which occupied the orthosteric binding pocket in the apo state. Finally, a rotation of the Y241^6.48 side chain toward the centre of the 7TMs bundle was observed in the ligand-bound state. Moreover, previous literature suggested that distinctive features of TAS2Rs compared to GPCRs comprise a higher involvement of ECL2 in the binding of agonists and the absence of the ECL2-TM3 disulfide bridge, which is involved in GPCR stabilization. This implies an alternative mechanism for regulating conformational states in TAS2Rs, potentially resulting in a less stabilized inactive state (Di Pizio et al., 2016). Nevertheless, the mechanism underlying the activation of TAS2R receptors remains incompletely elucidated, necessitating further investigation to delineate the distinctive features of these receptors in comparison to class A GPCRs. Moreover, in the absence of molecular structures, it has thus far been infeasible to gain insight into the conformational dynamics of TAS2Rs, thereby hindering an understanding of the principal dynamic molecular mechanisms and structural features underlying receptor function.

In light of the scientific context described above and based on the recent experimental structures of the TAS2R46 receptor (Xu et al., 2022), this study aims to explore the conformational changes and characterise the allosteric networks of this bitter taste receptor induced by the presence or absence of an agonist in the ligand-binding pocket through computational modelling methodologies. Computational techniques are rapidly becoming paramount techniques to highlight the key molecular features linked to the recognition of small molecules by taste receptors, especially for bitter taste, and identify subsequent molecular events and activation mechanisms (Born et al., 2013; Levit et al., 2014; Di Pizio and Niv, 2015; Di Pizio et al., 2017; Nowak et al., 2018; Dagan-Wiener et al., 2019; DiPizio et al., 2020; Pallante et al., 2021; Fierro et al., 2022; Malavolta et al., 2022). Moreover, the analysis of MD simulations through graph-based approaches is becoming an elective method to study the intra-protein structural communication and crucial residues for protein functions (Fanelli et al., 2016; Melo et al., 2020; Cannariato et al., 2023), as in the case of GPCRs (Siemers et al., 2019; Bertalan et al., 2020; Bondar, 2022). Using MD simulations coupled with network-based techniques, the present study aims to characterize the conformational states of the human TAS2R46 bitter taste receptor and investigate the molecular mechanisms at the basis of the allosteric mechanical communication from the EC to the IC regions. We simulated different receptor conditions, i.e., strychnine-bound, ligand-free, and transition states, to pinpoint major differences in the receptor dynamics and allosteric networks. Consequently, this work, through in silico simulations of the dynamic behaviour of the TAS2R46 receptor, offers new insights into the molecular structural mechanisms governing the function of this receptor and establishes the foundation for a comprehensive understanding of the unique functioning features of TAS2Rs.

2 Materials and methods

2.1 System setup

The molecular structures of human TAS2R46 were retrieved from the RCSB Protein Data Bank, using the PDB codes 7XP6 and 7XP4 for strychnine-bound and Apo states, respectively (Xu et al., 2022). Only the receptor and the ligand (if present) were preserved from the original PDB files. The missing residues (157–172) of the experimental strychnine-bound TAS2R46 structure were modelled using the corresponding model from the AlphaFold Protein Structure Database (P59540 entry) (Jumper et al., 2021) after root-mean-square fitting on the alpha carbons of the experimental structure. Instead, missing residues in disordered regions at the N-terminus (residues 1–2) and C-terminus (residues 302–309) were not included in the present study, as they are absent in the experimental structures and exhibit low pLDDT scores (below 70) in the corresponding AlphaFold model. The stereochemical quality of the obtained model was validated using the PROCHECK software (Laskowski et al., 1993), showing that more than 90% of the residues fall within the most favoured regions of the Ramachandran plot (see Supplementary Figure S1). Additionally, the TM-score (Zhang and Skolnick, 2004; Xu and Zhang, 2010) between the obtained model and the original experimental structure is 0.97, indicating a high degree of similarity between the structures. We also double-checked the obtained structure by comparing it with models generated using alternative procedures. Specifically, we created new models of TAS2R46 through both SWISS-MODEL (Waterhouse et al., 2018) and MODELLER (Webb and Sali, 2016), and then compared them with our original structure to assess any remarkable differences. For the SWISS-MODEL, we followed a standard template-based approach, using the AlphaFold model (code: P59540) as the template. For MODELLER, we employed a multi-template modelling strategy, utilizing the Cryo-EM structure (PDB: 7XP6) as the primary template and integrating the AlphaFold model (code: P59540) to reconstruct the missing residues (157–172). The obtained models exhibited only minor differences in their structural arrangement, and our model demonstrated acceptable stereochemical properties in comparison to those generated by the other modelling methods (see Supplementary Figure S1).

The structure of strychnine was refined using MOE (Molecular Operating Environment (MOE), 2022.02 Chemical Computing Group ULC, 1010 Sherbooke St. West, Suite #910, Montreal, QC, Canada, H3A 2R7, 2023., 2022), predicting the protonation state at neutral pH. As reported in previous literature (Xu et al., 2022), the tertiary amine of the molecule is protonated and its total charge is +1, as also reported in the DrugBank database for physiological conditions (https://go.drugbank.com/drugs/DB15954).

Therefore, two structures were obtained: (i) the receptor bound to strychnine (Holo) and (ii) in the absence of the ligand (Apo). A third structure, which will be referred to as Trans, was then defined by removing strychnine from the binding pocket of the Holo state.

For each of the three structures, the protein-membrane complex has been built using CHARMM-GUI (Feng et al., 2023) as follows. The homogeneous bilayer membrane was composed of phosphatidylcholine (POPC, 16:0/18:1 acyl chains), as reported in previous literature (Hénin et al., 2006; Zou et al., 2019). The system has been inserted in a rectangular box with dimensions of 8 × 8 × 11 nm³. This choice allowed us to reach a good compromise between acceptable computational cost and the need to ensure the minimum image convention. Then, the system was solvated by using the TIP3P water model before adding an appropriate number of Na⁺ and Cl⁻ ions to reach a physiological salt concentration of 0.15 M and neutralize the overall system charge. The AMBER19SB force field (Lee et al., 2020; Tian et al., 2020) was used to describe the protein, ions, and water, the Lipid-21 forcefield (Dickson et al., 2022) was used for lipids, and the General Amber Force Field (GAFF2) forcefield (Wang et al., 2004) to obtain the topology for strychnine. System preparation and topology definition were performed directly in CHARMM-GUI as done in previous literature (Born et al., 2013; Sengupta et al., 2018; Zhou et al., 2019; Zou et al., 2019; Pallante et al., 2023; 2024). The protein-membrane complexes in the Holo, Trans, and Apo states are shown in Figure 1.

Figure 1

Figure 1. Visual representation of the Holo, Trans and Apo systems. The 7 TMs of TAS2R46 are highlighted with different colours, the hydrophobic tails and hydrophilic head of POPC are coloured in light green and red, respectively.

2.2 Molecular dynamics simulation

MD simulations were performed using the simulation engine GROMACS 2022 (Bauer et al., 2023) starting from the models described above. For each system, the same simulation protocol was followed as detailed below. First, energy minimization was performed through the steepest descent method for 5,000 steps. Then, three replicas were performed following the CHARMM-GUI protocol of equilibration before a production phase of 500 ns long. In detail, six equilibration steps were performed gradually reducing the position restraints on the lipids and protein-heavy atoms (from 1,000 to 0 kJmol^-1 nm^-1 for lipids, from 4,000 to 50 kJmol^-1nm^-1 for the protein backbone and from 2,000 to 0 kJmol^-1nm^-1 for protein side-chain heavy atoms). The equilibration phase started with two NVT simulations, performed for 125 ps with a timestep of 1 fs. Then, four NPT equilibration steps were performed, the first one for 125 ps with a timestep of 1 fs, the remaining three for 500 ps with a timestep of 2 fs. The total time of equilibration of each replica was 1.875 ns. The NVT simulations were performed at a reference temperature of T = 303,15 K (τ = 1 ps), which is above the phase-transition temperature for POPC, using the Berendsen thermostat (Berendsen et al., 1984), while NPT simulations were carried out at 1.0 bar using the Berendsen barostat with semi-isotropic coupling (τ = 5 ps). Finally, the unrestrained production phase was carried out in the NPT ensemble with a Nose-Hoover thermostat and Parrinello-Rahman barostat for 500 ns. The leapfrog integrator was used, using a time step of 2 fs. The PME algorithm was used for electrostatic interactions with a cut-off of 0.9 nm. A reciprocal grid of 72 × 72 × 96 cells was used with 4th-order B-spline interpolation. A single cut-off of 0.9 nm was used for Van der Waals interactions. LINCS (LINear Constraint Solver) algorithm for h-bonds (Hess et al., 1997) was applied in each simulation step. Three simulation replicas were performed for each investigated system to increase the statistics of the data and ensure the repeatability of the results. Therefore, a total of 4.5 µs of simulation during the unrestrained simulation phase were performed.

2.3 Analysis

2.3.1 Conformational analysis

The conformational changes of the investigated states (Holo, Trans, Apo) were assessed by calculating the root-mean-squared deviation (RMSD) from the initial configuration of backbone atom positions throughout the trajectory. The root-mean-squared fluctuation (RMSF) of alpha carbons was computed in the last 400 ns of each replica, allowing the evaluation of the fluctuations of each residue during the simulation time. The final 400 ns of simulation were concatenated obtaining a single 1.2 µs trajectory for each system under investigation. Then, the following analysis was performed considering a sampling time of 50 ps, unless otherwise specified.

To further examine the most frequent conformational states, cluster analysis was performed. Specifically, we analysed the concatenated final 400 ns of each replica, sampled at a 1 ns time step, using the linkage clustering algorithm implemented in GROMACS. The analysis was performed with an RMSD cutoff of 0.1 nm, considering only the backbone atoms. The selected cutoff threshold accounts for the substantial fluctuations typically observed in GPCRs within a single conformational state and has been chosen according to previous literature (Latorraca et al., 2017).

The Protein-Ligand Interaction Profiler (PLIP) (Salentin et al., 2015) tool was used to evaluate the specific interactions between strychnine and TAS2R46, underlining the most important residues involved in strychnine binding and the types of interactions established (i.e., hydrogen bonds, hydrophobic interactions, salt bridges, etc.). In detail, the probability of a specific interaction between the receptor and strychnine has been evaluated by considering the interactions on each frame and then averaging the number of occurrences of the interactions on the total number of frames, as done previously (Miceli et al., 2022).

The binding pocket volume was evaluated using the Epock tool (Laurent et al., 2015). The maximum englobing region (MER), i.e., the region of space delimiting the binding pocket, was defined as a sphere of radius 1.3 nm located at the centre of mass of strychnine, considering the starting configuration of Holo replicas. Before the binding pocket volume evaluation, the concatenated trajectories were RMS fitted on the configuration used to define the MER. Moreover, in the analysis, the residues 153 to 176 and, in the Holo system, strychnine were not considered to allow a better comparison between the three systems under study.

Then, the conformation of Y241^6.48, which is characterized by a different position in the experimental structures of human TAS2R46, was assessed by calculating its dihedral angle, defined by the atoms CA-CB-CG-CD1. Additionally, the interactions and relative distances of Y241^6.48 with other residues in the protein were evaluated using the PLIP tool for interaction analysis and GROMACS pairdist functionality for minimum distance measurements.

The calculation of RMSD, RMSF, and cluster analysis was performed using GROMACS, whereas the secondary structure was evaluated through the STRIDE (Heinig and Frishman, 2004) software package. Conformational analysis was performed through custom-made scripts in Python using the MDAnalysis module. All plots were generated using the Matplotlib (Hunter, 2007) and Seaborn (Waskom, 2021) libraries, whereas the three-dimensional representations of receptor structures were rendered in Visual Molecular Dynamics (VMD) software (Humphrey et al., 1996).

2.3.2 Generalized correlation analysis

A correlation analysis was also employed to identify the regions of the receptor more correlated with each other in the Holo, Trans, and Apo states. In particular, the generalized correlation coefficient (r_MI) was computed as it takes into account linear and non-linear contributions to correlations (Lange and Grubmüller, 2005). This analysis was performed since correlated motions are essential for the biomolecular function of several systems, such as orthosteric and structural signal transduction in GPCRs (Scheer and Cotecchia, 1997).

The generalized correlation coefficient between residues i and j has been computed as:

$r_{MI}\left[i,j\right]={\left(1-e^{-\frac{2}{3}I\left[i,j\right]}\right)}^{\frac{1}{2}}$

where $I\left[i,j\right]$ is the mutual information between the positions of residues i and j, where the position of one residue was defined as the position of its alpha carbon. The mutual information has been computed using the density estimator described by Kraskov et al. (Kraskov et al., 2004) with neighbour parameter k of 6 as done in previous literature (Lange and Grubmüller, 2005; Melo et al., 2020):

$I\left[i,j\right]=\psi \left(k\right)-\frac{1}{k}-〈\psi \left(n_i\right)+\psi \left(n_j\right)〉+\psi \left(N\right)$

where N is the total number of simulation frames, ψ(x) is the digamma function, n_i is the number of frames in which residue i is close to the one in the reference, and $〈〉$ stands for the average by varying the reference frame over the trajectory. For the generalized correlation coefficient calculation, a trajectory sampling time of 1 ns was considered to avoid considering correlated frames in the analysis as done in previous literature (Cannariato et al., 2023; Manrique et al., 2023).

2.3.3 Dynamic network analysis

To investigate structural communication within the TAS2R46 receptor, the MD simulations were analyzed through the Dynamical Network Analysis approach, using the dynetan library (Melo et al., 2020). This specific network analysis was selected as it is based on the generalized correlation coefficient, thus considering nonlinear contributions to amino acid dynamical correlations. In detail, in the Dynamical Network Analysis, each protein residue was represented by a node located in its alpha carbon and strychnine was modelled through a single node in the closest atom to its centre of mass. Nodes were linked with edges if the frequency of contacts between the corresponding residues was greater or equal to 0.75, considering two residues to be in contact at a simulation frame if the shortest distance between their heavy atoms is lower than 0.45 nm. These thresholds were chosen based on the original development of dynamic network analysis in protein systems (Melo et al., 2020). The edges of the obtained graph were weighted using the generalized correlation coefficient (Lange and Grubmüller, 2005), computed using a sampling time of 1 ns (Manrique et al., 2023). The network was characterized in terms of betweenness centrality and eigenvector centrality. The betweenness centrality of an edge or node is described as the fraction of the shortest paths in which the considered edge or node is involved and highlights the importance of edges and nodes for the connection of distant parts of the network. In this study, the shortest path between two nodes was defined as the path that maximizes the sum of correlations between the nodes involved in the path. The betweenness centrality is computed as:

$b_i=\frac{1}{C}{\displaystyle \sum _{s,t\in V}}\frac{\sigma \left(s,t|i\right)}{\sigma \left(s,t\right)}$

where $\sigma \left(s,t\right)$ is the number of shortest paths between nodes s and t, $\sigma \left(s,t|e\right)$ is the number of such paths passing through node e, V is the ensemble of graph nodes, and C is a normalization factor to allow the comparison of networks with different numbers of nodes. In particular, for a graph of n nodes, is equal to

$\frac{2}{\left(n-1\right)\left(n-2\right)}\text{\hspace{0.17em}or\hspace{0.17em}}\frac{2}{n\left(n-1\right)}$

for node and edge betweennesses, respectively.

The eigenvector centrality of a network measures the influence of a node in a graph considering the topology of the graph itself, such that the centrality of a node depends on the centrality of its neighbours. Therefore, if a node has many connections with nodes of small influence, it will also have a small centrality in the network. On the other hand, if a node has few connections with nodes very influential in the network, it might have high centrality because of its indirect influence within the graph. Mathematically, the eigenvector centrality of node i is the i-th entry of the eigenvector (x) of the adjacency matrix weighted on the correlation between the nodes (A):

$Ax=\lambda x$

3 Results

From the visual inspection of the RMSD plots, the conformational states appear to deviate from the starting structure during the first 100 ns of simulation (Supplementary Figure S2). Following this period, the RMSD values oscillate around their mean values with standard deviations below 0.04 nm, providing a quantification of the exploration of the investigated conformational state. Moreover, the receptor structure was analyzed in terms of secondary structure probability. This analysis highlighted that the TMs maintained their structure across the simulations without any remarkable alteration of their secondary structure regardless of the receptor’s configuration state (Supplementary Figure S1).

3.1 Conformational changes of TAS2R46 and structural effects of strychnine

This paragraph describes the conformational analysis conducted to identify remarkable differences among the three analyzed states - Holo, Trans, and Apo. The objective of these analyses is to identify specific structural characteristics associated with the different states of TAS2R46.

The Holo, Trans, and Apo states were initially characterized in terms of RMSF, which was computed to evaluate the flexibility of the different receptor regions. As expected, the most flexible regions of TAS2R46 were the unstructured ones and some differences have been observed between the three receptor states for ICL3 and ECL3 (Supplementary Figure S4). The Trans state is characterized by higher fluctuations in the ECL3, but the flexibility of the ICL3 in this state is the same as the Holo state. Interestingly, the ICL3 region is more stable and less flexible in the Apo state compared to the other two. On the other side, the ECL2 displays similar fluctuations in the three systems, although it is characterized by a different conformation in the Apo state compared to the Holo and Trans states (Figure 1). The most frequent structural rearrangements of the three systems were analyzed through a cluster analysis (see Materials and Methods for details). This analysis revealed 22 clusters for the Holo state and 12 clusters for the Apo state, while the Trans system exhibited 105 clusters, indicating greater structural fluctuations compared to the other states (see Supplementary Figure S5). For further details regarding the cluster analysis results, please refer to the Supplementary Material.

Then, the interaction between strychnine and the receptor was analyzed in terms of type and stability using PLIP. This analysis revealed three main interactions with a probability greater than 0.5, namely, two hydrophobic interactions with residues Y85^3.29 and W88^3.32 and a salt bridge interaction with E265^7.39 (Supplementary Figure S6). This confirms that the main interactions detected in the experimental structure (Xu et al., 2022) are conserved and involve the TM3 and TM7, highlighting the stability of the ligand inside the receptor’s binding pocket during the simulations.

Previous literature reported a rotation of the Y241^6.48 side chain from pointing outward the 7TMs bundle to pointing into the core of TAS2R46 in the ligand-bound state (Xu et al., 2022). Therefore, the rotation of the Y241^6.48 side chain was evaluated using the dihedral angle defined by the atoms CA-CB-CG-CD1 (Figure 2A, see also Supplementary Video S1). The results revealed higher absolute values of the dihedral angles in the Holo state compared to the Apo state, indicating that the residue is pointing towards the centre of the TM bundle. Interestingly, the dihedral angle for the Trans system was characterized by a wider distribution compared to the other states. The position of residue Y241^6.48 was also characterized in terms of its interactions and distances with the other residues of the protein. Specifically, the interactions were evaluated using PLIP over the concatenated 1.2 μs trajectory (last 400 ns of the three replicas) for each system. The interaction probability between the residues was determined by calculating the number of frames in which they formed an interaction, divided by the total number of frames in the trajectory. We considered only interactions for residues characterised by probability higher than 0.8 in at least one of the three systems (Holo, Trans, Apo). All these interactions are hydrogen bonds, with the most notable differences between the systems being the interactions of Y241^6.48 with residues N92^3.36 and T274^7.48 (see Supplementary Figure S7). More in detail, Y241^6.48 and N92^3.36 formed a hydrogen bond primarily in the Holo state, where the interaction was present in 94.7% of the frames. In the Trans state, the probability of interaction decreased to 31.1%, and no hydrogen bond was observed in the Apo state (0% probability). Conversely, the hydrogen bond between Y241^6.48 and T274^7.48 is found only in the Apo system (probability to 85.9%), while no interactions were detected in the Holo or Trans states. Interestingly, all systems demonstrated hydrogen bonds with a probability higher than 85% with residues L237^6.44, S244^6.51 and I245^6.52. The minimal distances of Y241^6.48 from residues N92^3.36 and T274^7.48 over the concatenated 1.2 μs trajectory were also evaluated (Figure 2B). In agreement with the interaction probability analysis, residues Y241^6.48 and N92^3.36 are in close contact in the Holo state (0.21 ± 0.05 nm), while demonstrating higher distances in the Apo system (0.98 ± 0.06 nm). The Trans state, instead, showed a wider distribution of distances with the main peak between the Apo and Holo states. Conversely, the distances between residues Y241^6.48 and T274^7.48 assumed lower values in the Apo state (0.20 ± 0.04 nm) compared to the Holo (0.76 ± 0.05 nm) and Trans (0.73 ± 0.04 nm) states, which showed similar values. However, the different conformations of Y241^6.48 are not related to outward movements of the IC region of TM6, which is one of the hallmarks of the class A GPCRs activation process (Zhou et al., 2019): the distance between the IC areas of TM3 and TM6 does not show remarkable differences between the three states (see Supplementary Figure S8), strengthening the previously observed different behaviour compared to class A GPCRs (Xu et al., 2022).

Figure 2

Figure 2. Conformational analysis of TAS2R46 in the three investigated states. (A) Probability distribution of the dihedral angle (defined by the atoms CA-CB-CG-CD1) of residue Y241^6.48. Orange and blue dashed lines indicate the dihedral angles observed in the experimental structures of the Holo (PDB: 7XP6) and Apo states (PDB: 7XP4), respectively. (B) Probability distribution of minimal distances of Y241^6.48 from residues T274^7.48 (left) and N92^3.36 (right). (C) Rendering of representative structures of the three receptor states, highlighting the relative positions of residues Y241^6.48 (white), T274^7.48 (red) and N92^3.36 (pink). Hydrogen bonds are depicted as black dotted lines. The Holo structure is shown in orange, Trans in green, and Apo in blue.

Then, the volume of the TAS2R46 orthosteric binding pocket was evaluated in the three systems. It is worth mentioning that, in this analysis, the ECL2 was not considered to prevent any alterations in the estimation, as it is located within the investigated pocket in the Apo state. In this way, we evaluated only the change in the volume of the orthosteric binding pocket due to the rearrangements of the helices forming the ligand pocket, regardless of the initial location of the ECL2, which varies in the ligand-free or bound states. The results showed that the volume of the orthosteric pocket is higher in the absence of strychnine (Figure 3). Interestingly, the Trans state displayed a wide volume distribution with the main peak between the Apo and Holo states.

Figure 3

Figure 3. (A) Visual representation of the binding pocket volume (in red) of TAS2R46 for one selected frame of the Apo system. (B) Probability distribution of the binding pocket volume for the Holo, Trans, and Apo states, represented in orange, green and blue, respectively. The ECL2, not shown in this representation, was excluded from the calculation of the binding pocket volume to prevent any alterations in the estimation, as it is located within the investigated pocket in the Apo state.

3.2 The ligand-bound state is characterized by an increased dynamical correlation

In this section, we assessed intra-structure correlations and examined how they were altered by the presence or absence of strychnine to identify major differences. The correlation between the different receptor regions in the three states was analyzed using the generalized correlation coefficients as described in the Material and Methods section. We evaluated the correlation between residues for the Holo, Trans, and Apo states, which showed remarkable differences in their dynamic behaviour. In the Holo system, high correlations were observed between the EC region of the receptor, whereas the IC regions of TM5 and especially TM6 were less correlated with the rest of the 7TM bundle (Figure 4A). Interestingly, the ICL3 displayed higher correlations than the neighbouring regions of TM5 and TM6, especially with the ECL1. Our results also showed that the removal of strychnine from the orthosteric binding pocket induced a remarkable and general loss of dynamic correlation (more bluish areas) within the 7TM bundle (Figure 4B). High correlations for the Trans state were observed only for the EC regions of TM1, TM6, and TM7, the ECL2, and, less markedly, the ECL1. On the other hand, the regions showing the lowest correlation with the 7TM bundle were the IC regions of TM3, TM6, and TM7. Finally, in the Apo state, the results highlighted remarkable correlations of the ECL2 with the EC regions of the receptor, while lower values could be observed for the IC regions (Figure 4C). To easily compare the receptor states, we also calculated the differences of the generalized correlation coefficient between the Trans and Holo states, i.e., Δr_MI(Trans, Holo) = r_MI(Trans) - r_MI(Holo), and between the Apo and Holo states, i.e., Δr_MI(Apo, Holo) = r_MI(Apo) - r_MI(Holo) (Figure 4D). Comparing the observed correlations for the three states, it was observed that, in general, the receptor evolved dynamically in a more correlated way in the presence of strychnine, whereas its absence was associated with the decorrelation of the IC regions. Moreover, it is worth noticing that the Trans state overall demonstrated lower structural correlation values compared to the other two. Furthermore, Figure 4D highlighted that the ICL3 and the IC portion of TM3 are less correlated with the rest of the 7TM bundle in Trans and Apo states than in the Holo system. Interestingly, for both Trans and Apo systems, the IC portion of TM6 is more correlated with the rest of the receptor compared to the Holo state. Therefore, the Apo and Trans systems behaved in the same way if compared to the Holo state.

Figure 4

Figure 4. Matrices showing the generalized correlation coefficient between the residues of TAS2R46 in the (A) Holo, (B) Trans, and (C) Apo states. In panel (D), the differences of the generalized correlation coefficients are shown: Δr_MI(Trans, Holo) = r_MI(Trans) - r_MI(Holo) in the lower triangle and Δr_MI(Apo, Holo) = r_MI(Apo) - r_MI(Holo) in the upper triangle. Matrices on panels from (A to B) are coloured according to the left colorbar, while the matrix on panel (D) is colored according to the right colorbar.

3.3 Holo structural network is mediated by the TM3-TM6 connection

As a first step in the analysis of the dynamic networks of the Holo, Trans and Apo systems, we extracted information regarding node and edge centralities from the graph topology: after plotting the node betweenness centrality (NBC) and edge betweenness centrality (EBC) versus the node and edge rank, respectively, we observed that both NBC and EBC were similar in the three networks (Supplementary Figure S9). We then identified the knees of these curves for the three systems, and used the lowest values, related to the Holo state, as thresholds. Only nodes or edges with NBC or EBC higher than the identified threshold were considered in the subsequent analysis. Firstly, we considered the distribution of NBC in the different helices of the receptor with specific attention to the upper tails of the distributions and outliers: these highlight residues of particular importance in the information flow inside the receptor. We hereby observed that for the Holo state the tails of TM3 and TM6 extended toward higher centralities, while other isolated residues were more central in other helices like TM1. On the other hand, for the Trans and Apo states, residues belonging to the TM3 were in general higher than the TM6 ones (p < 0.05 with Wilcoxon-signed-rank test), especially in the case of the Trans system (Figure 5A). Moreover, we observed how the tail of the distribution for TM4 extends toward higher values for Holo and Trans states compared to the Apo state. Lastly, concerning the regions TM2 and TM5, we noted how their NBC distributions were similar for the Holo and Apo systems, but differed in the Trans system, whose tails showed higher centralities. Considering that the NBC measures how influential a node is in the flow of information inside the network, these results highlight that in TAS2R46 such flow is mainly conveyed by the TM3 in the absence of strychnine, as opposed to both TM3 and TM6 in the presence of the ligand. Similar information can be obtained from the visual inspection of the networks after highlighting the EBC (Supplementary Figure S10).

Figure 5

Figure 5. Dynamical Network Analysis of the Holo, Trans, and Apo systems. (A) Boxplots representing the distribution of node betweenness centrality, with nodes grouped into the respective TMs. Data emphasize that in the absence of the ligand TM3 assumed a primary role in information transfer, whereas in the presence of strychnine TM6 nodes also demonstrated high values of betweenness centrality. (B) Eigenvector centrality of TAS2R46 nodes. The TMs and loops are highlighted as well as residues pinpointed in previous experimental analysis. Results showed higher centrality of TM6 with strychnine, especially for H205^5.68, H224^6.31, and Y241^6.48, and decreased centralities in the Trans state. (C) Visual representation of the dynamic networks, where the receptor is viewed from the IC region. The edges of the network are red cylinders with a radius proportional to the betweenness centrality of the edge.

As a second step, we analyzed the influence of the single residues on the overall network in terms of eigenvector centrality (EVC): in this case, we observed similar values of EVC for the EC and IC portions of the TM3 in the Apo and Holo systems, while the TM6 was characterized by remarkably higher EVC in the presence of strychnine (Figure 5B). In particular, on one hand, the EVCs of W88^3.32, involved in strychnine binding, and of Y106^3.50, involved in G-protein interaction at the TM3 level (Xu et al., 2022), were similar in the Holo and Apo states. On the other hand, H205^5.68, H224^6.31, and Y241^6.48 were more central in the presence of the ligand. Moreover, except for W88^3.32, all mentioned residues showed lower centralities in the Trans systems. Lastly, we focused on the difference, between the three systems in terms of connections between TMs. Interestingly, the main difference appears related to the TM3 (Supplementary Figure S11): in the presence of the ligand, there is a connection between TM3 and TM6 mediated by the edge between N92^3.36 and Y241^6.48, whereas in the ligand-free state TM3 is connected to TM7 through the N96^3.40-Y271^7.45 edge, and in the Trans state there is no connection between TM3 and TM6 or TM7 (Figure 5C, see also Supplementary Figure S11).

We finally investigated the effects of the highlighted changes in the dynamic network properties on the structural communication inside TAS2R46. To this end, we computed the optimal paths linking strychnine binding residues, in particular W88^3.32 and E265^7.39, to H224^6.31 and Y106^3.50, which were identified as G-protein binding residues by (Xu et al., 2022). Coherently with the network topology analysis, TM3 drove the communication between W88^3.32 and H224^6.31 in all systems. Conversely, we observed differences when considering the paths from W88^3.32 to H224^6.31 and from E265^7.39 to Y106^3.50 (Figure 6, see also Supplementary Table S1): in greater detail, in the presence of strychnine, both paths involve an edge between N92^3.36 and Y241^6.48, which creates a bridge between TM3 and TM6 conveying information between the EC and IC regions. However, in the Apo state, the path from W88^3.32 to H224^6.31 reaches the IC region of the receptor mainly through the TM3, then passes to the TM5, and finally to H224^6.31. Interestingly, we noted the same path in the Trans state, with the only difference being that the path reaches TM6 directly at the H224^6.31 residue, instead of A227^6.34 as in the Apo state. On the other hand, the connection between E265^7.39 and Y106^3.50 for the Apo system is mediated by the Y271^7.45-N96^3.40 edge linking TM3 and TM7, whereas in the Trans and Holo systems, the paths are mediated by the TM2 and TM6, respectively (Figure 6, see also Supplementary Table S1). Therefore, the change in the network topology is reflected in remarkable differences in the communication between the EC and IC regions.

Figure 6

Figure 6. Visual representation of the optimal path connecting strychnine and G protein-binding residues. The left column represents the communication path between W88^3.32 in TM3, highlighted as a green sphere, and H224^6.31 in TM6, represented as a cyan sphere. The right column represents the communication path between E265^7.39 in TM7, highlighted as a green sphere, and Y106^3.50 in TM3, represented as a cyan sphere. The edges in the paths are represented by cylinders whose radius is proportional to the generalized correlation coefficient between residues.

4 Discussion

TAS2Rs constitute the molecular basis of bitter taste perception. Concurrently, there is mounting evidence that various extra-oral tissues express bitter taste receptors, and their activation elicits diverse signals and cellular responses essential for metabolism and homeostasis (Dotson et al., 2008; Liggett, 2014; Pan et al., 2017). Extra-oral TAS2Rs are associated with several diseases, and they could represent a promising target for pharmacological intervention. However, the lack of experimental structures of TAS2Rs represented an important hurdle for understanding the mechanisms underlying bitter receptor activation. The recent release of the experimental structures of TAS2R46 both coupled with strychnine and in a ligand-free state (Xu et al., 2022) paves the way for the structural characterization of this receptor and for the definition of specific dynamic hallmarks underlying the receptor state. Strychnine is of important interest as a ligand because it is experimentally known to target not only TAS2R46 but also other bitter taste receptors, such as TAS2R10, and even other proteins (Jensen et al., 2006; Brockhoff et al., 2010; Born et al., 2013; Sandal et al., 2015; Xue et al., 2018). In this context, this work aims to contribute to the overall characterization of local conformational changes and global allosteric networks in experimentally determined structures of TAS2R46. By doing so, it seeks to identify unique structural and dynamic features of TAS2R46 that may be linked to its activation mechanisms. In detail, the conformational and dynamical features of the receptor in the presence and absence of strychnine were considered and then the dynamic network analysis was employed to investigate how such features could be related to the allosteric activation, which implies the transfer of information from the EC region to the IC region. Indeed, network analysis of protein dynamics has already proven their potential for analyzing the structural communication in macromolecular structures, including GPCRs (Fanelli et al., 2016; Bertalan et al., 2020; Melo et al., 2020; Sullivan et al., 2020; Bondar, 2022).

Initially, conformational features typical of class A GPCRs or highlighted by previous experimental studies were considered. Regarding the interaction between strychnine and TAS2R46, the PLIP analysis showed that the initial interactions with W88^3.32 and E265^7.39 remained stable during the MD simulation, and at the same time pinpointed a stable hydrophobic interaction with Y85^3.29 (Supplementary Figure S6). It is worth mentioning that these residues have been highlighted in previous literature to be involved in contacts with strychnine in TAS2R46 (Sandal et al., 2015; Xue et al., 2018) or agonist selectivity in other TAS2Rs, including TAS2R38 (Marchiori et al., 2013), TAS2R10 (Born et al., 2013), TAS2R16 (Sakurai et al., 2010), TAS2R43 (Pronin et al., 2004; Brockhoff et al., 2010), TAS2R44 (Pronin et al., 2004; Brockhoff et al., 2010), TAS2R47 (Pronin et al., 2004), TAS2R1 (Singh et al., 2011), and TAS2R4 (Pydi et al., 2012). Additionally, residue W88^3.32 is involved in the activations of TAS2R43 and TAS2R30 (Pronin et al., 2004), whereas residue E265^7.39 was previously reported to be involved in a salt-bridge interaction with strychnine (Xue et al., 2018) and its mutation implicates a reduced responsiveness to the compound (Brockhoff et al., 2010). However, while residue W88^3.32 is highly conserved among TAS2Rs (84%), residues E265 and Y85 show a lower conservation (32% and 24%, respectively), suggesting a different impact on receptor selectivity. Furthermore, residues in positions 3.32 and 7.39 are known to be involved in the binding of diverse ligands in class A GPCRs (Brockhoff et al., 2010; Venkatakrishnan et al., 2013).

Nevertheless, in our study, we were unable to identify the specific conformational activation hallmarks typically observed in class A GPCR. Specifically, our simulations did not reveal a remarkable outward movement of TM6 (Supplementary Table S8) or an increase in the A¹⁰⁰ index (Supplementary Table S12). This indicates that the activation mechanisms of TAS2Rs remain not fully understood and suggests that this class of receptors may possess unique structural characteristics that differentiate them from other GPCR families, as highlighted in previous literature (Ibrahim et al., 2019; Zhou et al., 2019; Bock and Bermudez, 2021). Therefore, additional evidence is needed to investigate the conformational activation hallmarks of TAS2Rs and to support the recent classification of TAS2Rs into the class T of GPCRs (Pándy-Szekeres et al., 2023), which is characterized by distinct features that differ from those of class A receptors (Di Pizio et al., 2016; Topin et al., 2021). On the contrary, as was observed for the class A GPCRs activation process (Dalton et al., 2015), in the TAS2R46 strychnine-bound state the volume of the orthosteric binding pocket is lower than the one in the Apo state and the removal of strychnine from the Holo system was associated with an increase in the pocket volume (Figure 3B). Notably, the Trans state exhibited a broad volume distribution, with its main peak positioned between those of the Apo and Holo states.

The overall dynamic of TAS2R46 was influenced by the presence or absence of strychnine in the orthosteric binding site. This was highlighted by the analysis of the intra-receptor correlations, which showed higher overall correlation values for the Holo state, while both the Trans and the Apo states were characterized by lower correlations (Figure 4). In particular, in Trans and Apo states, reduced correlation for the IC region of TM3 and for the ICL3, which are important regions for G-protein binding, was observed. Therefore, the Apo and Trans systems were characterized by similar behaviour in terms of intra-receptor correlations if compared to the Holo state. The main difference between the Apo and Trans states could be linked to the different conformation of the ECL2. Indeed, in the initial state, the ECL2 is placed within the orthosteric binding site in the Apo state and outside in the Trans state. During the MD simulation, the ECL2 showed remarkable correlations with other regions of the receptor only for the Apo structure. Furthermore, it is noteworthy that the removal of strychnine in the Trans state resulted in a globally less correlated state if compared to the Apo and Holo states. Therefore, the whole receptor in the Trans state seemed to evolve in time in a more decorrelated way.

The different conformational and dynamical features of the receptor in the three states were then related to different behaviours of the receptor in terms of structural communication using the dynamic network analysis. The betweenness centrality pointed out the importance of TM3 in the information transfer from the EC to the IC region of the receptor independently of its state (Figure 5). This finding is of particular interest given the pivotal role of TM3 in class A GPCRs, serving as a structural and functional hub responsible for maintaining the receptor scaffold in both active and inactive states (Venkatakrishnan et al., 2013). Thus, these results support the hypothesis that TM3 represents a central and critical structural element that defines the overall structure of both class A GPCRs and TAS2Rs in a similar manner. On the other hand, the TM6 was characterized by increased importance and influence inside the network in the presence of strychnine (Holo state) (Figures 5A, B). This information is particularly relevant as it highlights that high values of the centrality of TM6 in the presence of strychnine are related to the above-mentioned connection to TM3 through the network edge (and a hydrogen bond interaction) between Y241^6.48 and N92^3.36 (Figure 5C). Moreover, it is worth noting that the central role of TM6 in the activation process has also been highlighted in previous literature related to class A GPCRs (Trzaskowski et al., 2012; Venkatakrishnan et al., 2013; Zhou et al., 2019). Moreover, as shown in Figure 5C, residues H205^5.68, H224^6.31, and Y241^6.48 exhibited higher eigenvector centrality values in the Holo state. This is particularly noteworthy as these residues are involved in G-protein interactions at the TM5 and TM6 levels, respectively, and given the previously mentioned significance of Y241^6.48 in the activation of TAS2R46 (Xu et al., 2022). Regarding the network connecting the TMs, the dynamic network analysis revealed distinct correlations for the Holo, Trans, and Apo structures (Figure 4C and see also Supplementary Figure S11). Specifically, the Holo structure exhibited a high correlation between TM3 and TM6, while the Apo structure showed an edge connecting TM3 and TM7. Interestingly, the dynamic network of the Trans structure did not show any edges connecting TM3-TM6 or TM3-TM7. These results once more indicate a possible remarkable distinction between class A GPCRs and TAS2Rs. Specifically, following activation in class A receptors, the outward movement of the TM6 helix (not relevant for TAS2Rs as noted above) results in a reduction in contacts between the TM3 and TM6 helices and the formation of contacts between the TM3 and TM7 helices (Venkatakrishnan et al., 2016; Zhou et al., 2019), thus showing a different trend to that observed during the dynamics of TAS2R46.

We specifically focused our attention on the residue at position 6.48, which is considered crucial in the activation process of class A GPCRs (Tokmakova et al., 2023) and has been suggested to be involved in the “toggle switch” mechanism in previous work related to the experimental resolution of the TAS2R46 structure (Xu et al., 2022). To investigate side-chain rearrangements during the MD simulations, we assessed the dihedral angle of Y241^6.48 to quantify its rotation (Figure 2A, see also Supplementary Video S1), examined its specific interactions with the rest of the protein (Supplementary Figure S1), and analyzed the mutual distances with key interacting residues (Figures 2B, C, see also Supplementary Video S2). This analysis revealed that, in the presence of strychnine, Y241^6.48 is oriented toward the centre of the 7TM bundle, forming a stable hydrogen bond with N92^3.36 on helix TM3 in most simulation frames. In contrast, in the Apo state, Y241^6.48 interacts with helix TM7 via a hydrogen bond with T274^7.48. Notably, in the Trans state (after strychnine removal), residue Y241^6.48 was unable to form stable interactions with either N92^3.36 or T274^7.48. The interaction of Y241^6.48 and N92^3.36 was also reported in previous literature regarding TAS2R46 bound to strychnine (Sandal et al., 2015). Additionally, Y241^6.48 and N92^3.36 demonstrate high conservation among TAS2Rs (64% and 84%, respectively) and the mutations of such residues lead to reduced sensitivity and response levels, also in other TAS2Rs, confirming their importance in the activation process (Pronin et al., 2004; Sandal et al., 2015). Recent literature has also pointed out a similar interaction between residues Y^6.48 and N^3.36 through a hydrogen bond in the predicted active state of TAS2R14 (Tokmakova et al., 2023). Interestingly, in the same study, one of the agonists, Flufenamic acid (FFA), was shown to interact with residues at positions 6.48 and 3.36 (or 3.37), thereby preventing these residues from forming a hydrogen bond with each other. Despite this differing interaction profile of FFA and strychnine, the stabilizing mechanism by which these agonists maintain proximity between TM3 and TM6 appears to be similar. Indeed, in our study, strychnine was essential for preserving the hydrogen bond between residues Y241^6.48 and N92^3.36 (Supplementary Figure S7). In contrast, in class A GPCRs, the side chain of the highly conserved residue W^6.48 alternates between Gauche+ and trans conformations for the active and inactive states, respectively (Tokmakova et al., 2023). This highlights a distinct difference in behaviour between class A receptors and TAS2Rs. Specifically, residue W^6.48 in class A GPCRs is able to create a bridge between TM6 and TM3 through an interaction with residue in position 3.36 in the inactive state (Trzaskowski et al., 2012; Venkatakrishnan et al., 2013).

Concerning the dynamic network analysis involving residue Y241^6.48, we pointed out the formation in the Apo state of an edge connecting TM3 and TM7 due to the orientation of residue Y241^6.48 towards the TM7. At the same time, the influence of TM6 in the network is reduced. On the other hand, the intermediate localization of the Y241^6.48 side chain for the Trans system did not allow the connection of TM3 either with TM6 or TM7, again placing the behaviour of the Trans structure between Holo and Apo states. These network analysis results mirror the positioning of residue Y241^6.48 within the TM bundle across different states, highlighting its ability to form stable or unstable interactions with TM3 or TM7, depending on the system being analysed, as previously discussed in the interaction analysis (see also Supplementary Figure S7 and Figure 2). The pivotal position of Y241^6.48 was also confirmed by its implication in the optimal paths linking strychnine binding residues (W88^3.32 and E265^7.39) to previously identified G-protein binding residues (H224^6.31 and Y106^3.50) (Xu et al., 2022), as represented in Figure 6 and reported in Supplementary Table S1. It is worth mentioning that residues H224^6.31 and Y106^3.50 are highly conserved among TAS2Rs (88% and 92%, respectively), underscoring their importance for the majority of the bitter taste receptors. In all systems (Apo, Trans and Holo) the communication between W88^3.32 and H224^6.31 is mediated by TM3, whereas the connection between E265^7.39 and Y106^3.50 is facilitated by TM6 only in the Holo state. In the presence of strychnine (Holo state), both paths involve the edge between N92^3.36 and Y241^6.48, which creates a direct bridge between TM3 and TM6. On the contrary, in the absence of the ligand (Apo and Trans states) the path from W88^3.32 to H224^6.31 involves the creation of an edge connecting TM3 and TM5, before reaching TM6. On the other hand, the connection between E265^7.39 and Y106^3.50 is mediated by a direct connection linking TM3 and TM7 (Y271^7.45-N96^3.40) in the Apo system, whereas this path involves the formation of the network between TM7-TM2-TM3 in the Trans state. In summary, it is noteworthy that the dynamic network analysis revealed that the mechanical information is transmitted from the orthosteric binding site to the intracellular (IC) domain via a pathway involving TM3-TM6 in the presence of strychnine, directly implicating residue Y241^6.48 in the network edge. As also highlighted by the interaction and contact analysis (see also Figure 2; Supplementary Figure S7), the presence of strychnine is fundamental to preserve the hydrogen bond between residues Y241^6.48 and N92^3.36. Additionally, the removal of the ligand from the binding site caused the Trans structure to form pathways resembling those of the Apo state, connecting residues W88^3.32 to H224^6.31.

In conclusion, this study suggests that the “toggle-switch” mechanism, represented by the side chain rotation of Y241^6.48, enables the formation of an allosteric network that links the extracellular (EC) and intracellular (IC) domains, without triggering important conformational changes in the IC region of TAS2R46. More in detail, the rotation of Y241 in TM6 towards the centre of the 7TM bundle allows for the formation of an interaction with N92^3.36 in TM3, which forms a bridge through which the mechanical information can be transferred between the two helices whose IC regions are involved in G-protein binding. This is also associated with different dynamical behaviour of the receptor, which, in the presence of strychnine, is characterized by higher correlations between the IC and EC regions. In this way, the allosteric network generated by ligand binding is directly transferred to the G-protein binding sites.

It is essential to emphasize that the present study includes some limitations to be addressed in future works. First, the simulated structures did not include the G-protein within the models and therefore cannot directly quantify the effects of the presence or absence of the bitter agonist on the G-protein. Future studies that incorporate the G-protein into the simulated system could extend the analysis of generalized correlations and dynamic networks directly onto the G-protein, providing deeper insights into the critical communication pathways linking TAS2R46 and the associated G-protein both in the presence and absence of a bitter agonist. Moreover, the TAS2R46 structures were embedded into a homogeneous POPC bilayer, despite it is known that the membrane composition can modulate GPCRs function, stability, and signalling (Gimpl, 2016; Sengupta et al., 2018). Therefore, further studies will be fundamental to expand the results of the present work towards a complete characterization of the TAS2R46 structure and dynamics. However, this study marks a significant step forward in the characterization of TAS2Rs, shedding light on specific structural and conformational features of TAS2R46 by means of molecular dynamics simulations and network-based analysis of experimental structures. Furthermore, the workflow employed in the present study could be extended in future works aimed at understanding the functional effects of specific structural variations, such as mutations. Notably, specific mutations in TAS2Rs have been shown to either enhance or reduce receptor activation, resulting in altered bitterness perception and subsequently influencing individual dietary and beverage preferences (Roudnitzky et al., 2011; Born et al., 2013; Yamaki et al., 2023).

In conclusion, through a computational approach involving molecular dynamics simulations and network-based analysis of experimentally derived structures, we provided key insights into conformational properties, global intra-structural correlations and allosteric networks associated with TAS2R46 bound or unbound with a bitter agonist. The results highlighted that TM3 and TM6 are the main helices involved in the allosteric network when the receptor is bound to strychnine, while TM6 reduces its influence when strychnine is absent. Moreover, molecular simulations confirmed previous evidence regarding the importance of Y241^6.48 side-chain localization and provided additional insights into the effect of such localization on the allosteric network of the receptor. Lastly, we emphasised that the presence of the bitter agonist increased the overall intra-structural correlations, and the removal of the ligand resulted in a loss of these correlations. The results of the present work also highlighted the similarities and differences between class A GPCRs and TAS2Rs, supporting the hypothesis that bitter taste receptors should be classified into a distinct family, characterised by unique structural and dynamic features. Hence, the proposed methodology could potentially be expanded to characterise other TAS2Rs or other bitter agonists, aiming to identify both common and specific features, which may vary given the wide variability in this fascinating class of receptors. Therefore, this work serves as a crucial starting point for a comprehensive understanding of the functional mechanisms of TAS2Rs and their agonists, laying the groundwork for characterizing their activation processes and exploring how to specifically and effectively target these receptors with tailored agonists.

Data availability statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://github.com/lorenzopallante/TAS2R46.

Author contributions

MC: Conceptualization, Data curation, Formal Analysis, Methodology, Software, Validation, Visualization, Writing–original draft, Writing–review and editing. RF: Conceptualization, Data curation, Formal Analysis, Investigation, Validation, Visualization, Writing–original draft, Writing–review and editing. EZ: Supervision, Writing–original draft, Writing–review and editing. MM: Supervision, Visualization, Writing–original draft, Writing–review and editing. GDB: Supervision, Writing–original draft, Writing–review and editing. Marco Agostino MD: Conceptualization, Data curation, Funding acquisition, Project administration, Resources, Supervision, Validation, Visualization, Writing–original draft, Writing–review and editing. LP: Conceptualization, Data curation, Formal Analysis, Methodology, Software, Supervision, Validation, Visualization, Writing–original draft, Writing–review and editing.

Funding

The author(s) declare that financial support was received for the research, authorship, and/or publication of this article. The present work has been developed as part of the VIRTUOUS project, funded by the European Union’s Horizon 2020 research and innovation program under the Marie Sklodowska-Curie-RISE Grant Agreement No. 872181 (https://www.virtuoush2020.com/). This work has been conducted as part of the GALATEA project, which is funded by the European Union’s Horizon Europe research and innovation program under the Marie Skłodowska-Curie Actions (HORIZON-MSCA-2023-SE-01-01, Grant Agreement No. 101183057).

Conflict of interest

Author GDB was employed by 7HC s.r.l.

The remaining authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

The author(s) declared that they were an editorial board member of Frontiers, at the time of submission. This had no impact on the peer review process and the final decision.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmolb.2024.1473675/full#supplementary-material

References

Ballesteros, J. A., and Weinstein, H. (1995). “[19] Integrated methods for the construction of three-dimensional models and computational probing of structure-function relations in G protein-coupled receptors,” in Methods in Neurosciences (Elsevier), 366–428. doi:10.1016/S1043-9471(05)80049-7

CrossRef Full Text | Google Scholar

Bauer, P., Hess, B., and Lindahl, E. (2023). GROMACS 2022.5 manual. doi:10.5281/ZENODO.7586765

CrossRef Full Text | Google Scholar

Behrens, M., and Lang, T. (2022). Extra-oral taste receptors—function, disease, and perspectives. Front. Nutr. 9, 881177. doi:10.3389/fnut.2022.881177

PubMed Abstract | CrossRef Full Text | Google Scholar

Behrens, M., and Meyerhof, W. (2009). Mammalian bitter taste perception. Results Problems Cell Differ. 47, 203–220. doi:10.1007/400_2008_5

CrossRef Full Text | Google Scholar

Behrens, M., and Meyerhof, W. (2013). Bitter taste receptor research comes of age: from characterization to modulation of TAS2Rs. Seminars Cell and Dev. Biol. 24, 215–221. doi:10.1016/j.semcdb.2012.08.006

CrossRef Full Text | Google Scholar

Behrens, M., and Ziegler, F. (2020). Structure-function analyses of human bitter taste receptors—where do we stand? Molecules 25, 4423. doi:10.3390/molecules25194423

PubMed Abstract | CrossRef Full Text | Google Scholar

Berendsen, H. J. C., Postma, J. P. M., Van Gunsteren, W. F., DiNola, A., and Haak, J. R. (1984). Molecular dynamics with coupling to an external bath. J. Chem. Phys. 81, 3684–3690. doi:10.1063/1.448118

CrossRef Full Text | Google Scholar

Bertalan, É., Lešnik, S., Bren, U., and Bondar, A.-N. (2020). Protein-water hydrogen-bond networks of G protein-coupled receptors: graph-based analyses of static structures and molecular dynamics. J. Struct. Biol. 212, 107634. doi:10.1016/j.jsb.2020.107634

PubMed Abstract | CrossRef Full Text | Google Scholar

Bhattacharya, S., and Vaidehi, N. (2014). Differences in allosteric communication pipelines in the inactive and active states of a GPCR. Biophysical J. 107, 422–434. doi:10.1016/j.bpj.2014.06.015

CrossRef Full Text | Google Scholar

Bock, A., and Bermudez, M. (2021). Allosteric coupling and biased agonism in G protein-coupled receptors. FEBS J. 288, 2513–2528. doi:10.1111/febs.15783

PubMed Abstract | CrossRef Full Text | Google Scholar

Bondar, A.-N. (2022). Graphs of hydrogen-bond networks to dissect protein conformational dynamics. J. Phys. Chem. B 126, 3973–3984. doi:10.1021/acs.jpcb.2c00200

PubMed Abstract | CrossRef Full Text | Google Scholar

Born, S., Levit, A., Niv, M. Y., Meyerhof, W., and Behrens, M. (2013). The human bitter taste receptor TAS2R10 is tailored to accommodate numerous diverse ligands. J. Neurosci. 33, 201–213. doi:10.1523/JNEUROSCI.3248-12.2013

PubMed Abstract | CrossRef Full Text | Google Scholar

Brockhoff, A., Behrens, M., Massarotti, A., Appendino, G., and Meyerhof, W. (2007). Broad tuning of the human bitter taste receptor hTAS2R46 to various sesquiterpene lactones, clerodane and labdane diterpenoids, strychnine, and denatonium. J. Agric. Food Chem. 55, 6236–6243. doi:10.1021/jf070503p

PubMed Abstract | CrossRef Full Text | Google Scholar

Brockhoff, A., Behrens, M., Niv, M. Y., and Meyerhof, W. (2010). Structural requirements of bitter taste receptor activation. Proc. Natl. Acad. Sci. U. S. A. 107, 11110–11115. doi:10.1073/pnas.0913862107

PubMed Abstract | CrossRef Full Text | Google Scholar

Cannariato, M., Zizzi, E. A., Pallante, L., Miceli, M., and Deriu, M. A. (2023). Mechanical communication within the microtubule through network-based analysis of tubulin dynamics. Biomech. Model Mechanobiol. 23, 569–579. doi:10.1007/s10237-023-01792-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Chandrashekar, J., Mueller, K. L., Hoon, M. A., Adler, E., Feng, L., Guo, W., et al. (2000). T2Rs function as bitter taste receptors. Cell 100, 703–711. doi:10.1016/S0092-8674(00)80706-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Dagan-Wiener, A., Di Pizio, A., Nissim, I., Bahia, M. S., Dubovski, N., Margulis, E., et al. (2019). BitterDB: taste ligands and receptors database in 2019. Nucleic Acids Res. 47, D1179–D1185. doi:10.1093/nar/gky974

PubMed Abstract | CrossRef Full Text | Google Scholar

Dalton, J. A., Lans, I., and Giraldo, J. (2015). Quantifying conformational changes in GPCRs: glimpse of a common functional mechanism. BMC Bioinforma. 16, 124. doi:10.1186/s12859-015-0567-3

CrossRef Full Text | Google Scholar

Dickson, C. J., Walker, R. C., and Gould, I. R. (2022). Lipid21: complex lipid membrane simulations with AMBER. J. Chem. Theory Comput. 18, 1726–1736. doi:10.1021/acs.jctc.1c01217

PubMed Abstract | CrossRef Full Text | Google Scholar

Di Pizio, A., Kruetzfeldt, L. M., Cheled-Shoval, S., Meyerhof, W., Behrens, M., and Niv, M. Y. (2017). Ligand binding modes from low resolution GPCR models and mutagenesis: chicken bitter taste receptor as a test-case. Sci. Rep. 7, 1–11. doi:10.1038/s41598-017-08344-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Di Pizio, A., Levit, A., Slutzki, M., Behrens, M., Karaman, R., and Niv, M. Y. (2016). “Comparing Class A GPCRs to bitter taste receptors,” in Methods in cell biology (Elsevier), 401–427. doi:10.1016/bs.mcb.2015.10.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Di Pizio, A., Niv, M., Waterloo, L., Brox, R., Löber, S., Weikert, D., et al. (2020). Rational design of agonists for bitter taste receptor TAS2R14: from modeling to bench and back. Cell. Mol. Life Sci. 77, 531–542. doi:10.1007/s00018-019-03194-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Di Pizio, A., and Niv, M. Y. (2015). Promiscuity and selectivity of bitter molecules and their receptors. Bioorg. and Med. Chem. 23, 4082–4091. doi:10.1016/j.bmc.2015.04.025

CrossRef Full Text | Google Scholar

Dotson, C. D., Zhang, L., Xu, H., Shin, Y.-K., Vigues, S., Ott, S. H., et al. (2008). Bitter taste receptors influence glucose homeostasis. PLoS ONE 3, e3974. doi:10.1371/journal.pone.0003974

PubMed Abstract | CrossRef Full Text | Google Scholar

Fanelli, F., Felline, A., Raimondi, F., and Seeber, M. (2016). Structure network analysis to gain insights into GPCR function. Biochem. Soc. Trans. 44, 613–618. doi:10.1042/BST20150283

PubMed Abstract | CrossRef Full Text | Google Scholar

Feng, S., Park, S., Choi, Y. K., and Im, W. (2023). CHARMM-GUI membrane builder: past, current, and future developments and applications. J. Chem. Theory Comput. 19, 2161–2185. doi:10.1021/acs.jctc.2c01246

PubMed Abstract | CrossRef Full Text | Google Scholar

Fierro, F., Peri, L., Hübner, H., Tabor-Schkade, A., Waterloo, L., Löber, S., et al. (2022). Inhibiting a promiscuous GPCR: iterative discovery of bitter taste receptor ligands. Bioinformatics. doi:10.1101/2022.11.24.517821

CrossRef Full Text | Google Scholar

Gimpl, G. (2016). Interaction of G protein coupled receptors and cholesterol. Chem. Phys. Lipids 199, 61–73. doi:10.1016/j.chemphyslip.2016.04.006

PubMed Abstract | CrossRef Full Text | Google Scholar

Heinig, M., and Frishman, D. (2004). STRIDE: a web server for secondary structure assignment from known atomic coordinates of proteins. Nucleic Acids Res. 32, W500–W502. doi:10.1093/nar/gkh429

PubMed Abstract | CrossRef Full Text | Google Scholar

Hénin, J., Maigret, B., Tarek, M., Escrieut, C., Fourmy, D., and Chipot, C. (2006). Probing a model of a GPCR/ligand complex in an explicit membrane environment: the human cholecystokinin-1 receptor. Biophysical J. 90, 1232–1240. doi:10.1529/biophysj.105.070599

CrossRef Full Text | Google Scholar

Hess, B., Bekker, H., Berendsen, H. J. C., and Fraaije, J. G. E. M. (1997). LINCS: a linear constraint solver for molecular simulations. J. Comput. Chem. 18, 1463–1472. doi:10.1002/(SICI)1096-987X(199709)18:12<1463::AID-JCC4>3.0.CO;2-H

CrossRef Full Text | Google Scholar

Humphrey, W., Dalke, A., and Schulten, K. (1996). VMD: visual molecular dynamics. J. Mol. Graph. 14, 33–28. doi:10.1016/0263-7855(96)00018-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Hunter, J. D. (2007). Matplotlib: a 2D graphics environment. Comput. Sci. Eng. 9, 90–95. doi:10.1109/MCSE.2007.55

CrossRef Full Text | Google Scholar

Ibrahim, P., Wifling, D., and Clark, T. (2019). Universal activation index for class A GPCRs. J. Chem. Inf. Model. 59, 3938–3945. doi:10.1021/acs.jcim.9b00604

PubMed Abstract | CrossRef Full Text | Google Scholar

Jensen, A. A., Gharagozloo, P., Birdsall, N. J. M., and Zlotos, D. P. (2006). Pharmacological characterisation of strychnine and brucine analogues at glycine and alpha7 nicotinic acetylcholine receptors. Eur. J. Pharmacol. 539, 27–33. doi:10.1016/j.ejphar.2006.04.010

PubMed Abstract | CrossRef Full Text | Google Scholar

Jumper, J., Evans, R., Pritzel, A., Green, T., Figurnov, M., Ronneberger, O., et al. (2021). Highly accurate protein structure prediction with AlphaFold. Nature 596, 583–589. doi:10.1038/s41586-021-03819-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Kraskov, A., Stögbauer, H., and Grassberger, P. (2004). Estimating mutual information. Phys. Rev. E 69, 066138. doi:10.1103/PhysRevE.69.066138

CrossRef Full Text | Google Scholar

Lange, O. F., and Grubmüller, H. (2005). Generalized correlation for biomolecular dynamics. Proteins 62, 1053–1061. doi:10.1002/prot.20784

CrossRef Full Text | Google Scholar

Laskowski, R. A., MacArthur, M. W., Moss, D. S., and Thornton, J. M. (1993). PROCHECK: a program to check the stereochemical quality of protein structures. J. Appl. Crystallogr. 26, 283–291. doi:10.1107/s0021889892009944

CrossRef Full Text | Google Scholar

Latorraca, N. R., Venkatakrishnan, A. J., and Dror, R. O. (2017). GPCR dynamics: structures in motion. Chem. Rev. 117, 139–155. doi:10.1021/acs.chemrev.6b00177

PubMed Abstract | CrossRef Full Text | Google Scholar

Laurent, B., Chavent, M., Cragnolini, T., Dahl, A. C. E., Pasquali, S., Derreumaux, P., et al. (2015). Epock: rapid analysis of protein pocket dynamics. Bioinformatics 31, 1478–1480. doi:10.1093/bioinformatics/btu822

PubMed Abstract | CrossRef Full Text | Google Scholar

Lee, J., Hitzenberger, M., Rieger, M., Kern, N. R., Zacharias, M., and Im, W. (2020). CHARMM-GUI supports the Amber force fields. J. Chem. Phys. 153, 035103. doi:10.1063/5.0012280

PubMed Abstract | CrossRef Full Text | Google Scholar

Lee, S.-J., Depoortere, I., and Hatt, H. (2019). Therapeutic potential of ectopic olfactory and taste receptors. Nat. Rev. Drug Discov. 18, 116–138. doi:10.1038/s41573-018-0002-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Levit, A., Beuming, T., Krilov, G., Sherman, W., and Niv, M. Y. (2014). Predicting GPCR promiscuity using binding site features. J. Chem. Inf. Model. 54, 184–194. doi:10.1021/ci400552z

PubMed Abstract | CrossRef Full Text | Google Scholar

Liggett, S. B. (2014). Bitter taste receptors in the wrong place: novel airway smooth muscle targets for treating asthma. Trans. Am. Clin. Climatol. Assoc. 125, 64–75.

PubMed Abstract | Google Scholar

Lu, P., Zhang, C.-H., Lifshitz, L. M., and ZhuGe, R. (2017). Extraoral bitter taste receptors in health and disease. J. General Physiology 149, 181–197. doi:10.1085/jgp.201611637

CrossRef Full Text | Google Scholar

Malavolta, M., Pallante, L., Mavkov, B., Stojceski, F., Grasso, G., Korfiati, A., et al. (2022). A survey on computational taste predictors. Eur. Food Res. Technol. 248, 2215–2235. doi:10.1007/s00217-022-04044-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Manrique, P. D., Chakraborty, S., Henderson, R., Edwards, R. J., Mansbach, R., Nguyen, K., et al. (2023). Network analysis uncovers the communication structure of SARS-CoV-2 spike protein identifying sites for immunogen design. iScience 26, 105855. doi:10.1016/j.isci.2022.105855

PubMed Abstract | CrossRef Full Text | Google Scholar

Marchiori, A., Capece, L., Giorgetti, A., Gasparini, P., Behrens, M., Carloni, P., et al. (2013). Coarse-grained/molecular mechanics of the TAS2R38 bitter taste receptor: experimentally-validated detailed structural prediction of agonist binding. PLoS ONE 8, e64675. doi:10.1371/journal.pone.0064675

PubMed Abstract | CrossRef Full Text | Google Scholar

Melo, M. C. R., Bernardi, R. C., De La Fuente-Nunez, C., and Luthey-Schulten, Z. (2020). Generalized correlation-based dynamical network analysis: a new high-performance approach for identifying allosteric communications in molecular dynamics trajectories. J. Chem. Phys. 153, 134104. doi:10.1063/5.0018980

PubMed Abstract | CrossRef Full Text | Google Scholar

Miceli, M., Deriu, M. A., and Grasso, G. (2022). Toward the design and development of peptidomimetic inhibitors of the Ataxin-1 aggregation pathway. Biophysical J. 121, 4679–4688. doi:10.1016/j.bpj.2022.10.021

CrossRef Full Text | Google Scholar

Molecular Operating Environment (MOE) (2022). 02 chemical computing Group ULC, 1010 Sherbooke St. West, suite #910. Montreal, QC, Canada. H3A 2R7, 2023.

Google Scholar

Nivedha, A. K., Tautermann, C. S., Bhattacharya, S., Lee, S., Casarosa, P., Kollak, I., et al. (2018). Identifying functional hotspot residues for biased ligand design in G-protein-coupled receptors. Mol. Pharmacol. 93, 288–296. doi:10.1124/mol.117.110395

PubMed Abstract | CrossRef Full Text | Google Scholar

Nowak, S., Di Pizio, A., Levit, A., Niv, M. Y., Meyerhof, W., and Behrens, M. (2018). Reengineering the ligand sensitivity of the broadly tuned human bitter taste receptor TAS2R14. Biochimica Biophysica Acta - General Subj. 1862, 2162–2173. doi:10.1016/j.bbagen.2018.07.009

CrossRef Full Text | Google Scholar

Nygaard, R., Frimurer, T. M., Holst, B., Rosenkilde, M. M., and Schwartz, T. W. (2009). Ligand binding and micro-switches in 7TM receptor structures. Trends Pharmacol. Sci. 30, 249–259. doi:10.1016/j.tips.2009.02.006

PubMed Abstract | CrossRef Full Text | Google Scholar

Pallante, L., Cannariato, M., Androutsos, L., Zizzi, E. A., Bompotas, A., Hada, X., et al. (2024). VirtuousPocketome: a computational tool for screening protein–ligand complexes to identify similar binding sites. Sci. Rep. 14, 6296. doi:10.1038/s41598-024-56893-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Pallante, L., Cannariato, M., Vezzulli, F., Malavolta, M., Lambri, M., and Deriu, M. A. (2023). Machine learning aided molecular modelling of taste to identify food fingerprints. Chem. Eng. Trans. 102, 283–288. doi:10.3303/CET23102048

CrossRef Full Text | Google Scholar

Pallante, L., Malavolta, M., Grasso, G., Korfiati, A., Mavroudi, S., Mavkov, B., et al. (2021). On the human taste perception: molecular-level understanding empowered by computational methods. Trends Food Sci. and Technol. 116, 445–459. doi:10.1016/j.tifs.2021.07.013

CrossRef Full Text | Google Scholar

Pan, S., Sharma, P., Shah, S. D., and Deshpande, D. A. (2017). Bitter taste receptor agonists alter mitochondrial function and induce autophagy in airway smooth muscle cells. Am. J. Physiology-Lung Cell. Mol. Physiology 313, L154–L165. doi:10.1152/ajplung.00106.2017

CrossRef Full Text | Google Scholar

Pándy-Szekeres, G., Caroli, J., Mamyrbekov, A., Kermani, A. A., Keserű, G. M., Kooistra, A. J., et al. (2023). GPCRdb in 2023: state-specific structure models using AlphaFold2 and new ligand resources. Nucleic Acids Res. 51, D395–D402. doi:10.1093/nar/gkac1013

PubMed Abstract | CrossRef Full Text | Google Scholar

Pronin, A. N., Tang, H., Connor, J., and Keung, W. (2004). Identification of ligands for two human bitter T2R receptors. Chem. Senses 29, 583–593. doi:10.1093/chemse/bjh064

PubMed Abstract | CrossRef Full Text | Google Scholar

Pydi, S. P., Bhullar, R. P., and Chelikani, P. (2012). Constitutively active mutant gives novel insights into the mechanism of bitter taste receptor activation. J. Neurochem. 122, 537–544. doi:10.1111/j.1471-4159.2012.07808.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Rasmussen, S. G. F., DeVree, B. T., Zou, Y., Kruse, A. C., Chung, K. Y., Kobilka, T. S., et al. (2011). Crystal structure of the β2 adrenergic receptor-Gs protein complex. Nature 477, 549–555. doi:10.1038/nature10361

PubMed Abstract | CrossRef Full Text | Google Scholar

Roper, S. D. (2017). “Taste: mammalian taste bud physiology,” in Reference module in neuroscience and biobehavioral psychology (Elsevier). doi:10.1016/B978-0-12-809324-5.02908-4

CrossRef Full Text | Google Scholar

Roudnitzky, N., Bufe, B., Thalmann, S., Kuhn, C., Gunn, H. C., Xing, C., et al. (2011). Genomic, genetic and functional dissection of bitter taste responses to artificial sweeteners. Hum. Mol. Genet. 20, 3437–3449. doi:10.1093/hmg/ddr252

PubMed Abstract | CrossRef Full Text | Google Scholar

Sakurai, T., Misaka, T., Ishiguro, M., Masuda, K., Sugawara, T., Ito, K., et al. (2010). Characterization of the beta-D-glucopyranoside binding site of the human bitter taste receptor hTAS2R16. J. Biol. Chem. 285, 28373–28378. doi:10.1074/jbc.M110.144444

PubMed Abstract | CrossRef Full Text | Google Scholar

Salentin, S., Schreiber, S., Haupt, V. J., Adasme, M. F., and Schroeder, M. (2015). PLIP: fully automated protein–ligand interaction profiler. Nucleic Acids Res. 43, W443–W447. doi:10.1093/nar/gkv315

PubMed Abstract | CrossRef Full Text | Google Scholar

Sandal, M., Behrens, M., Brockhoff, A., Musiani, F., Giorgetti, A., Carloni, P., et al. (2015). Evidence for a transient additional ligand binding site in the TAS2R46 bitter taste receptor. J. Chem. Theory Comput. 11, 4439–4449. doi:10.1021/acs.jctc.5b00472

PubMed Abstract | CrossRef Full Text | Google Scholar

Scheer, A., and Cotecchia, S. (1997). Constitutively active G protein-coupled receptors: potential mechanisms of receptor activation. J. Recept. Signal Transduct. 17, 57–73. doi:10.3109/10799899709036594

PubMed Abstract | CrossRef Full Text | Google Scholar

Sengupta, D., Prasanna, X., Mohole, M., and Chattopadhyay, A. (2018). Exploring GPCR–lipid interactions by molecular dynamics simulations: excitements, challenges, and the way forward. J. Phys. Chem. B 122, 5727–5737. doi:10.1021/acs.jpcb.8b01657

PubMed Abstract | CrossRef Full Text | Google Scholar

Siemers, M., Lazaratos, M., Karathanou, K., Guerra, F., Brown, L. S., and Bondar, A.-N. (2019). Bridge: a graph-based algorithm to analyze dynamic H-bond networks in membrane proteins. J. Chem. Theory Comput. 15, 6781–6798. doi:10.1021/acs.jctc.9b00697

PubMed Abstract | CrossRef Full Text | Google Scholar

Singh, N., Pydi, S. P., Upadhyaya, J., and Chelikani, P. (2011). Structural basis of activation of bitter taste receptor T2R1 and comparison with class A G-protein-coupled receptors (GPCRs). J. Biol. Chem. 286, 36032–36041. doi:10.1074/jbc.M111.246983

PubMed Abstract | CrossRef Full Text | Google Scholar

Sullivan, H.-J., Tursi, A., Moore, K., Campbell, A., Floyd, C., and Wu, C. (2020). Binding interactions of ergotamine and dihydroergotamine to 5-hydroxytryptamine receptor 1B (5-HT _1b) using molecular dynamics simulations and dynamic network analysis. J. Chem. Inf. Model. 60, 1749–1765. doi:10.1021/acs.jcim.9b01082

PubMed Abstract | CrossRef Full Text | Google Scholar

Tian, C., Kasavajhala, K., Belfon, K. A. A., Raguette, L., Huang, H., Migues, A. N., et al. (2020). ff19SB: amino-acid-specific protein backbone parameters trained against quantum mechanics energy surfaces in solution. J. Chem. Theory Comput. 16, 528–552. doi:10.1021/acs.jctc.9b00591

PubMed Abstract | CrossRef Full Text | Google Scholar

Tokmakova, A., Kim, D., Guthrie, B., Kim, S.-K., Goddard, W. A., and Liggett, S. B. (2023). Predicted structure and cell signaling of TAS2R14 reveal receptor hyper-flexibility for detecting diverse bitter tastes. iScience 26, 106422. doi:10.1016/j.isci.2023.106422

PubMed Abstract | CrossRef Full Text | Google Scholar

Topin, J., Bouysset, C., Pacalon, J., Kim, Y., Rhyu, M., Fiorucci, S., et al. (2021). Functional molecular switches of mammalian G protein-coupled bitter-taste receptors. Cell. Mol. Life Sci. 6, 7605–7615. doi:10.1007/s00018-021-03968-7

CrossRef Full Text | Google Scholar

Trzaskowski, B., Latek, D., Yuan, S., Ghoshdastider, U., Debinski, A., and Filipek, S. (2012). Action of molecular switches in GPCRs - theoretical and experimental studies. CMC 19, 1090–1109. doi:10.2174/092986712799320556

CrossRef Full Text | Google Scholar

Venkatakrishnan, A. J., Deupi, X., Lebon, G., Heydenreich, F. M., Flock, T., Miljus, T., et al. (2016). Diverse activation pathways in class A GPCRs converge near the G-protein-coupling region. Nature 536, 484–487. doi:10.1038/nature19107

PubMed Abstract | CrossRef Full Text | Google Scholar

Venkatakrishnan, A. J., Deupi, X., Lebon, G., Tate, C. G., Schertler, G. F., and Babu, M. M. (2013). Molecular signatures of G-protein-coupled receptors. Nature 494, 185–194. doi:10.1038/nature11896

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, J., Wolf, R. M., Caldwell, J. W., Kollman, P. A., and Case, D. A. (2004). Development and testing of a general amber force field. J. Comput. Chem. 25, 1157–1174. doi:10.1002/jcc.20035

PubMed Abstract | CrossRef Full Text | Google Scholar

Waskom, M. (2021). seaborn: statistical data visualization. JOSS 6, 3021. doi:10.21105/joss.03021

CrossRef Full Text | Google Scholar

Waterhouse, A., Bertoni, M., Bienert, S., Studer, G., Tauriello, G., Gumienny, R., et al. (2018). SWISS-MODEL: homology modelling of protein structures and complexes. Nucleic Acids Res. 46, W296–W303. doi:10.1093/nar/gky427

PubMed Abstract | CrossRef Full Text | Google Scholar

Webb, B., and Sali, A. (2016). Comparative protein structure modeling using MODELLER. CP Bioinforma. 54, 5.6.1–5.6.37. doi:10.1002/cpbi.3

CrossRef Full Text | Google Scholar

Workman, A. D., Palmer, J. N., Adappa, N. D., and Cohen, N. A. (2015). The role of bitter and sweet taste receptors in upper airway immunity. Curr. Allergy Asthma Rep. 15, 72. doi:10.1007/s11882-015-0571-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Xu, J., and Zhang, Y. (2010). How significant is a protein structure similarity with TM-score = 0.5? Bioinformatics 26, 889–895. doi:10.1093/bioinformatics/btq066

PubMed Abstract | CrossRef Full Text | Google Scholar

Xu, W., Wu, L., Liu, S., Liu, X., Cao, X., Zhou, C., et al. (2022). Structural basis for strychnine activation of human bitter taste receptor TAS2R46. Science 377, 1298–1304. doi:10.1126/science.abo1633

PubMed Abstract | CrossRef Full Text | Google Scholar

Xue, A. Y., Di Pizio, A., Levit, A., Yarnitzky, T., Penn, O., Pupko, T., et al. (2018). Independent evolution of strychnine recognition by bitter taste receptor subtypes. Front. Mol. Biosci. 5, 9. doi:10.3389/fmolb.2018.00009

PubMed Abstract | CrossRef Full Text | Google Scholar

Yamaki, M., Saito, H., Mimori, T., Suzuki, Y., Nagasaki, M., Suzuki, K., et al. (2023). Analysis of genetic polymorphism of bitter taste receptor TAS2R38 and TAS2R46, and its relationship with eating and drinking habits in Japanese ToMMo subjects. J. Nutr. Sci. Vitaminol. 69, 347–356. doi:10.3177/jnsv.69.347

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, Y., and Skolnick, J. (2004). Scoring function for automated assessment of protein structure template quality. Proteins 57, 702–710. doi:10.1002/prot.20264

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhou, Q., Yang, D., Wu, M., Guo, Y., Guo, W., Zhong, L., et al. (2019). Common activation mechanism of class A GPCRs. eLife 8, e50279. doi:10.7554/eLife.50279

PubMed Abstract | CrossRef Full Text | Google Scholar

Zou, E., Ewalt, J., and Ng, H. L. (2019). Recent insights from molecular dynamics simulations for G protein-coupled receptor drug discovery. IJMS 20, 4237. doi:10.3390/ijms20174237

PubMed Abstract | CrossRef Full Text | Google Scholar