Basal Histamine H4 Receptor Activation: Agonist Mimicry by the Diphenylalanine Motif.

Histamine H4 receptor (H4 R) orthologues are G-protein-coupled receptors (GPCRs) that exhibit species-dependent basal activity. In contrast to the basally inactive mouse H4 R (mH4 R), human H4 R (hH4 R) shows a high degree of basal activity. We have performed long-timescale molecular dynamics simulations and rigidity analyses on wild-type hH4 R, the experimentally characterized hH4 R variants S179M, F169V, F169V+S179M, F168A, and on mH4 R to investigate the molecular nature of the differential basal activity. H4 R variant-dependent differences between essential motifs of GPCR activation and structural stabilities correlate with experimentally determined basal activities and provide a molecular explanation for the differences in basal activation. Strikingly, during the MD simulations, F16945.55 dips into the orthosteric binding pocket only in the case of hH4 R, thus adopting the role of an agonist and contributing to the stabilization of the active state. The results shed new light on the molecular mechanism of basal H4 R activation that are of importance for other GPCRs.

Introduction

GPCR crystal structures, coupled with molecular dynamics (MD) simulations, site‐directed mutagenesis, and NMR studies, have led to increased knowledge about the activation mechanism of GPCRs. Agonists that bind to the orthosteric binding site located between transmembrane helices (TM) III, V, VI and VII favor the active receptor state, eliciting the maximal possible response and inducing G‐protein‐ and/or β‐arrestin‐mediated signaling. Inverse agonists exert the opposite effect to agonists and favor the inactive receptor state. Neutral antagonists do not affect the basal activity. However, little is known, both on the pharmacological and molecular levels, about basal (constitutive) GPCR activation, which describes the activity of a GPCR in the absence of a ligand. The existence of basal activity suggests that conformational transitions leading to the active receptor state may occur spontaneously. Such dynamic oscillations between active and inactive conformations have also been posited as a mechanism for partial activation, suggesting a mechanistic link between ligand‐independent and ‐dependent GPCR activation. Alternatively, a partially active receptor may adopt a static conformation that is intermediate between active and inactive. Furthermore, single‐molecule studies have indicated that differential ligand dwell‐times and/or influences on nucleotide exchange rates in the G‐protein add further layers of complexity to GPCR activation. Understanding the mechanisms and pathways of GPCR activation will help design new GPCR‐targeted drugs with tailored pharmacological responses and fewer side effects.

The histamine H4 receptor (H4R), a class A GPCR expressed mainly on immune cells, plays a fundamental role in processes such as cytokine release and chemotaxis. H4R offers a unique opportunity towards understanding the mechanisms and pathways of GPCR activation in a ligand‐independent manner. H4R orthologues represent ideal naturally occurring GPCRs with different degrees of basal activity. In contrast to mouse H4R (mH4R), human H4R (hH4R) shows a high level of basal activity (Figure 1).[ , ] In a previous molecular pharmacological study, we characterized the basal activity of a series of hH4R variants, comprising hH4R‐S179M, hH4R‐F169V, hH4R‐F169V+S179M, and hH4R‐F168A (Figure 1).[ , ] While the hH4R‐S179M variant exhibits a basal activity similar to that of hH4R, replacing F16945.55, situated in extracellular loop 2 (ECL2, numbering scheme according to GPCRdb ), by valine significantly decreased the basal activity. The basal activity of the hH4R‐F169V+S179M and hH4R‐F168A variants is even comparable to that of mH4R. These observations thus identified residues that account for the high basal activity of hH4R (Figure 1). However, the molecular mechanisms of basal hH4R activation, and how the equilibrium between inactive and active receptor states is shifted towards the inactive state in the other H4R variants, are still unknown.

Figure 1

Amino acids that act as key players in basal H4R activation. Fundamental amino acids for basal H4R activation (shown as spheres) were previously determined by site‐directed mutagenesis studies.[ , ] The respective basal activities of hH4R, mH4R, and the hH4R variants containing mutations in either F16845.54, F16945.55, and/or S1795.43 are depicted in the upper box. Further amino acids found to be involved in basal H4R activation are shown as ball and sticks. Amino acid numbering corresponds to hH4R.

Here, we performed long‐timescale MD simulations (2 μs) of each H4R variant to monitor conformational transitions in the ECL2 region, the orthosteric binding pocket, the transmission region, and the G‐protein binding site (Figure 1). Most notably, hH4R showed the most pronounced binding‐pocket contraction and TM VI outward movement, which reflects its highest basal activity. To complement the analyses of structural dynamics, we performed rigidity analyses to probe the structural stability of key H4R elements involved in GPCR activation. Our results highlight differential activation patterns of the H4R variants in different regions shown in Figure 1 and, hence, provide a molecular explanation for basal activity.

Results and Discussion

Dipping of F16945.55 into the orthosteric binding pocket

ECL2 has been proposed to be involved in ligand binding, selectivity, recognition, and to act as a gatekeeper for ligand entry. However, little is known about its contribution to basal GPCR activity. Among others, the M2 muscarinic receptor (M2R), the β2‐adrenergic receptor (β2AR), the histamine H3 receptor and hH4R, which all show basal activity, share the same F45.54 –F45.55 motif in ECL2 (F168–F169 for hH4R). Furthermore, 15 other GPCRs feature analogous F–Y, Y–F or Y– motifs in place of F–F. Remarkably, our MD simulations of hH4R revealed that F16945.55 dips into the orthosteric binding pocket and mainly interacts with the surrounding hydrophobic and aromatic residues W903.28, L913.29, Y953.33, P16645.52, W1725.36, L1755.39, Y3196.51, L3266.58, Y3407.35 and F3447.39 (Figure 2, numbering of TMs I–VII according to the Ballesteros–Weinstein nomenclature ). The dipping of F16945.55 was observed in four out of seven (approximately 60 %) of our wild‐type hH4R simulations but not at all for other variants (simulations performed in singlet). Strikingly, F16945.55 occupies the same area of the binding pocket as assumed for the imidazole moiety of histamine. It may thus act as a surrogate ligand in the highly basally active hH4R. Interestingly, GPCR crystal structures show two different orientations of the diphenylalanine (F–Y in case of H1R) motif: the two amino acids either point to opposite (β2AR, H1R ) or the same (M2R ) directions compared with hH4R simulations. Hence, though our simulations suggest the dipping of F16945.55 into the binding pocket and the direction of F16845.54 towards the cytoplasm, the two amino acids could potentially swap their orientations on longer timescales.

