Understanding Ligand Binding to G-Protein Coupled Receptors Using Multiscale Simulations

Human G-protein coupled receptors (GPCRs) convey a wide variety of extracellular signals inside the cell and they are one of the main targets for pharmaceutical intervention. Rational drug design requires structural information on these receptors; however, the number of experimental structures is scarce. This gap can be filled by computational models, based on homology modeling and docking techniques. Nonetheless, the low sequence identity across GPCRs and the chemical diversity of their ligands may limit the quality of these models and hence refinement using molecular dynamics simulations is recommended. This is the case for olfactory and bitter taste receptors, which constitute the first and third largest GPCR groups and show sequence identities with the available GPCR templates below 20%. We have developed a molecular dynamics approach, based on the combination of molecular mechanics and coarse grained (MM/CG), tailored to study ligand binding in GPCRs. This approach has been applied so far to bitter taste receptor complexes, showing significant predictive power. The protein/ligand interactions observed in the simulations were consistent with extensive mutagenesis and functional data. Moreover, the simulations predicted several binding residues not previously tested, which were subsequently verified by carrying out additional experiments. Comparison of the simulations of two bitter taste receptors with different ligand selectivity also provided some insights into the binding determinants of bitter taste receptors. Although the MM/CG approach has been applied so far to a limited number of GPCR/ligand complexes, the excellent agreement of the computational models with the mutagenesis and functional data supports the applicability of this method to other GPCRs for which experimental structures are missing. This is particularly important for the challenging case of GPCRs with low sequence identity with available templates, for which molecular docking shows limited predictive power.

Human G-protein coupled receptors (GPCRs) convey a wide variety of extracellular signals inside the cell and they are one of the main targets for pharmaceutical intervention. Rational drug design requires structural information on these receptors; however, the number of experimental structures is scarce. This gap can be filled by computational models, based on homology modeling and docking techniques. Nonetheless, the low sequence identity across GPCRs and the chemical diversity of their ligands may limit the quality of these models and hence refinement using molecular dynamics simulations is recommended. This is the case for olfactory and bitter taste receptors, which constitute the first and third largest GPCR groups and show sequence identities with the available GPCR templates below 20%. We have developed a molecular dynamics approach, based on the combination of molecular mechanics and coarse grained (MM/CG), tailored to study ligand binding in GPCRs. This approach has been applied so far to bitter taste receptor complexes, showing significant predictive power. The protein/ligand interactions observed in the simulations were consistent with extensive mutagenesis and functional data. Moreover, the simulations predicted several binding residues not previously tested, which were subsequently verified by carrying out additional experiments. Comparison of the simulations of two bitter taste receptors with different ligand selectivity also provided some insights into the binding determinants of bitter taste receptors. Although the MM/CG approach has been applied so far to a limited number of GPCR/ligand complexes, the excellent agreement of the computational models with the mutagenesis and functional data supports the applicability of this method to other GPCRs for which experimental structures are missing. This is particularly important for the challenging case of GPCRs with low sequence identity with available templates, for which molecular docking shows limited predictive power.

