Multi-scale simulations of membrane proteins: The case of bitter taste receptors Journal of Science: Advanced Materials and Devices

Human bitter taste receptors (hTAS2Rs) are the second largest group of chemosensory G-protein coupled receptors (25 members). hTAS2Rs are expressed in many tissues (e.g. tongue, gastrointestinal tract, respiratory system, brain, etc.), performing a variety of functions, from bitter taste perception to hormone secretion and bronchodilation. Due to the lack of experimental structural information, computations are currently the methods of choice to get insights into ligand e receptor interactions. Here we review our efforts at predicting the binding pose of agonists to hTAS2Rs, using state-of-the-art bioinformatics approaches followed by hybrid Molecular Mechanics/Coarse-Grained (MM/CG) simulations. The latter method, developed by us, describes atomistically only the agonist binding region, including hydration, and it may be particularly suited to be used when bioinformatics predictions generate very low- resolution models, such as the case of hTAS2Rs. Our structural predictions of the hTAS2R38 and hTAS2R46 receptors in complex with their agonists turn out to be fully consistent with experimental mutagenesis data. In addition, they suggest a two-binding site architecture in hTAS2R46, consisting of the usual orthosteric site together with a “ vestibular ” site toward the extracellular space, as observed in other GPCRs. The presence of the vestibular site may help to discriminate among the wide spectrum of bitter ligands. © 2017 The Authors. Publishing services by Elsevier B.V. on behalf of Vietnam National University, Hanoi. This