Figure 2

H4R variant‐dependent locations and orientations of key amino acids for basal H4R activation: F/A16845.54, F/V16945.55, S/M1795.43 for a) hH4R and the hH4R variants b) S179M, c) F169V, d) F169V+S179M, e) F168A, and f) F17045.54, V17145.55, M1815.43 of mH4R. Wild‐type and mutated residues of pre‐aligned H4R variant structures (a–f), showing the respective predominating cluster, are illustrated as sticks and semi‐transparent surface representation. To enhance comparability, the hH4R cartoon is highlighted in red (semi‐transparent) for all H4R variants (a–f). Carbon atoms of aromatic and hydrophobic hH4R residues, capable of interacting with F169, are colored in gray (a).

Contraction of the orthosteric binding pocket

In all MD simulations, the binding‐pocket dimension between TM III and VII, represented by the distances between the Cα atoms of D943.32 and Q3477.42 (hH4R variants) or Q3497.42 (mH4R), was initially approximately 11 Å, close to the corresponding distance in the inactive hH1R structure used for homology modelling. This distance increases by a maximum of approximately 3 Å for the basally inactive hH4R variants F169V, F169V+S179M and F168A, and for mH4R (Figure 3 b and Table 1). In contrast, it remained nearly constant for the hH4R‐S179M variant and decreased by 2 to 3 Å for wild‐type hH4R, both of which exhibit high basal activity. This distance is approximately 11 Å in aminergic GPCR crystal structures (Table S1 in the Supporting Information), both with bound inverse agonist and antagonist ((10.9±0.2) Å; N=30) and bound agonist ((10.9±0.5) Å; N=20). Thus, in our simulations, TM III approaches TM VII more closely in apo‐hH4R (≈8 Å) than previously observed in any ligand‐bound GPCR or found for the other H4R variants. The decreased distance between TM III and TM VII enables D943.32, which is involved in ligand binding and receptor activation,[ , ] to form additional interactions to residues in TM VII. The simulations of hH4R, its variants S179M (basally active) and F169V+S179M (weakly basally active) exhibited hydrogen bonds between D943.32, R3417.36 and W3487.43 in addition to a water‐mediated contact to Q3477.42 (Figures 3 a, 4, and S1 a and c in the Supporting Information). By contrast, in all variants with lower basal activity and mH4R, D943.32 is involved in a hydrogen‐bond network with S682.57 (and with W903.28 in hH4R‐F169V+S179M) (Figures 3 d, 4, and S1 b–d in the Supporting Information). A network comparable to the latter has been described for the crystal structures of M1 and M4 receptor‐antagonist complexes, and may be an indicator of the inactive receptor conformation of H4R. In a further hydrogen‐bonding network, all H4R variants showed differential, highly variant‐specific interactions between the key residue[ , ] E1825.46 (E1845.46 in mH4R) and Y953.33 or T1785.42 (Figures 3 a, d, 4, and S1 in the Supporting Information). While E1825.46 established hydrogen‐bond contacts to both Y953.33 and T1785.42 in hH4R and its variants S179M and F169V+S179M, it only displayed a single interaction to T1785.42 in hH4R‐F169V, and Y953.33 in hH4R‐F168A and mH4R. Consequently, the results emphasize that the key amino acids D3.32 and E5.46, which are usually involved in ligand‐dependent H4R activation, formed intra‐receptor interactions in the absence of a ligand. Depending on the amino acid composition, H4R variant‐specific active receptor conformations are stabilized to a greater or lesser extent.

Figure 3

Structural differences within the orthosteric binding pocket that relate to differences in basal activity. Insights into a) hH4R and d) mH4R binding pockets, represented by the corresponding predominant cluster (for cluster sizes see Table S2 in the Supporting Information). Residues involved in a hydrogen‐bond network with D943.32 (a,d) and E1825.46 (a) or E1845.46 (d) in at least one H4R variant are highlighted in dark purple and green, respectively. W3166.48 (a) or W3186.48 (d), and the semi‐transparent surface are shaded in red (a) or magenta (d). b) Time‐evolution of distances between Cα atoms of D943.32 and Q3477.42 (hH4R variants) or Q3497.42 (mH4R). c) Time course of W3166.48 (hH4R variants) or W3186.48 (mH4R) χ 2 rotational angles. For comparison, the corresponding distances of the inactive state hH1R (b) and the χ 2 torsional angles of both the inactive state M2R and the active state β2AR (c) are illustrated with straight lines. Colors referring to the H4R variants in (b,c) are specified in the legend.

Table 1

Analysis of selected distances and torsional angles (mean ± standard deviation) of wild‐type hH4R, its four variants, and mH4R.[a]