INTRODUCTION
G-protein coupled receptors (GPCRs) are one of the largest protein superfamilies, with more than 800 (4%) genes in humans (Venter et al., 2001;Fredriksson et al., 2003;Lagerstrom and Schioth, 2008;Tikhonova and Fourmy, 2010). They detect a wide variety of extracellular signals (from photons to hormones and neurotransmitters) and trigger a myriad of intracellular transduction cascades (using different G-proteins and second messengers) (Alexander et al., 2017). These pleiotropic receptors are involved in many physiological functions, from vision to chemical sensing and neurotransmission, and, hence, they are attractive targets for pharmaceutical intervention. Approximately 34% of currently FDA-approved drugs bind to GPCRs (Hauser et al., 2018) and they are used to treat disorders as diverse as pain, hypertension, diabetes, cancer or neurological diseases . Given the physiological and pharmacological relevance of GPCRs, unraveling their ligand binding determinants can be extremely useful both for understanding receptor function and for designing new drugs.
Based on phylogenetic and sequence conservation analyses, GPCRs can be classified in 5 different families or classes (Fredriksson et al., 2003;Schioth and Fredriksson, 2005): rhodopsin (class A), secretin (class B1), adhesion (class B2), glutamate (class C), and frizzled/taste2 (class F). Nonetheless, taste 2 (or bitter taste) receptors have also been proposed to form part of class A (Nordstrom et al., 2011) or even constitute a sixth, additional family (class T) (Munk et al., 2016a). Since the appearance of the first crystal structure of rhodopsin in 2000, experimental structural characterization of GPCRs is blossoming (Munk et al., 2019). As of February 2019, there are 59 unique receptor structures solved (https:// gpcrdb.org/structure/statistics), most of them corresponding to the rhodopsin (or class A) family (Figure 1). Molecular dynamics (MD) simulations started from these experimental structures have provided very important insights into ligand binding and receptor activation (Miao and McCammon, 2016;Sengupta et al., 2016;Latorraca et al., 2017;Marino and Filizola, 2018;Torrens-Fontanals et al., 2018;Velgy et al., 2018).
Nonetheless, the experimental structural coverage is still very far from the total of 800 GPCRs. In particular, there are no experimental structures available for three receptor groups: olfactory receptors (ORs, which constitute half of class A), taste 2 receptors (TAS2Rs, which represent the third largest GPCR family) and adhesion (class B2) receptors. In silico modeling can help to fill this gap of ∼87% structurally uncharacterized GPCRs (Pándy-Szekeres et al., 2017). Indeed, the communitywide GPCR Dock assessment (Michino et al., 2009: Kufareva et al., 2011 has shown that homology modeling and ligand docking are able to provide valuable information on receptorligand interactions, in particular for those GPCR targets that have templates with sequence identity higher than 35-40% (Kufareva et al., 2011;Beuming and Sherman, 2012). Subsequent refinement of the bioinformatics-based models through molecular dynamics simulations (Kufareva et al., 2014;Cavasotto and Palomba, 2015;Lupala et al., 2018) and integration of experimental (mutagenesis FIGURE 1 | GPCR statistics. The number of members is based on reference (Munk et al., 2016a) and the number of experimental structures was taken from the GPCRdb database (https://gpcrdb.org/structure/statistics, accessed on January 2019). GPCRs are grouped according to the class A-F nomenclature (Fredriksson et al., 2003). Within class A, two groups are differentiated: non-olfactory and olfactory receptors. The taste 2 receptors have been proposed to belong to either class A (Nordstrom et al., 2011) or class F (Fredriksson et al., 2003), or even constitute a novel, sixth class (Munk et al., 2016a). Legend labels for those groups without experimental structures are in gray. and ligand structure-activity relationship) data (Munk et al., 2016b) further increases the model quality to values close to experimental accuracy. However, approximately half of GPCRs do not have a close template (i.e., an experimental structure of a receptor from the same family with a similar ligand). For instance, the sequence identity of 90% GPCRs with the rhodopsin template (representative of the largest GPCR family, class A) is lower than 20% (Zhang et al., 2006). Therefore, in most cases the in silico modeling approach needs further improvement, typically using molecular dynamics (Kufareva et al., 2014;Cavasotto and Palomba, 2015;Lupala et al., 2018).
Chemosensory receptors (olfactory and bitter taste receptors) are among the GPCRs without close templates. Increasing evidence shows that these receptors are expressed not only in the nose and the tongue, but also in other parts of the body (Foster et al., 2014;Abaffy, 2015;Ferrer et al., 2016;Shaik et al., 2016;Lu et al., 2017;Behrens and Meyerhof, 2019;Lee et al., 2019) and thus they have become attractive novel targets for drug design campaigns (Lee et al., 2019). However, chemosensory receptors represent a major challenge for computational modeling. Their sequence identity with the available GPCR templates is lower than 20%  and thus only low resolution homology models can be generated (Kufareva et al., 2011;Fierro et al., 2017). Hence, our lab has made a major effort to attempt at improving such low resolution homology models and at making valuable predictions of the ligand binding determinants of these receptors. In this review, we first explain the computational approach used in our group to study low resolution GPCR models, based on the combination of stateof-the-art bioinformatics techniques and multiscale molecular dynamics simulations, as well as its validation on a class A GPCR (the β2-adrenergic receptor) with a solved crystallographic structure. Then, we show that, although bioinformatics-based models can be a good starting point to study receptor-ligand interactions, multiscale simulations significantly improve the quality of the models for which MM/CG simulations have been run so far. A perspective on this multiscale approach concludes this review.

Bioinformatics
Given the lack of experimental structures, the initial structures of the receptor/ligand complexes are generated using bioinformatics approaches. Although there are several webservers specialized in GPCR modeling (Launay et al., 2012;Zhang et al., 2015;Busato and Giorgetti, 2016;Esguerra et al., 2016;Pándy-Szekeres et al., 2017;Worth et al., 2017;Miszta et al., 2018), here we used the GOMoDo webserver (Sandal et al., 2013), which combines in a single pipeline homology modeling and molecular docking for GPCRs.
Since the sequence identity of any given olfactory or bitter taste receptor with the available GPCR templates is lower than 20% , special care needs to be taken in the sequence alignment step. Hence, the alignment was done using profile Hidden Markov Models (HMMs) of the corresponding target receptor family and the GPCR template(s), which were generated with HHPred (Soding et al., 2005). This approach has been shown to improve the target-template alignment for distant homologs (Soding et al., 2005), in particular for GPCRs (Kufareva et al., 2014). This alignment was further improved by manual curation, taking advantage of the conserved seven transmembrane (7TM) helix topology and the presence of common conserved features across GPCRs (Lagerstrom and Schioth, 2008;Venkatakrishnan et al., 2013;Pydi et al., 2014Pydi et al., , 2016Tehan et al., 2014;de March et al., 2015;Di Pizio et al., 2016;Fierro et al., 2017). Moreover, since template selection is difficult with such low sequence identity, several models based on different templates were built using MODELLER (Webb and Sali, 2016), and the best model was selected considering also structural quality parameters (Melo et al., 2002;Shen and Sali, 2006). The receptor structural model thus generated was then submitted to molecular docking using HADDOCK (Dominguez et al., 2003). Although other docking approaches were tested [based on AutoDock Vina (Trott and Olson, 2010) or Glide (Friesner et al., 2004)], no significant improvement in the quality of the models was observed . The location of the ligand binding pocket inside the 7TM bundle is conserved (Venkatakrishnan et al., 2013), despite the low sequence identity among GPCRs. Moreover, the results of the GPCR Dock competitions (Michino et al., 2009;Kufareva et al., 2011Kufareva et al., , 2014Cavasotto and Palomba, 2015;Munk et al., 2016b) seem to indicate that incorporating information about putative binding residues (from experimental data or computational predictions) helps to improve the docking results. Therefore, an informationdriven approach was taken, in which the computationally predicted binding residues [using fpocket (Le Guilloux et al., 2009)] were used to guide the docking. Nonetheless, the fine details of the ligand binding site are expected to be highly variable across GPCRs (Venkatakrishnan et al., 2013), due to the chemical diversity of the GPCR ligands. Hence, in our HADDOCK-based docking approach both receptor and ligand were considered fully flexible in order to allow mutual readjustments. Other flexible docking approaches have also been successfully employed by other groups to predict the binding determinants of chemosensory receptors [see for instance (Di Pizio and Niv, 2014;Di Pizio et al., 2017;Xue et al., 2018)].

Multiscale Molecular Dynamics Simulations
The results of the GPCR Dock competitions [reviewed in references (Cavasotto and Palomba, 2015) and (Ranganathan et al., 2017)] showed that refinement of the docked complexes using molecular dynamics simulations can significantly improve the prediction of receptor/ligand interactions. This is particularly important for GPCR models based on low sequence identity, as it is the case for chemosensory receptors, where the low accuracy of the side chain prediction and the limited sampling of the docking algorithms may undermine the quality of the bioinformaticsbased models. There are several studies in the literature applying molecular dynamics simulations to chemosensory receptors (Gelis et al., 2012;Lai and Crasto, 2012;Charlier et al., 2013;Lai et al., 2014;Chen et al., 2018;Jaggupilli et al., 2018;Liu et al., 2018;Bushdid et al., 2019). Here we focus on a hybrid, multiscale approach developed in our group (Neri et al., 2005(Neri et al., , 2008Leguèbe et al., 2012;Giorgetti and Carloni, 2014;Musiani et al., 2015;Tarenzi et al., 2017), which is tailored to study ligand binding in GPCRs.
As shown in Figure 2, the ligand, the surrounding protein residues (typically the extracellular half of the receptor) and water molecules are described with molecular mechanics (MM) using the GROMOS united-atom force field (Schuler and Van Gunsteren, 2000;Schuler et al., 2001;Oostenbrink et al., 2004). Instead, the rest of the protein (i.e., the intracellular half of the receptor) is treated at the coarse grained (CG) level using a Gō model (Go and Abe, 1981). Each amino acid is mapped into a single coarse grained bead corresponding to the alpha carbon atom and native contacts are mimicked by introducing two new potential terms. The bonded interactions between consecutive CG beads are taken into account using a quartic potential, whereas the non-bonded interactions between non-consecutive CG beads are described through a Morse-like potential. The MM and CG regions are connected by an interface (I) region, which ensures the continuity of the protein backbone by coupling the two levels of resolution. The MM-I interaction is treated at the atomistic level using the GROMOS force field, whereas the I-CG interaction is described using the Gō model. Namely, bonded interactions are calculated between the Cα atoms of the I residues and the consecutive CG beads, whereas non-bonded interactions are computed using both the Cα and Cβ atoms of the I residues and the non-consecutive CG beads.
The presence of the lipid bilayer is modeled implicitly, using three wall potentials: a "coating surface" wall that simulates the FIGURE 2 | Simulation setup of the Molecular Mechanics/Coarse-Grained (MM/CG) approach. The receptor (GPCR) is divided in three regions: the extracellular MM part (in blue), the intracellular CG part (in yellow) and the connecting interface (I, in green). The ligand (shown in red) and the surrounding receptor residues and water molecules (in blue and white, respectively) are described with MM. The presence of the lipid bilayer is modeled implicitly by incorporating three different wall potentials (upper and lower membrane planes and coating surface) and two additional hemispheric walls are included to cap the ends of the protein and prevent water evaporation. effect of the lipid hydrophobic tails embracing the protein surface and two "membrane plane" walls that mimic the presence of the lipid head groups. In addition, two "hemispheric" wall potentials are included to cap the extracellular and cytoplasmic ends of the protein and to prevent water evaporation. Water molecules, Cα atoms and aromatic residues Phe, Trp, and Tyr (the so-called "anchor residues") are affected by these boundary potentials, which are added to the MM/CG potential energy function as functions of the distance of an atom to the closest wall. Recently, a reservoir of CG water has been introduced around the MM water cap, permitting water molecules to freely diffuse between the MM and CG regions, changing on the fly their resolution. This allows to carry out simulations in a statistically well-defined (grand canonical) ensemble in the higher-resolution MM region, resulting in a further improved description of the binding poses and the binding site flexibility (Tarenzi et al., 2019).
Compared to docking, these multiscale simulations allow to (i) sample protein flexibility and protein/ligand interactions more extensively (∼1 µs timescale) and (ii) include explicit water molecules, which may be involved in ligand binding in GPCRs (Pardo et al., 2007;Angel et al., 2009;Venkatakrishnan et al., 2019). Moreover, the use of the Go model in the intracellular half of the receptor prevents possible unfolding problems due the initial wrong orientation of the side chains in the low resolution homology model. For further details on the MM/CG implementation, we refer the reader to some recent reviews (Giorgetti and Carloni, 2014;Musiani et al., 2015;Schneider et al., 2018).

Validation of the Molecular Mechanics/Coarse Grained (MM/CG) Approach
The reliability of the MM/CG approach was assessed using the β2-adrenergic receptor (β2-AR) in complex with either its inverse agonist S-carazolol or its agonist R-isoprenaline (Leguèbe et al., 2012;Marchiori et al., 2013). The availability of a crystal structure of the receptor (for the first complex) (Cherezov et al., 2007), as well all-atom (AA) molecular dynamics simulations (for both complexes) (Vanni et al., 2011) allows to compare the results of the MM/CG simulations with both static and dynamical data. Three different types of tests were carried out (Leguèbe et al., 2012;Marchiori et al., 2013), started from different initial structures: (i) the same initial structures of the β2-AR/Scarazolol and β2-AR/R-isoprenaline complexes as the atomistic simulations, (ii) an alternative initial structure of the β2-AR/Scarazolol complex built by displacing the ligand to a position where none of the crystallographic receptor/ligand interactions was present, and (iii) a low resolution model of the β2-AR/Scarazolol complex built using bioinformatics. Each of the test simulations were ∼0.8 µs long.
The first test (Leguèbe et al., 2012) showed that the MM/CG approach is able to preserve the receptor/ligand complex structure observed in the crystal structure, as well as to provide dynamical and hydration information similar to the AA simulations, but at a lower computational cost. Complementarily, the second test (Leguèbe et al., 2012) confirmed that the agreement between the MM/CG and AA simulations observed in the first test is not due to the use of a common initial structure and, furthermore, demonstrated the predictive power of the MM/CG approach even when starting from a wrong binding pose. Nonetheless, the two previous tests can be considered as redocking experiments: even though the system was converted from AA into hybrid MM/CG [test (i)] or the ligand was moved out of place [test (ii)], the binding residues were already positioned as in the correct binding pose. Instead, the third test (Marchiori et al., 2013) validated the reliability of the MM/CG approach applied to low resolution models, as the ones used for the bitter taste receptors discussed in the next section. In such models, the orientation of the side chains is expected to be hardly accurate, due to the low sequence identity with the template used in the homology modeling (Chothia and Lesk, 1986;Baker and Sali, 2001;Eramian et al., 2008;Olivella et al., 2013;Piccoli et al., 2013;Busato and Giorgetti, 2016). Indeed, the homology model of the β2-adrenergic receptor (built using as template the experimental structure of squid rhodopsin) shares only a 20% sequence identity with the target and thus docking of the ligand S-carazolol resulted in a wrong binding pose. However, the ∼0.8 µs MM/CG simulation is able to yield a binding pose showing receptor/ligand interactions similar to the crystallographic ones (Marchiori et al., 2013).

Predictive Power of the Computational Models of Chemosensory Receptors
In order to investigate the performance of bioinformatics and multiscale simulations in predicting receptor/ligand interactions in chemosensory receptors, we selected those receptor/ligand pairs for which experimental data are available . As of August 2017, these included seven olfactory receptor/odorant complexes and fifteen bitter taste receptor/bitter tastant complexes with available sitedirected mutagenesis data and functional assays, typically agonist dose-response curves. The docked receptor/ligand complexes were obtained using the bioinformatics protocol described in the Materials and Methods section, whereas the three MM/CG simulations analyzed (for the complexes TAS2R38/6-n-propylthiouracil, TAS2R38/phenylthiocarbamide and TAS2R46/strychnine) were taken from previous studies from our group (Biarnés et al., 2010;Marchiori et al., 2013;Sandal et al., 2015).
In order to compare the computational models with the experimental data, we defined "computational binding" and "computational non-binding" residues, as well as "experimental binding" and "experimental non-binding" residues . Computational binding and non-binding residues were determined based on the receptor/ligand distance (using a cutoff of 5.5 Å) and the presence or absence (respectively) of an actual chemical interaction (such as hydrogen bonds, salt bridges, hydrophobic or aromatic interactions, etc.). Experimental binding residues were inferred from experiments based on (i) the effect of the mutation on the half maximal effective concentration (EC 50 ) value and (ii) their position in the upper extracellular part of the receptor, where the canonical binding site of class A GPCRs is located (Venkatakrishnan et al., 2013). Residues whose mutation does not change EC 50 and/or that are located in the lower intracellular part of the receptor are considered as experimental non-binding residues. Obviously, this is a simplistic definition of binding residue, as from the experimental data we cannot discard that these residues might also be involved in receptor activation (see reference Fierro et al., 2017 for further discussion).
Comparison of the computational and experimental residues yielded four different possible test outcomes. "True positives" (TP) were amino acids identified as binding residues by both experiment and computation, "false positives" (FP) were amino acids identified as non-binding residues by experiment but as binding residues in computation, "true negatives" (TN) were amino acids identified as non-binding residues by both experiment and computation, and "false negatives" (FN) were amino acids identified as binding residues by experiment but not in computation. In order to assess the agreement of the computational models with the experimental data, two statistical parameters, precision (PREC) and recall (REC), were calculated: These parameters are close to 1 when the computational predictions were consistent with the experimental data, and zero when they were not. Precision and recall values were calculated for the best docking poses of the twenty-two complexes investigated and for a representative snapshot of each of the three MM/CG simulations analyzed .
We found that the predictive power of the bioinformatics approach varied from complex to complex. Nonetheless, the general agreement between the binding residues identified in the docking poses and those inferred from experiments was low, with only 36% of the predictions consistent with experiment . Residues shown experimentally to be important for binding were not observed in the docked complexes (i.e., low recall) and/or residues not involved in protein/ligand interactions according to experiments were predicted as binding residues by computation (i.e., low precision). Most likely, this is due, among other factors, to the low sequence identity (<20%) between the chemosensory receptor targets and the available GPCR templates, as well as the limited sampling of the docking algorithms. Therefore, although the bioinformaticsbased models are a good starting point to study ligand binding determinants in chemosensory receptors, they appear to require further refinement . This finding is consistent with the results of the GPCR Dock competitions, which indicated that models based on sequence identity below 30% need substantial improvement in order to reach accuracy comparable to experimental structures (Kufareva et al., 2011(Kufareva et al., , 2014. Next, we compared the performance of molecular dynamics for the three bitter taste receptor complexes studied so far with (∼0.8-1 µs long) MM/CG simulations (Marchiori et al., 2013;Sandal et al., 2015). We found that the predictive power of the computational models improved dramatically, with 96-100% of the predictions in agreement with experiment . Most residues shown to be involved in binding by experiments are captured by the MM/CG simulations and the number of wrong predictions was minimal (i.e., both recall and precision increased to values near or equal to one, see Table 1).
Considering the nearly 20 mutants tested experimentally for either TAS2R38 (Biarnés et al., 2010: Marchiori et al., 2013 or TAS2R46 Born et al., 2013;Sandal et al., 2015), the agreement of the computational models with experiments seems really remarkable. Moreover, although in our analysis we used all the available mutagenesis data to validate a posteriori the MM/CG results, simulations were also able to predict new binding residues. Indeed, the simulations of the TAS2R38 and TAS2R46 complexes suggested several binding residues not tested previously and these predictions were subsequently verified by performing additional mutagenesis and functional assays (Marchiori et al., 2013;Sandal et al., 2015). Altogether, multiscale simulations seem to be a robust approach for identifying ligand binding residues in olfactory and 1 | Evaluation of the performance of the bioinformatics models and the multiscale MM/CG simulations for the three bitter taste receptors studied so far, i.e. TAS2R38/6-n-propylthiouracil (PROP), TAS2R38/phenylthiocarbamide (PTC) and TAS2R46/strychnine . bitter taste receptors, at least for the three bitter taste receptor complexes studied so far . This is consistent with the conclusions of the GPCR Dock competitions, where molecular dynamics simulations and integration of experimental data (such as site-directed mutagenesis or ligand structureactivity relationships) were shown to improve the computational predictions (Kufareva et al., 2011(Kufareva et al., , 2014Cavasotto and Palomba, 2015;Ranganathan et al., 2017). Nonetheless, the application of MM/CG simulations to other chemosensory receptor complexes with available experimental data is still needed to firmly establish the reliability and transferability of this method.

Insights Into Ligand Selectivity Determinants in Bitter Taste Receptors
There are around 1,000 compounds characterized as bitter, which vary significantly in size, polarity and chemical structure Behrens and Meyerhof, 2018;Dagan-Wiener et al., 2019). To make things even more puzzling, three receptors (TAS2R10, TAS2R14, and TAS2R46) out of the 25 bitter taste receptors are able to recognize about half of these bitter compounds (Behrens and Meyerhof, 2018). In contrast to this broad agonist spectrum, there are two receptors, TAS2R38 and TAS2R16, that are specialized in detecting a specific chemical group (thiourea/isothiocyanate and β-Dglucopyranoside, respectively) (Behrens and Meyerhof, 2018).
Here we discuss in detail the structural predictions described above to investigate whether they can help understand the molecular basis of this disparate ligand selectivity. In particular, the MM/CG approach has been applied so far to one receptor of each group, i.e., TAS2R46 (Sandal et al., 2015) and TAS2R38 (Biarnés et al., 2010;Marchiori et al., 2013). The microsecond-long simulations of TAS2R46 in complex with its agonist strychnine (Sandal et al., 2015) showed that the ligand can explore not only one but two different binding cavities (Figure 3). The first one coincides with the canonical binding site of class A GPCRs (i.e., the so-called orthosteric site), whereas the second is located further toward the extracellular side and thus has been denoted as "vestibular." The mutagenesis data is compatible with this two-site architecture, as the residues experimentally inferred to be involved in binding Born et al., 2013;Sandal et al., 2015) are distributed between the two sites (Figure 3). Moreover, the identified vestibular binding cavity overlaps with the extracellular allosteric binding site observed for class A GPCRs (Dror et al., 2011(Dror et al., , 2013Kruse et al., 2012;Abdul-Ridha et al., 2014;Latorraca et al., 2017;Thal et al., 2018), further supporting its existence. This two-step binding architecture may constitute the molecular basis of the "access control" mechanism proposed by Meyerhof and coworkers  and would help TAS2R46 to discriminate the wide range of ligands recognized by this promiscuous receptor (Sandal et al., 2015). Moreover, a bioinformatics analysis of the binding residues predicted for TAS2R46 across the bitter taste receptor family showed that half of these functionally relevant positions are conserved in two or more TAS2Rs, suggesting that the vestibular site might also be present in other receptors of this family (Sandal et al., 2015). However, the ∼0.8 µs simulations of TAS2R38 in complex with either PTC or PROP showed the ligand bound in a single site, corresponding to the orthosteric one (Marchiori et al., 2013). This hints at the possibility that the vestibular site is not as crucial for a group specific receptor such as TAS2R38 or even that the two-site architecture is not required for a more selective receptor . Naturally, given the crudeness of our models, further simulations and experimental studies on other members of the bitter taste receptor family are needed in order to confirm this proposal.

CONCLUSIONS
Given the scarcity of experimental structural data (Munk et al., 2019), computational modeling of GPCRs is essential to understand ligand binding and design new drugs targeting this biologically and pharmacologically relevant family (Michino et al., 2009;Kufareva et al., 2011Kufareva et al., , 2014Cavasotto and Palomba, 2015;Ranganathan et al., 2017;Lupala et al., 2018). These computational approaches (Figure 4) include homology modeling and molecular docking, often supplemented with experimental (mutagenesis and ligand structure-activity relationship) data. Subsequent refinement with molecular dynamics simulations has been shown to further improve the computational predictions (Kufareva et al., 2014;Cavasotto and Palomba, 2015;Ranganathan et al., 2017;Lupala et al., 2018). The accuracy of the models thus generated might reach values near the experimental ones for those GPCRs with a close structural template (i.e., with sequence identity larger than 35-40% and a chemically similar ligand) (Kufareva et al., 2011;Beuming and Sherman, 2012). However, for most GPCRs the closest structural template has sequence identity below this threshold, and thus computational predictions become challenging. This the case for olfactory and bitter taste receptors, which constitute the first and third largest GPCR groups, respectively, as their sequence identity with the available GPCR templates is below 20%.
In this review, we have shown that molecular dynamics simulations, in particular the multiscale molecular mechanics / coarse grained approach developed in our group (Neri et al., 2005(Neri et al., , 2008Leguèbe et al., 2012;Giorgetti and Carloni, 2014; FIGURE 3 | Two binding site architecture of TAS2R46. The agonist strychnine (in licorice representation) can bind in either the orthosteric site (left panel) or the vestibular site (right panel). Receptor residues interacting with the ligand in the orthosteric site, the vestibular site or both sites are shown with blue, yellow or green spheres, respectively. The central panel displays the distribution of the ligand center-of-mass z coordinate for the two ∼1 µs MM/CG simulations, showing that strychnine stabilized either in the orthosteric site or in a second (vestibular) site, located further toward the extracellular side.
FIGURE 4 | Proposed protocol to study ligand binding in low resolution GPCR models. The initial model generated with homology modeling and molecular docking (e.g., using the GOMoDo webserver, Sandal et al., 2013) is further refined using molecular dynamics (e.g., using the molecular mechanics/coarse grained or MM/CG approach developed in our group, Leguèbe et al., 2012). The receptor/ligand interactions predicted by the simulations have to be validated by extensive comparison with experimental data (typically site-directed mutagenesis and dose-response functional assays). Musiani et al., 2015;Tarenzi et al., 2017), can overcome, at least in part, these limitations  and successfully predict residues involved in ligand binding for the three bitter taste receptor complexes studied so far (Biarnés et al., 2010;Marchiori et al., 2013;Sandal et al., 2015). The natural extension of these previous works would be to other bitter taste and olfactory receptors for which experimental data are available. In addition, MM/CG simulations could be easily applied to other GPCRs. Although this approach has been used so far for a limited number of GPCR/ligand complexes (Leguèbe et al., 2012;Marchiori et al., 2013;Sandal et al., 2015), the excellent agreement of the computationally predicted binding poses with the experimental mutagenesis data [for the aforementioned three bitter taste receptor complexes (Marchiori et al., 2013;Sandal et al., 2015)] or the crystal structure [for the β2-adrenergic receptor (Leguèbe et al., 2012)] further supports the applicability of the MM/CG method to other GPCR/ligand complexes. Indeed, MM/CG simulations have been recently used to model the synthetic agonist diphenyleneiodonium chloride (DPI) bound to its target receptor GPR3 (Capaldi et al., 2018). Two of the predicted DPI binding residues were successfully validated a posteriori using mutagenesis and functional assays, as previously done for TAS2R38 (Marchiori et al., 2013) and TAS2R46 (Sandal et al., 2015).

DATA AVAILABILITY
The raw data supporting the conclusions of this manuscript will be made available by the authors, without undue reservation, to any qualified researcher.

AUTHOR CONTRIBUTIONS
All authors listed have made a substantial, direct and intellectual contribution to the work, and approved it for publication.

FUNDING
The authors acknowledge the Ernesto Illy Foundation (Trieste, Italy) for financial support. We are also grateful to the Jülich-Aachen Research Alliance High Performance Computing for computer time grants JARA0023, JARA0082, and JARA0165. We also thank the Central Library of Forschungszentrum Jülich for the open access publications fees.