Introduction
The 25 human bitter taste receptors (hTAS2Rs) [1,2] constitute the second largest group of chemosensory G-protein coupled receptors (GPCRs), in turn the largest membrane protein superfamily, with about 850 members in humans. hTAS2Rs are found in many different tissues of the human body [3e5]. These include the plasma membrane of the type II taste receptor cells (from which their name, TAS2Rs, comes from), located in the taste buds of the tongue [1,6e8], the respiratory system [9e11], the gastrointestinal tract [12,13] the endocrine system [13] and the brain [14]. Hence, hTAS2Rs play different roles, ranging from perception of bitter taste, to detection of toxins [15], to bronchodilation [16], and to hormone secretion [17]. hTAS2Rs can recognize hundreds of structurally diverse agonists using a combinatorial coding scheme [18,19]. One hTAS2R is able to recognize more than one agonist [20,21], and one agonist can be recognized by more than one hTAS2R [20]. Understanding the details of hTAS2Rseagonists interactions may provide important hints on the effect of genetic variability on bitter taste perception, and new opportunities for designing more subtype-specific ligands [22,23] and novel therapies against diseases related to hTAS2Rs' dysfunction, e.g. asthma or chronic rhinosinusitis [4,5].
hTAS2Rs, as all GPCRs, are made up of seven transmembrane helices arranged in a helix-bundle shape and connected by three extracellular loops (ECLs) and three intracellular loops (ICLs) [24]. Agonist binding to the receptor's binding site (called orthosteric site) facilitates conformational changes towards an "active state". The latter allows the activation of downstream effectors [7,25]. The location of the orthosteric site in hTAS2Rs is similar to other receptors of the largest GPCR class, class A [26e31]. However, because of the low sequence similarity between hTAS2Rs and the other GPCRs, it has not been clearly established yet if hTAS2Rs belong to class A [32e34] or class F [27,35,36] GPCRs. These proteins could even constitute a different family [27].
At present no experimental structural information is available for hTAS2Rs. Therefore, any attempt at understanding the hTAS2Ragonist complexes has to rely on computational approaches. Bioinformatics techniques, such as homology modeling [37], along with molecular docking [38e40], could in principle provide insights into agonist/antagonist binding. Unfortunately, however, the sequence identity between bitter taste receptors and the possible templates is extremely low (~10e17% with any of the 42 unique Xray structures as of February 2017 (http://blanco.biomol.uci.edu/ mpstruc/)). As a consequence, the construction of reliable alignments between the target sequence and the available structural templates is challenging [41e43]. Moreover, even with a good sequence alignment, the orientation of the side chains in the orthosteric binding site, which is key for proteineligand interactions, is not accurately modeled [44,45]. This hinders the correct prediction of docking poses. In addition, current bioinformatics and docking algorithms face at times limitations (such as the lack of protein flexibility and hydration [46,47]), which may further limit the power of the predictions, especially in light of the fact that factors such as conformational dynamics [48] and water molecules [49,50] play a crucial role for ligand binding and receptor activation. A way to overcome these difficulties is to combine these static computational approaches with molecular simulation techniques, such as molecular dynamics (MD) and enhanced sampling [51e53]. These methods may explore efficiently the conformational space, including hydration and ligandeprotein interactions. All-atom MD has been successfully used in high quality homology models (i.e. based on a template with sequence identity above 60%) [48,54,55]; however, it may provide far less satisfying results when the protein structure is a homology model based on a low sequence identity template (as it is the case for hTAS2Rs). Here, the side chains' rotamers are poorly predicted and often their relaxation requires longer time scales that cannot be reached with atomistic MD. Coarse-grained (CG)-based MD can be used to sample longer timescales [56e59], yet it cannot describe in detail the molecular recognition events between protein and ligand. A way to overcome these limitations is represented by the combination of the two aforementioned techniques [60e67]. In this context, our group has developed a hybrid "Molecular Mechanics/Coarse-Grained" (MM/ CG) method for refinement of GPCRs homology models [63,68,69]. Here, the system is modeled at two different resolutions. While ligand, binding site residues and surrounding water molecules are treated using an atomistic force field, the rest of the protein is described at a CG level. A coupling scheme is then used to connect the two regions at the boundary. This MM/CG method maintains the atomistic resolution needed to describe correctly the proteineligand interactions at the binding site, while allowing a larger conformational sampling and a reduced computational cost compared to an all-atom simulation. The presence of the membrane is mimicked by introducing five repulsive walls. Two planar walls coincide with the height of the head groups of the membrane lipids, two hemispheric walls set a limit on the extracellular and intracellular ends of the protein and the last wall follows the initial shape of the interface between protein and membrane [70e72].
The accuracy of the MM/CG method in reproducing binding poses and protein fluctuations was established in our early work [68]. Here, we will present more recent predictions for widely studied hTAS2Rs, which were successfully validated against extensive mutagenesis data [73e75]. Specifically, we investigate hTAS2R46 [76], a promiscuous bitter taste receptor [20,73,77] that can detect bitter molecules belonging to several different chemical classes, and hTAS2R38 [74,75], a receptor that recognizes agonists containing an isothiocyanate or thiourea group [20,77,78]. Given their different receptive range, these two receptors constitute excellent contrasting test cases to assess the applicability of the MM/CG methodology to study ligand binding in human bitter taste receptors.

Materials and methods
Our web-server GOMoDO [79] performs automatically both the homology modeling and molecular docking steps, by combining state-of-the-art bioinformatics tools for GPCRs. In particular, GOMoDO uses the profileeprofile HMM algorithm (for database search and target-template alignment) and the MODELLER program [80] (for protein homology model construction), followed by information-driven flexible docking of ligands through the HADDOCK program [81]. This protocol was used to produce the initial model of hTAS2R46 in complex with one of its agonists, strychnine, as well as the models of hTAS2R38 in complex with its two agonists, namely propylthiouracil (PROP) and phenylthiocarbamide (PTC). Specifically, the MODELLER algorithm [80] was used to generate 200 models of hTAS2R46 and hTAS2R38, applying a single-template or multiple-template approach, respectively [74e76]. Then, a clustering analysis was performed to identify "representative" receptor models, using as criteria both the MOD-ELLER quality scores and available experimental site-directed mutagenesis data. In the case of hTAS2R46, one single model was taken as representative, whereas for hTAS2R38 two models were selected, which mainly differ in the conformation of the ECL2. The agonists, strychnine for hTAS2R46 and PROP and PTC for hTAS2R38, were docked into the modeled receptor structures using HADDOCK [81]. Information about the putative binding residues was used to drive the docking. For hTAS2R46 the putative binding residues were predicted using FPOCKET [82], whereas for hTAS2R38, they were selected based on previous bioinformatics and site-directed mutagenesis studies [74]. In the docking protocol, first 1000 structures were generated in the rigid body step and, then, the top scoring 200 complexes were further optimized using a flexible simulated annealing step, followed by a final refinement step in explicit water. Next, a clustering analysis was performed to identify the best initial model, that is, the structure of the most populated cluster with the lowest binding energy. The best docking models then underwent MM/CG simulations [63,68,69]. In these multiscale approach, ligand, binding site residues and surrounding water molecules were treated using the GROMOS96 atomistic force field [83], whereas the rest of the protein was described at a CG level, including only the Ca atoms of the amino acids and using a Go-like model [84]. For hTAS2R46, the model of the receptorestrychnine complex was used to set up three replicas, with different initial velocities; a 1 ms-long MM/CG simulation was run for each [76]. For hTAS2R38, the models for each receptor-agonist complex (PROP and PTC) were submitted to MM/CG simulations; for each complex, two replicas, differing only in the initial velocities, were run for 0.6 ms [75].

Results and discussion
3.1. hTAS2R46 in complex with strychnine hTAS2R46 is a receptor involved not only in bitter taste perception, but also in ciliary beat frequency and clearance of microorganisms in the airways [85] and in blood pressure control in vascular smooth muscle cells [86]. hTAS2R46 is a promiscuous receptor [20,73]: it can detect bitter molecules belonging to several different chemical classes. How hTAS2R46 can discriminate this wide range of agonists from other bitter molecules is still an open question. Prof. Meyerhof and coworkers suggested the existence of an "access control" within the extracellular opening of the receptor [73] that may act as a selectivity filter. In an effort at providing a molecular basis of such "access control", we carried out bioinformatics and MM/CG calculations of hTAS2R46 in complex with its agonist strychnine [76].
Interestingly, our simulations identify two different binding poses. In the first pose (Fig. 1a), the ligand is localized in a region that coincides with the orthosteric site identified in the X-ray structures of class A GPCRs in complex with their corresponding agonists. Moreover, like in hTAS2R38 (see below), our MM/CG simulations predict several binding pocket residues, which are subsequently validated through experiments [76]. In particular, Tyr241 and Asn92, which are also highly conserved in the hTAS2R family, are identified in the orthosteric cavity. Tyr241 forms a pstacking interaction with the aromatic ring of the strychnine, as well as a H-bond with Asn92. Consistently, the mutations Asn92Gln and Tyr241Phe in hTAS2R46 reduce the receptor activation levels or abolish the signal completely, respectively, whereas the Tyr241Trp lowers the EC 50 value. Thus, according to these findings, the interaction between Tyr241 and Asn92 could play a role in receptor activation more than in ligand selectivity. In this regard, the latter residue has been shown to be crucial for receptor activation also in hTAS2R43 [87].
In the second pose (Fig. 1b), the ligand is positioned in the extracellular region, in a site that resembles the allosteric binding cavity in class A GPCRs [26,28,88e95], which we called "vestibular" site. Our simulations identify Leu71 and Asn176 as part of the vestibular site and provide a molecular level explanation of previous mutagenesis experimental data [73]. Therefore, the decreased receptor activation for the Leu71Phe mutant is most likely due to a reduction of the volume of the vestibular cavity, whereas, for the Asn176Ala mutant, it is probably caused by the disruption of a Hbond network involving Asn176 and ECL2, a loop known to be involved in ligand binding and receptor activation in GPCRs [89,96,97].
Importantly, some residues known experimentally as functionally important, i.e. Leu71 and Asn176 [73], interact with strychnine only in the vestibular cavity. Hence, the experimental mutagenesis data cannot be rationalized by taking into account only the canonical orthosteric binding site, and, instead, the two topographically distinct ligand binding cavities need to be considered. Therefore, hTAS2R46 features two binding sites (orthosteric and vestibular), and both cavities may contribute to the selectivity of the receptor. In this regard, hTAS2R46 has been found to recognize at least 28 different agonists [77], belonging to diverse chemical classes. Given this promiscuity, it is unlikely the orthosteric binding site alone could discriminate this wide variety of compounds. We hypothesize that the presence of a second, vestibular site can provide additional proteineligand contacts that will help to filter the appropriate agonists out of the pool of more than 100 bitter compounds known [20].
In order to assess whether this two-step mechanism could also apply to other bitter taste receptors, we performed a bioinformatics analysis of the conservation across the hTAS2R family of the residues identified for hTAS2R46 as functionally important (Fig. 2). We found that more than 50% of these residues were conserved in at least two hTAS2Rs. Interestingly, while four of the conserved residues (positions 2.65, 3.26, 3.29 and 5.39, following the generic GPCR numbering [98]) were found to be localized only in the putative vestibular binding site (in red in Fig. 2), five other residues (3.33, 3.36, 3.37, 3.40 and 7.42) were placed only in the orthosteric binding site (in green in Fig. 2). These analyses thus suggest that the two-site architecture may also be present in other human bitter taste receptors, besides hTAS2R46. This could be related to the ability of most hTAS2Rs to detect more than one agonist (see Table 1). Two sites can offer more ligand recognition points than a single one, thus helping to select the appropriate agonists. Fig. 1. Human TAS2R46 in complex with its agonist strychnine bound in the vestibular (a) and orthosteric (b) binding sites, together with the 2D structure of the strychnine agonist (c). In panels A and B, strychnine is shown in cyan and the regions of the protein treated at different resolutions in our MM/CG simulations are displayed with different colors: in yellow, the atomistic (MM) part of the receptor and, in blue, the coarse-grained (CG) part (see the Materials and methods section). A water droplet (in red) surrounds the extracellular part of the protein, in order to explicitly account for the hydration of the binding site. The two horizontal black lines delimit the transmembrane part of the receptor or the position of the lipid bilayer, mimicked by two planar walls in our MM/CG setup. For the sake of clarity, the two hemispheric walls setting a limit on the extracellular and intracellular ends of the protein, as well as the wall following the initial shape of the interface between protein and membrane, are not shown.
Nonetheless, further in silico and wet lab experiments are necessary to confirm whether the two-site architecture is present across the whole hTAS2R family.

hTAS2R38 in complex with its agonists propylthiouracil (PROP) and phenylthiocarbamide (PTC)
hTAS2R38 is a receptor involved in bitter taste perception in the tongue, as well as other extra-oral roles, such as anti-microbial response in the sinonasal cavity [10,15] and activation of transcription factors in pancreatic tumor cells [100]. Supporting the putative ectopic roles of hTAS2R38, different hTAS2R38 polymorphisms have been associated with several pathologies, such as predisposition to chronic rhinosinusitis [101], risk of dental caries [102], alteration of alcohol intake [103], alteration of body mass index [104], and ingestive behavior in women [105]. hTAS2R38 is a chemical group-specific bitter taste receptor [77], since it detects bitter agonists containing an isothiocyanate or thiourea group.
Here we investigate the receptor in complex with two typical agonists, PTC and PROP, by MM/CG simulations (Fig. 3a). The calculations turned out to be consistent with functional data for nine mutants [74,75]. In particular, we observed that the Asn103 sidechain forms a H-bond with both ligands (see Fig. 3b). This is consistent with the experimental data showing that Asn103Ala, Asn103Val and Asn103Asp mutations result in EC 50 larger values than the WT for both agonists: the first two mutations impair the formation of the H-bond, whereas the presence of the Asp in position 103 causes a repulsive electrostatic interaction with the partially negatively charged sulfur atom of the two ligands (see Fig. 3b). The simulations also show that Ser259 is in close proximity to the ligand without any direct interaction. Therefore, the larger EC 50 value of Ser259Val mutant compared to the WT is probably due to the presence of a bulkier residue that could hinder binding, rather than to the loss of a H-bond. Indeed, mutation of Ser259 into Ala, a residue similar in size, maintains EC 50 values similar to the WT. The EC 50 values of Trp99Ala, Trp99Val and Met100Ala are similar to those of WT for both agonists, while those of the Met100Val mutant are larger than the WT. The simulations suggest that Trp99 and Met100 do not interact directly with the ligand, though they are located close to the binding pocket and, thus, when Met100 is mutated into a branched amino acid, Val, it may occlude the binding site.
In conclusion, our MM/CG simulations results on hTAS2R38 are consistent with more than 20 mutagenesis data. These predictions would have been impossible to achieve using the bioinformatics approach only. In particular, the poses predicted by bioinformatics lack key H-bond and pep stacking ligand/protein interactions. This points to the relevance of molecular dynamics simulations for the structural refinement of these receptors' models. The fact that our simulations were not able to capture the vestibular binding site in hTAS2R38 may imply either that only one cavity is needed for the less promiscuous hTAS2R38 receptor or that more simulations are needed.

Conclusion
Our MM/CG-based predictions provide a rather detailed description of hTAS2R46-and hTAS2R38-agonist interactions, consistent with mutagenesis data [74e76]. They also allow us to hypothesize that hTAS2R46 features a two-site architecture, with an orthosteric and a vestibular binding site, similar to what has been already suggested for other members of the class A GPCRs [26,28,88e95]. The existence of a second binding site may be crucial to recognize the wide variety of hTAS2R46 agonists, by providing a two-step authentication mechanism for this promiscuous receptor. In contrast, the vestibular site was not captured by our simulations of hTAS2R38, perhaps because it is not required for a more selective receptor [77]. Nonetheless, a conservation analysis of the binding residues across the whole hTAS2R family suggests that this two-site architecture might also be present in other hTAS2Rs. Therefore, further simulations and mutagenesis studies are necessary to clarify this point. A water droplet (in red) surrounds the extracellular part of the protein, in order to explicitly account for the hydration of the binding site. The two horizontal black lines delimit the transmembrane part of the receptor or the position of the lipid bilayer, mimicked by two planar walls in our MM/CG setup. For the sake of clarity, the two hemispheric walls setting a limit on the extracellular and intracellular ends of the protein, as well as the wall following the initial shape of the interface between protein and membrane, are not shown. (b) 2D structures of PROP (top) and PTC (bottom).