Parameter	hH4R	hH4R‐S179M	hH4R‐F169V	hH4R‐F169V+S179M	hH4R‐F168A	mH4R
distance[b] (Cα) D943.32–Q3477.42 [Å]	9.5±0.7	10.7±0.5	10.8±0.8	11.2±1.2	12.1±0.5	12.7±0.8
torsional angle[c] for χ 2 W3166.48 [°]	117.3±13.6	89.5±14.9	69.7±14.7	85.3±10.7	94.5±11.4	71.5±26.4
torsional angle[d] for χ 2 Y3587.53 [°]	102.8±21.1	87.3±31.4	129.7±30.3	126.9±20.2	105.8±20.4	139.1±20.7
distance[e] (Cα) R1123.50–A2986.30 [Å]	10.1±0.6	9.1±0.6	9.4±0.4	9.3±0.5	9.5±0.4	8.5±0.5

[a] Amino acid names and numbers correspond to hH4R. [b] Inactive‐state hH1R : 10.9, inactive‐state β2AR : 10.9, active‐state β2AR : 11.4 Å. [c] Inactive‐state hH1R : 105.9, inactive‐state M2R : 52.2, active‐state β2AR : 115.4°. [d] Inactive‐state hH1R : 128.7, inactive‐state β2AR : 138.4, active‐state β2AR : 86.3°. For torsion angle calculation, a periodicity of 180° was considered. [e] Inactive‐state hH1R : 8.0, inactive‐state β2AR : 11.1, active‐state β2AR : 17.2 Å.

Figure 4

Hydrogen‐bond (straight and water‐mediated) analysis of the MD simulation trajectories performed for the orthosteric binding pocket, the D612.50 cluster (transmission region), and the R1123.50 cluster (G‐protein binding site). The % occupancy over the simulation time values represent cumulated values taking into account all possible hydrogen‐bond donor and acceptor atom‐to‐atom combinations of the given amino acid pairs (excluding the backbone). The amino acid names and numbers given on the y axis correspond to hH4R.

The rotamer toggle switch contributes to basal H4R activation

The binding pocket connects to structures closer to the cytoplasm, so that interactions formed within the binding pocket affect distal rearrangements of the receptor. In this context, the proposed W6.48 rotamer toggle switch has been suggested to drive the P6.50 kink, and thus the outward movement of TM VI at the intracellular side. Analysis of rotameric W6.48 states of aminergic GPCR crystal structures (Table S1 in the Supporting Information) resulted in χ 2 torsional angles ranging from 29.7 to 113.7° for antagonist/inverse agonist‐bound structures ((94.7±29.5)°, N=30), from 80.1 to 86.1° for ligand‐free structures ((82.2±2.8)°, N=3) and from 97.3 to 124.4° for agonist‐bound structures ((112.4±6.5)°, N=20). In the hH4R simulations, this angle ((117.3±13.6)°) was comparable to that of active‐state aminergic structures such as, for example, β2AR (115.4°) (Table 1, Figure 3 c). Simulations of both mH4R and hH4R‐F169V gave values ((71.5±26.4) and (69.7±14.7)°) comparable to those of inactive‐state aminergic structures such as, for example, M2R (52.2°). The corresponding angles of the other hH4R variants ranged between these two values. These findings thus support the hypothesis that particularly apo‐hH4R features more highly populated active states than the other H4R variants. Hence, although the rotameric toggle switch has been discussed controversially in literature,[ , ] our results suggest that it is involved in basal H4R activation.

Structural rearrangements in the transmission region

D2.50, whose fundamental role for agonist‐dependent GPCR activation and allosteric sodium ion binding has been demonstrated, is located in the central core of the receptor. Hydrogen bonds between D612.50 and S1013.39 were more pronounced for both the highly basally active hH4R and its hH4R‐S179M variant than for the other H4R variants (see Figures 4, 5 a, b, and S2 in the Supporting Information; for further interactions formed by D2.50, please also see Figure S3a, Supporting Information). Consequently, the hydrogen‐bond network of the D2.50 cluster is expanded, and thus, TM III could undergo concerted helix movements with TM I, II and VII. The fact that active‐state crystal structures,[ , , ] unlike inactive[ , , ] ones, usually show this hydrogen bond confirms its important role for basal GPCR activation. Y3587.53, which is part of the NPxxY motif and lies adjacent to D2.50, acts as a further key player in GPCR activation. It conducts the effect of the tyrosine toggle switch after the so‐called hydrophobic barrier has been opened and thus, the water‐mediated hydrogen bonding network is expanded towards the DRY motif. Consistently with our interpretation of the simulations, Y3587.53 (Y3607.53 in mH4R) revealed overall χ 2 dihedral angles comparable to those of the inactive and active state β2AR (Table 1, Figure 5 c) in mH4R and hH4R, respectively. The mean χ 2 angles increased in the order hH4R‐S179M<hH4R≈hH4R‐F168A< hH4R‐F169V+S179M≈hH4R‐F169V<mH4R (Table 1). Consequently, in hH4R and its variants S179M and F168A, Y7.53 prevalently occupies a conformation that corresponds to a (basally) activated receptor state. In contrast, Y7.53 conformations of the hH4R variants F169V and F169V+S179M are comparable to that of mH4R and thus reflect an inactive receptor state. Driven by such conformational changes, Y7.53 only formed a prevalent hydrogen bond with D612.50 in mH4R (Figures 4, 5 a, b and S2 in the Supporting Information). Water‐mediated interactions, however, were present both in the other H4R variants (Figure S3 a, Supporting Information) and identified in both inactive and active aminergic GPCR crystal structures.[ , ] This result further supports the hypothesis that, of the H4R variants investigated, the inactive receptor state is most predominant in mH4R.

Figure 5

Structural differences within the transmission region that relate to differences in basal activity. Predominating cluster shown for a) hH4R and b) mH4R. Carbon atoms and the semi‐transparent surface of D612.50 are shaded in dark purple, and carbon atoms of the other residues are colored red (a) or magenta (b) (semi‐transparent surface shown for Y3587.53 (a) and Y3607.53 (b)). c) Time course of Y3587.53 (hH4R variants) and Y3607.53 (mH4R) χ 2 torsion angles (for color code, see legend). The χ 2 torsional angles of both inactive[15a]‐ and active[15b]‐state β2AR are shown with straight lines for comparison.

Outward movement of TM VI at the intracellular side

R3.50 has been proposed to play a key role in G‐protein activation. It forms the ionic lock with E6.30 that has been found in many GPCRs such as, for example, β2AR or H1R, and was suggested to play a role in the stabilization of the inactive receptor state. In contrast, H4R orthologues cannot form the ionic lock because E6.30 is replaced by non‐acidic amino acids (hH4R: Ala, mH4R: Gly).[ , ] The simulations suggest that R1123.50 forms more prevalent hydrogen‐bond contacts with D1113.49 in the H4R‐S5.43 variants (hH4R, hH4R‐F169V, hH4R‐F168A) than H4R‐M5.43 variants (hH4R‐S179M, hH4R‐F169V+S179M, mH4R) (Figures 4, 6 a, b, and S4 in the Supporting Information). However, only the latter showed hydrogen bonds between R3.50 and S6.36 that connect TM III and TM VI; these may replace the ionic‐lock interaction and therefore stabilize the receptor in the inactive conformation. Consequently, the substitution of S5.43 by Met leads to structural rearrangements that diminish the hydrogen‐bond formation between D3.49 and R3.50 and simultaneously enhance contacts between R3.50 and S6.36 (see Figures 4, 6 a, b, and S4 in the Supporting Information; for further interactions formed by D3.49, see Figure S3 b, Supporting Information). This result is in agreement with site‐directed mutagenesis studies, which show that the missing E6.30 alone is not the sole reason for the basal activity of H4R.

Figure 6

Structural differences within the G‐protein binding site that relate to differences in basal activity. Predominating cluster shown for a) the hH4R and b) the mH4R. Carbon atoms and the semi‐transparent surface of D1113.49 and R1123.50 are colored in dark purple. Carbon atoms of other residues and the respective cartoon are colored in red (a) or magenta (b). c) Time course of distances between R1123.50 and A2986.30 (hH4R variants) or G3006.30 (mH4R) (for color code of H4R variants, see legend). The corresponding distance of the inactive state hH1R is displayed with a straight line for comparison.

Depending on the amino acid composition of the H4R variant, all aforementioned motifs in concert contribute more or less to helical rearrangements at the intracellular face. In particular, TM VI is bent around P6.50 as part of the FxxCWxP motif, and thus forms a kink at this position. As a result, TM VI moves significantly outward upon receptor activation, enabling the binding of G‐proteins into the gap thus formed, and downstream G‐protein and/or β‐arrestin‐mediated signaling. As a measure of GPCR activation, initial mean Cα –Cα distances of around 8.0 Å (as in the hH1R template) between R1123.50 and A2986.30 (G3006.30 in mH4R) increased by a maximum of approximately 3 Å for WT‐hH4R during the last quarter of the simulation (Figure 6 c). In all other H4R variants the distance increase does not exceed 1.5 Å. This is consistent with the highest hH4R basal activity of all the H4R variants investigated, though not valid for hH4R‐S179M with the second highest basal activity. Comparison of the R3.50–E6.30 distance of inactive state β2AR (11.1 Å) with that of active state β2AR (17.2 Å) (see Table 1) revealed an outward movement of 6 Å during activation, suggesting that hH4R (≈3 Å) switched from an inactive to a partially active receptor conformation during MD simulation. Such an intermediate receptor conformation has also been suggested by Dror et al.

Mechanism of low‐level activation

To both gain insights into the mechanism of low‐level activation and to validate the reproducibility of our simulation results on basal hH4R activation, we performed six additional replica MD simulations of wild‐type hH4R (Figure 7). During the last 500 ns of the in total seven hH4R simulations, the dipping of F16945.55 into the binding pocket was observed in four simulations (1–4). Strikingly, the two of those four simulations that showed the most predominant dipping of F16945.55 into the binding pocket (>70 %) and either a binding pocket contraction (1) or only a subtle dilatation (2) also showed the most pronounced TM VI outward movement (>2.5 Å). By contrast, the other two simulations that either revealed lower dipping occupancies of F16945.55 (4) or higher binding pocket dilatations (3) showed a markedly reduced outward movement of TM VI (<2.0 Å). The three remaining simulations that showed no dipping of F16945.55 (5–7) resembled the H4R variants with a lower degree of basal activity (S179M, F169V, F169V+S179M or F168A) with respect to binding pocket contraction (<0.2 Å) and TM VI outward movement (0.9–1.3 Å). We thus observed two different activation states in the seven hH4R MD simulations, suggesting that hH4R switches between on (1, 2) and off (3–7) states in favor of adopting partially active conformations.

Figure 7

Activation patterns observed in MD simulations of hH4R (seven replica) and other H4R variants during the last 500 ns. The % dipping of F16945.55 into the binding pocket (blue, left y axis), measured as the occupancy of z distances <3.5 Å between the CZ atom of F16945.55 (CB in the case of Val) of the respective 500 simulation frames and the center of mass of the imidazole ring (non‐hydrogen atoms) of histamine (aligned docking pose). Contraction (+)/ dilatation (−) of the orthosteric binding pocket (dark magenta, right y axis) and TM VI outward (+)/inward (−) movement (green, right y axis), respectively, measured as the mean difference ± standard deviation of the D3.32–Q7.42 or R3.50–A/G6.30 distances between hH1R (PDB ID: 3RZE ) and the 500 simulation frames.

This is fully consistent with the concept of basal activity since basally active receptors transition spontaneously between inactive and active‐like conformations in the absence of a ligand. Note, however, that although the dipping and TM VI outward movement were reproducible and predominant for hH4R, we cannot entirely exclude their occurrences for other variants on a statistical basis. Indeed, such structural changes would have been expected for the still highly basally active hH4R‐S179M variant.

Shift in structural stability occurs with changing basal activity

To complement the analyses of differences in structural dynamics, we used Constraint Network Analysis (CNA) to compare the structural stability of WT‐hH4R, its four variants S179M, F169V, F169V+S179M, F168A, and mH4R. In CNA, biomolecules are modeled as a constraint network, in which atoms are represented as nodes, and covalent and non‐covalent binding interactions as constraints. The constraint network is efficiently decomposed into rigid clusters and connecting flexible hinge regions by applying rigidity theory. The inherent long‐range aspect of rigidity percolation (i.e., whether a region is flexible or rigid depends on distant structural details) makes rigidity analysis an attractive tool for studying changes in structural stability due to distant influences. Biomolecules display a hierarchy of structural stability that reflects the modularity of their structure. To identify this hierarchy, a “constraint dilution trajectory” of network states, obtained from an initial network topology by successively removing non‐covalent constraints, is analyzed (see Methods).[ , ] Figure 8 a shows a succession of five network states of hH4R and mH4R along such a constraint dilution trajectory, starting from the homology models (see Methods). In particular, regions known to be important in GPCR activation such as helices TM V and TM VI, become flexible earlier in hH4R than in mH4R. This qualitative analysis suggests that the lower structural stability of hH4R favors adaptations of these helices that allow the G‐protein to bind and, hence, that hH4R can be basally active, in contrast to mH4R.

Figure 8

Comparison of structural stability between hH4R and mH4R as well as four hH4R variants. a) Rigid cluster decomposition along the constraint dilution trajectory of hH4R (top) and mH4R (bottom). The five states show differences in the decay of structural stability in either species. At each state, all hydrogen bonds with an energy E HB>E cut are removed from the network. The asterisks indicate the position of TM VI. b) The heat map depicts ΔE values [Eq. (5)] of the four hH4R variants and mH4R relative to hH4R. A structural stabilization of either hH4R or hH4R variant/mH4R is color coded in blue and red, respectively (see color scale). The line plot shows the Spearman's rank correlation ρ between the computed, relative ΔE values and the rank of the basal activity of each system. Orange regions depict significant rank correlations ρ (p<0.05). c) Regions enclosing the residues F542.43–P702.59, W903.28–T923.30, and S1093.47–D1113.49 (V1845.48–V1905.54, A3066.38–G3106.42, and V3146.46–A3176.49) with a significant correlation between increasing (decreasing) structural stability and increasing levels of basal activity are mapped on the structure of H4R in blue (red).

To obtain more detailed information about structural stability in hH4R, the four hH4R variants and mH4R, we applied CNA to conformational ensembles of the last 500 ns of the 2 μs MD trajectories. We then computed the difference of the residue‐wise chemical potential energy [Eq. 1], which reflects changes in biomolecular stability between the system being considered and hH4R (Figures 8 b, S8, and S9 in the Supporting Information).

ΔE CNA has been shown to be a robust local stability measure for predicting the thermal stability of proteins. From the ΔE values, we identified regions in the receptors that correlate with the rank‐ordering of the systems with respect to basal activity (Figure 8 b): The structural stability of the systems increases with increasing levels of basal activity for the regions enclosing residues F542.43 –P702.59 in helix TM II as well as W903.28 –T923.30 and S1093.47 –D1113.49 in helix TM III. In contrast, the structural stability of the systems decreases with increasing levels of basal activity for residues V1845.48 –V1905.54 in TM V as well as A3066.38 –G3106.42 and V3146.46 –A3176.49 in TM VI. Notably, the analyses thus reveal focused, sequentially continuous but topologically distant residue clusters in the receptors in which changes in structural stability correlate with basal activity, rather than a broad distribution of stability changes across the structures.

Furthermore, the analyses show that, while regions critical for and known to move during GPCR activation (TM V and TM VI) become less structurally stable with increasing levels of basal activity, regions that are rather immobile during GPCR activation (TM II and TM III) become more structurally stable. Therefore, a local shift in structural stability within the GPCR occurs with changing basal activity, rather than a global, uniform change, so that overall structural stability is maintained in each case.

Interestingly, analyzing residues in TM VI adjacent to the G‐protein binding site (L2966.28 –F3126.44), which undergoes the largest structural changes during activation, in terms of the sum of their ΔE values yields a measure that significantly and very well correlates (R 2=0.94, p=0.001) with the level of basal activity (cf. Figure 9; a significant and good correlation is also found if the entire, 2 μs‐long trajectories are used for analysis (Figure S10 in the Supporting Information); likewise, a significant and very good (R 2=0.94, p=0.001) correlation is found if the last 500 ns of six additional replica MD simulations for hH4R are evaluated (Figures 9 a and S11 in the Supporting Information); both findings demonstrate the robustness of the CNA results): The hH4R‐S179M variant shows only a minor stabilization of −0.80 kcal mol−1 compared to hH4R, which is in agreement with the very similar basal activities of the two systems. A stronger stabilization is found for the hH4R‐F169V variant (−9.48 kcal mol−1), and the strongest effect (−14.20, −20.94, and −17.57 kcal mol−1, respectively) for hH4R‐F169V+S179M, hH4R‐F168A, and mH4R. The findings are again in agreement with hH4R‐F169V+S179M and hH4R‐F168A showing a decrease in basal activity comparable to that of mH4R.[ , ]

Figure 9

Correlation between structural stability of TM VI and experimental basal activity. a) The violin plot shows the distribution of structural stability in TM VI enclosing the residues L2966.28–F3126.44. The relative change in structural stability (ΔE CNA) for each species compared to hH4R is shown on the top. The statistical significance between pairs of mean E CNA values (dashed line) in increasing order is *: p<0.01, **: p<0.001, ***: p<0.005. For hH4R, open circles depict the mean ΔE CNA values of six additional replica MD simulations. b) The structural stability of TM VI correlates significantly and very well with experimental basal activity. Additionally, the mean ΔE CNA value over all seven hH4R replicas is shown, relative to the ΔE CNA value of the one hH4R trajectory shown as violin plot in Figure 9 a; error bars denote the SEM. The basal activity in percentage (0 (basally inactive) to 100 % (basally active)) was transformed into the logit scale.

Predictive value of the parameters

Table 2 summarizes further geometrical features observed in the simulations in relation to the basal activity of the six H4R variants. Although the binding pocket contraction is not as clear for the hH4R‐S179M mutant as for the wild‐type receptor, the D3.32 –Q7.42 distance is quite predictive. Both direct and water‐mediated hydrogen bonds between TM III and VII indicate activation, although they are also found for the inactive hH4R‐F169V+S179M variant. The hydrogen bond between D943.32 and S682.57 or W903.28, on the other hand, is only found in the inactive variants and mH4R, making its absence indicative of activation. The remaining hydrogen‐bond networks in the binding site and the W6.48 toggle switch are not diagnostic. More distally located motifs such as the D2.50 –S3.39 hydrogen bond, the Y7.53 toggle switch, and the R3.50 –A6.30 distance are quite predictive, but not the D2.50 –Y7.53 hydrogen bond. As an alternative measure, structural stabilities of the lower part of TM VI (close to the G‐protein binding site) can accurately predict basal activities (Figure 9 b).

Table 2

Summary[a] of the impact of wild‐type hH4R, its four variants, and mH4R[b] on GPCR‐activating motifs.

Parameter	hH4R	S179M	F169V	F169V+S179M	F168A	mH4R
basal activity[c] relative to hH4R [%]	100	86	45	20	0	0
binding pocket
dipping of F16945.55 into the binding pocket	√	×	×	×	×	×
distance (Cα) D943.32–Q3477.42 [%]	0	38	39	54	80	100
TM III–VII direct hydrogen bonds	√	√	×	√	×	×
TM III‐VII water‐mediated hydrogen bonds	√	√	×	√	×	×
hydrogen bond D943.32–S682.57 (or W903.28)	×	×	√	√	√	√
hydrogen bond Y953.33–E1825.46 [%]	97	24	0	79	100	64
hydrogen bond T1785.42–E1825.46	24	37	90	100	0	4
Y953.33, T1785.42, E1825.46 hydrogen‐bond network	√	√	×	√	×	×
torsional χ 2 angle for W3166.48 [%]	100	42	0	33	52	4
central region
hydrogen bond D612.50–S1013.39 [%]	100	58	0	2	1	19
hydrogen bond D612.50–Y3587.53	×	×	×	×	×	√
torsional χ 2 angle for Y3587.53 [%]	30	0	82	77	36	100
G‐protein binding site
hydrogen bond R1123.50–S3046.36 [%]	0	45	2	13	0	100
distance (Cα) R1123.50–A2986.30 [%]	100	32	57	48	61	0
structural stability L2966.28–F3126.44 [%]	0	4	45	68	100	84

[a] Effects of H4R variants on GPCR activating motifs are either illustrated with ticks (applicable) and crosses (not applicable) or in a range between 0 % (lowest value) and 100 % (highest value). The percentage values were calculated from mean distances, torsional angles, or from occupancy over the simulation time values (hydrogen bonds) as outlined in the Computational methods section and Table S3 in the Supporting Information. [b] Amino acid naming and numbering correspond to hH4R. [c] Basal activities (hH4R=100, mH4R=0 %) were calculated from experimentally determined[ , ] intrinsic activities of thioperamide.

√

×

×

×

×

×

√

√

×

√

×

×

√

√

×

√

×

×

×

×

√

√

√

√

√

√

×

√

×

×

×

×

×

×

×

√

Conclusion

Prior molecular pharmacological (site‐directed mutagenesis) studies have demonstrated the importance of the diphenylalanine (F16845.54 –F16945.55) motif for basal H4R activation: whereas hH4R‐F169V showed intermediate basal activity, the basal activities of both hH4R‐F169V+S179M and hH4R‐F168A were even comparable to that of the basally inactive mH4R. Therefore, the question arose, by which molecular mechanism the diphenylalanine motif contributes to basal H4R activation.

Our simulations unveiled the dipping of F16945.55 into the orthosteric binding pocket and its conformational stabilization by the neighboring F16845.54 as the trigger for high basal activity. The dipping was observed in approximately 60 % of the hH4R simulations but not at all for other variants (Table 2). hH4R, therefore, appears to be unique among the receptor variants studied. Strikingly, our hH4R simulations showed an overlap of the pharmacophoric groups of F16945.55 with those of hH4R ligands such as histamine. F16945.55 apparently resembles a ligand with an intrinsic activity ranging between that of a neutral antagonist and a full agonist. Depending on the ability of F16945.55 to interact with the surrounding hydrophobic and aromatic amino acids and, thus, to contract the binding pocket, TM VI at the intracellular side moves outward to a greater or lesser extent (Table 2). In this respect, hH4R prefers switching between on and off states rather than stabilizing partially active conformations. Basal activity thus increases to the extent that the basal equilibrium between inactive and active states is shifted towards the active state. An agonist binding to such a pre‐activated receptor can subsequently only shift an intermediate to an active receptor state. In the absence of basal activity, by contrast, an agonist could shift an inactive to an active state.

Complementarily, results from rigidity analysis provide an excellent measure for describing the degree of basal activity in H4R. First, we observed focused, sequentially continuous but topologically distant residue clusters in H4R with structural stability characteristics that correlate significantly either with increasing or decreasing system's basal activity. Second, the structural stability of TM VI (L2966.28 –F3126.44) alone, known to be important for GPCR activation, is an excellent indicator for the degree of a system's basal activity. Based on the almost perfect correlation between the computed structural stability and basal activity, a single‐parameter model can be derived to efficiently probe the potential basal activity of new H4R variants. Furthermore, this model might be transferable to probe for the efficacy of H4R ligands.

This work and its results pose many questions such as whether the results are solely applicable for H4R or more generally for class A GPCRs. Recently, the diphenylalanine motif gained high interest due to its central role for molecular recognition in peptide‐based supramolecular assemblies and regulatory influences on bacterial morphology. Notably, the class A GPCRs known to exhibit the most pronounced basal activity (H4R, H3R and β2AR) all contain a diphenylalanine motif in ECL2. Hence, it is tempting to speculate that the diphenylalanine motif exerts a general functional influence in which it drives the self‐regulatory effect in transmembrane receptors, and thus, the system's basal activity. Furthermore, the exact nature of basal receptor activation, i.e., the structural comparability to agonist‐induced activation, remains unknown. All these questions are the subject of further investigations.

Computational Methods

Homology modelling

A previous hH4R homology model, based on the inactive state hH1R (PDB ID: 3RZE) was used. It includes all extracellular (ECL) and intracellular (ICL) loops. Models of the other H4R variants (hH4R variants S179M, F169V, F169V+S179M, F168A, and mH4R) were prepared with the modeling suite SYBYL‐X 1.3 (Tripos Inc., St. Louis, MO USA) using the hH4R homology model as template. For the preparation of the mH4R model, the ECL2 loop upstream of C16445.50 was remodeled using the loop search module within SYBYL‐X 1.3. All models showed the disulfide bond between C873.25 and C16445.50 (hH4R variants) or C16645.50 (mH4R). In place of ICL3, eight alanine residues were introduced. All residues were simulated in their dominant protonation state at pH 7. The N and C termini were positively (NH3 +) and negatively (COO−) charged.

Molecular dynamics simulations

Topologies and coordinates of the solvated (explicit water) and ionized (net charge of the entire system was zero) wild‐type or mutant hH4R models, and the mH4R model, were prepared using the pdb2gmx module within GROMACS 4.5. The systems were submitted to both steepest‐descent and conjugated gradient energy minimization (without restraints) to remove bad van der Waals contacts of the amino acid side chains. The solvated GPCRs were then equilibrated in the NPT ensemble for 10 ns, applying harmonic restraints of 1000 kJ mol−1 nm−2 to protein main‐chain atoms.

For membrane simulations, a pre‐equilibrated dioleoylphosphatidylcholine (DOPC) lipid bilayer containing 72 DOPC lipids was used. This system was enlarged in the x and y directions, resulting in a lipid bilayer consisting of 253 DOPC molecules. This extended system was subjected to energy minimization (no restraints) and successively, equilibration runs were performed in the NVT (100 ps; weak harmonic restraints of 1000 kJ mol−1 nm−2 applied on DOPC atoms) and NPT (10 ns; no restraints) ensembles.

Each receptor was inserted into this augmented, fully hydrated and equilibrated DOPC lipid bilayer using g_membed. The net charge of the simulation systems was neutralized by addition of 5 (hH4R variants) and 11 (mH4R) Cl− ions. The systems additionally contained around 217 DOPC and 28,000 water molecules, accounting for a total of around 120,000 atoms in a box spanning approximately 90 Å×91 Å×143 Å. Consecutively, these systems were energy minimized (without restraints) and equilibrated in the NPT ensemble for 10 ns (weak harmonic restraints of 1000 kJ mol−1 nm−2 applied to protein main‐chain atoms). Production molecular‐dynamics (MD) simulations (2 μs) of each H4R variant were started from the final frame of the equilibration runs. Data were collected every 100 ps. Although in other cases GPCR activation/deactivation typically demands much longer than 2 μs of simulation[ , ] or enhanced sampling, we usually observe an induction period of approximately 500 ns with the AMBER force field, suggesting that 2 μs simulations are adequate for GPCR simulations in a membrane. The nature of basal activity might furthermore contribute to lower timescales needed for basal compared to ligand‐dependent GPCR activation.

All minimization and MD simulation steps were carried out with GROMACS 4.5 (GROMACS 2016/2018 in the case of six replicate hH4R simulations). Initial velocities of all MD simulation runs were randomly assigned to the atoms. The SPC/E water model, the Amber ff99SB protein force field, and for DOPC molecules the GAFF force field[ , ] were employed. A temperature of 310 K was maintained using a temperature bath and a time constant of 0.1 ps. Protein, DOPC, and water (including Cl− ions) atoms were coupled separately. The Berendsen barostat with a time constant of 5 ps and a surface tension of 22 dyn cm−1 was applied to maintain a pressure of 1 bar and to ensure membrane properties in agreement with experiment. The volume compressibility was chosen to be 4.5×10−5 bar−1. The nonbonded cutoff was set to 10 Å and long‐range electrostatic interactions were computed using the Particle Mesh Ewald (PME) method with an interpolation order of 4, an Ewald tolerance of 1×10−5 and FFT grid spacing of 1.2 Å. Periodic boundary conditions were applied in all directions. A cutoff of 14 Å was used for short‐range van der Waals interactions. Bonds involving hydrogen atoms were constrained using LINCS, enabling a time step of 2 fs.

Induced‐fit docking

“Flexible” docking of histamine to the hH4R model was essentially performed as described for muscarinic acetylcholine receptors in Pegoli et al. by using the Protein Preparation Wizard, and the LigPrep and Induced‐fit docking modules (Schrödinger LLC, Portland, OR, USA).

Data analysis

Data were analyzed by using GROMACS 4.5 (g_cluster, g_density, g_rms, g_rmsf) and GROMACS 5 (gmx distance) analysis tools every nanosecond. The gromos clustering method was applied for cluster analysis of the entire (2 μs) trajectories, setting an RMSD cutoff of 2.5 Å. Hydrogen bonds and torsional angles were analyzed using the CPPPTRAJ module of Amber 14 (University of California, San Francisco, CA USA). Hydrogen bond and xy plots (distance, torsion) were visualized by using the programing language R and the packages devEM,F plot3d, Plotrix, Peptides and CircStats. For the purpose of comparison, distances and torsional angles of aminergic GPCRs (Table S1 in the Supporting Information) were calculated using the UCSF Chimera package, version 1.11.2, and are given as mean ± standard deviation. % values given in Table 2 were calculated as follows [Eq. 2] from mean distances, torsional angles (Table 1), or occupancy over the simulation time values in the case of hydrogen bonds (Figure 4). For minimum and maximum values, please see Table S3 in the Supporting Information. Figure 7 was prepared with the Prism 5.01 software (GraphPad, San Diego, CA, USA), and all molecular figures with PyMOL Molecular Graphics system, version 1.8.2.1 (Schrödinger LLC, Portland, OR USA).

Rigidity analysis

Rigidity analyzes were performed using the CNA software package. CNA was applied on ensembles of network topologies generated from conformational ensembles. Each network of covalent and non‐covalent (hydrogen bonds including salt bridges and hydrophobic tethers) interactions was constructed with the FIRST (Floppy Inclusions and Rigid Substructure Topography) software (version 6.2), to which CNA is a front‐ and back‐end. The strengths of hydrogen bonds (including salt bridges) were assigned by the energy E HB computed by FIRST. Hydrophobic interactions between carbon or sulfur atoms were taken into account if the distance between these atoms was less than the sum of their van der Waals radii (C: 1.7 Å, S: 1.8 Å) plus D cut=0.25 Å. In order to elucidate the hierarchy of structural stability in a biomolecule, a “constraint dilution trajectory” of network states {σ} was analyzed, which was generated by successively removing hydrogen bond constraints (including salt bridges) in the order of increasing strength.[ , ] Thus, only those hydrogen bonds are retained in a network of state σ that have an energy E HB≤E cut(σ). Altered biomolecular stability along a constraint dilution trajectory was quantified based on neighbor stability maps (rc with i, j being residue numbers) [Eq. 3].

Here, only short‐range rigid contacts were considered that have at least one pair of heavy atoms of the residue pair R {, A {, separated by a distance ≤4.5 Å. A rigid contact rc between pairs of residues ceases to exist when both residues stop sharing the same rigid cluster c of a set of rigid cluster C (Figure S7, Supporting Information). The double sum [Eq. 4]

then represents the chemical potential energy (ECNA) due to non‐covalent bonding, obtained from the coarse‐grained, residue‐wise network representation of the underlying biomolecular structure.[ , ] A per‐residue decomposition of Equation (4), which is the chemical potential energy of residue i obtained by summation over all n short‐range rigid contacts the residue is involved in (Figure S7 in the Supporting Information), is computed according to Equation 5:

Conformational ensembles used as input for CNA were extracted from MD trajectories of 2 μs length starting from hH4R, the four hH4R variants S179M, F169V, F169V+S179M, and F168A, as well as mH4R. The first 1.5 μs of each trajectory were considered equilibration phase of the system, and conformational ensembles for CNA were extracted from the last 500 ns. To estimate the uncertainty in the CNA computations, we split the last 500 ns of each trajectory into five individual sets of 100 ns length. Because the correlation time for decay of fluctuations of E CNA [Eq. (4)] for all trajectories is about 0.1 ns, the five individual sets are statistically independent, and the standard error of the mean (SEM) is calculated from the standard deviation (SD) of the five mean E values [Eq. (5)] according to Equation 6:

To probe for a drift in the E values over the simulation time, we computed E for intervals of 100 ns along the whole trajectory of 2 μs length for residues in TM VI, lying next to the orthosteric binding site and the G‐protein binding site; these values show no correlation with simulation time, indicating that there is no drift (Figure S12, Supporting Information).

To probe for the robustness of the E values for residues in TM VI, lying next to the orthosteric binding site and the G‐protein binding site, we also analyzed entire, 2 μs‐long MD trajectories (Figure S10 in the Supporting Information) as well as the last 500 ns of six additional replica MD simulations for hH4R (Figures 9 and S11 in the Supporting Information).

Statistical analysis

For each system, the average E values and the corresponding SEMs from Eq. 6 were calculated using NumPy. According to the law of error propagation, the total SEM for differences in structural stability between wild‐type hH4R and the four hH4R variants or mH4R was computed according to Equation 7:

We used SciPy for performing the following statistical analyses: (I) A two‐sided Welch's t‐test to probe if ΔE values are significantly different from zero (p<0.05) and if independent samples of E values have significantly different average values; (II) Pearson's and Spearman's rank correlation between ΔE values and (the rank of) the basal activity of each hH4R variant and mH4R. For computing the Pearson correlation, we summed ΔE values for the residues in TM VI (L2966.28 –F3126.44) and then averaged over the last 500 ns for each species. The basal activity scale, ranging from 0 % (inactive) to 100 % (active), was transformed into the logit scale, with an adjustment of 2.5 % for the lower and upper end of the scale. For computing the rank correlation, we calculated a running average over the ΔE values with respect to i, with a windows size of ten and a step size of one.

Conflict of interest

The authors declare no conflict of interest.

As a service to our authors and readers, this journal provides supporting information supplied by the authors. Such materials are peer reviewed and may be re‐organized for online delivery, but are not copy‐edited or typeset. Technical support issues arising from supporting information (other than missing files) should be addressed to the authors.

Supplementary

Click here for additional data file.