Agonist Binding and G Protein Coupling in Histamine H2 Receptor: A Molecular Dynamics Study

The histamine H2 receptor (H2R) plays an important role in the regulation of gastric acid secretion. Therefore, it is a main drug target for the treatment of gastroesophageal reflux or peptic ulcer disease. However, there is as of yet no 3D-structural information available hampering a mechanistic understanding of H2R. Therefore, we created a model of the histamine-H2R-Gs complex based on the structure of the ternary complex of the β2-adrenoceptor and investigated the conformational stability of this active GPCR conformation. Since the physiologically relevant motions with respect to ligand binding and conformational changes of GPCRs can only partly be assessed on the timescale of conventional MD (cMD) simulations, we also applied metadynamics and Gaussian accelerated molecular dynamics (GaMD) simulations. A multiple walker metadynamics simulation in combination with cMD was applied for the determination of the histamine binding mode. The preferential binding pose detected is in good agreement with previous data from site directed mutagenesis and provides a basis for rational ligand design. Inspection of the H2R-Gs interface reveals a network of polar interactions that may contribute to H2R coupling selectivity. The cMD and GaMD simulations demonstrate that the active conformation is retained on a μs-timescale in the ternary histamine-H2R-Gs complex and in a truncated complex that contains only Gs helix α5 instead of the entire G protein. In contrast, histamine alone is unable to stabilize the active conformation, which is in line with previous studies of other GPCRs.


Introduction
Histamine is an important tissue hormone which is primarily detected by human cells via the four histamine receptors H 1 R, H 2 R, H 3 R and H 4 R [1,2]. All of them belong to the class A of G protein-coupled receptors (GPCRs) [3]. Binding of histamine favors the transition from inactive towards active subsets of receptor conformations which allows for the coupling to intracellular binding partners (IBPs) [4,5]. Each histamine receptor couples to a specific subset of IBPs including several types of G proteins and arrestins, which participate in a plethora of different pathways and thus lead to various changes within the behavior of the cell and the surrounding tissue [3]. While the H 1 R is especially relevant for inflammatory and allergic reactions [6], the main role of the H 2 R is regulation of gastrointestinal motility, intestinal and most notably gastric acid secretion [7]. Therefore, H 2 R antagonists are used in the treatment of gastroesophageal reflux or peptic ulcer disease [8,9].
Furthermore, H 2 R is expressed in myocardial cells [10]. Recent studies suggested that H 2 R antagonists may reduce the risk of heart failure (HF) [11]. Patients with HF that were subjected to treatment with H 2 R antagonists were shown to have a less alarming cardiac morphology and in general weaker symptoms [12]. Despite its pharmacologic importance, until now there is no 3D-structural information available for H 2 R hampering a mechanistic understanding and rational design of drugs to modulate H 2 R activity.
We have generated a computational model of the ternary histamine-H 2 R-G s complex and simulated its conformational dynamics. Both ligand binding and conformational changes of GPCRs frequently occur on timescales longer than microseconds [13,14] and can therefore only partly be assessed by conventional molecular dynamics (cMD) simulations. To enhance conformational sampling, we have applied metadynamics and Gaussian accelerated molecular dynamics (GaMD) simulations in this study.
Metadynamics is a method which accelerates sampling along one or more collective variables (CVs) by adding gaussian potentials to the free energy landscape at regular time intervals. These are centered at the evaluated CV values of the respective frame and aim to facilitate the exploration of so far unvisited configurations. This history-dependent potential reduces the probability of resampling the same position in CV space [15][16][17]. By using the multiple walker approach, which introduces several simultaneous simulations with shared potentials, metadynamics gets even more efficient. To date, metadynamics simulations have been successfully used to study GPCR-ligand interactions in a couple of cases including opioid [18,19], vasopressin [20], chemokine [21], or cannabinoid receptors [22].
In the present study, a multiple walker metadynamics simulation was applied to deduce the histamine binding mode in H 2 R, which was subsequently refined by 1 µs cMD simulations according to a previously established strategy [41]. GaMD was applied in addition to assess the role of different binding partners in the conformational stability of the active H 2 R. We compared dynamics of the histamine-H 2 R-G s ternary complex, the binary complex of histamine-H 2 R and the apo H 2 R without any binding partner. In addition, a truncated ternary complex that contains only the C-terminal G s -α5 helix instead of the entire heterotrimeric G s protein was simulated to investigate whether the α5 helix can sufficiently stabilize the active state. For the comparison of these four systems, a total of 8 cMD and 14 GaMD simulations (1 µs length each) were performed. Using the computational strategy described above, we were able to generate a model of the ternary histamine-H 2 R-G s complex and to characterize the H 2 R-histamine and H 2 R-G s interfaces in detail. In addition, simulations of H 2 R in complex with different binding partners provided important insights into structural dynamics of the active H 2 R and its conformational transitions in the absence of intracellular binding partners.

Modeling of the H 2 R in the Active State
To obtain a model of H 2 R in the active state, we included both an agonist (histamine) in the orthosteric pocket and a G protein at the cytosolic side. H 2 R efficiently couples to G s , but may also bind other types of G proteins [7,42]. Therefore, we used the crystal structure of the human β 2 adrenoceptor in complex with G s (PDB code 3SN6 [43]) as template to generate an H 2 R model in the active state (see Methods section for details). The resulting model ( Figure 1) still lacks the agonist histamine, for which the binding mode was determined in the second step by molecular dynamics simulations. For clarity, only the backbones are shown in ribbon representation. The H 2 R is colored in light blue, the G s α subunit is orange with the C-terminal α5 helix highlighted in red. The G s β and G s γ subunits are shown in green and violet, respectively. The camelid antibody fragment which was kept for stabilization is colored in white.
The simulation of histamine binding to the H 2 R was done using the previously established protocol from Söldner et al. [41] that combines multiple walker metadynamics simulations with a ligand clustering protocol and refinement by cMD. The free energy landscape derived from the respective metadynamics simulation is shown in Figure 2a. The shape of this curve is similar to those observed in previous ligand binding simulations [41,44] and indicates an energy minimum within the orthosteric pocket of the receptor. According to the strategy from Söldner [41] structures within 3 Å of the global minimum (red segment in Figure 2a) were grouped into five clusters. The representative structures from these clusters were all located in the orthosteric pocket but exhibited certain differences in their orientation and interactions in the pocket (Figure 2b). To assess the conformational stability and refine ligand binding mode of the H 2 R, subsequent 1-µs MD simulations were performed for the representative structures from all five clusters. According to the binding energy (Figure 2c), cluster5 exhibits the tightest interactions and the lowest fluctuations over simulation time. The high conformational stability of cluster5 also becomes apparent from inspection of the distance between the histamine ammonium group and the carboxyl group of D98 3.32 , which represents a key ionic anchor in aminergic GPCRs [45]. Again, cluster5 exhibits the tightest interaction and smallest fluctuations indicating the highest stability of this binding mode (Figure 2d). Based on the observations above, we considered cluster5 as the most favorable interaction mode and inspected the histamine-receptor interactions in more detail.
In addition to the salt bridge between the histamine ammonium group and D98 3.32 , polar interactions exist between one imidazole nitrogen and the sidechains of D186 5.54 and T190 5.461 (Figure 3a). This finding is in good agreement with data from site-directed mutagenesis: Mutation of the acidic residues (D98N or D186A) abolished both [ 3 H]-methyltiotidine binding and histamine-stimulated increases in cAMP content [46]. T190A or T190C mutations retained the ability to bind [ 3 H]-methyltiotidine, although with significantly reduced affinity [46]. The remaining interactions between histamine and H 2 R are mostly hydrophobic. Residues that form a large number of contacts over the simulation time include V99 3.33 , C102 3.36 , W247 6.48 , Y250 6.51 , F251 6.52 , and Y278 7.42 (Figure 3a,b). The receptor backbone is drawn as white ribbons. Histamine is displayed as sticks with the carbon atoms colored in green. All residues within 3 Å of the ligand are shown as sticks with black carbon atoms. (b) Number of contacts between histamine and individual H 2 R residues. The columns show per-frame averages for the refinement simulation of cluster5 whereas the error bars indicate the respective standard deviations. Only residues with an average of at least 1 contact per frame were taken into account.

Interface between the Receptor and the G s Protein
In addition to the histamine binding pocket analyzed above, the second key interface for GPCR function is formed on the intracellular side between the GPCR and the α-subunit of a G protein.
We analyzed the H 2 R-G s interactions formed in our modeled complex and compared it to the interactions present in the β 2 AR-G s template crystal structure used for modelling ( Figure 4).
In both GPCRs the major interaction sites are formed by the intracellular loops (ICL1-3) and cytosolic ends of helices H3, H5, and H6. The key interacting residues within these regions ( Figure 4) are frequently identical or at least exhibit conserved biophysical properties between H 2 R and β 2 AR. However, there are also differences in the G s interaction pattern between H 2 R and β 2 AR in the loop connecting helix H7 and H8. This loop is one residue longer in H 2 R compared to β 2 AR and residues L291-R293 form a significantly larger number of contacts than the respective sequence patch in β 2 AR ( Figure 4). This prompted us to compare the interaction of H 2 R and β 2 AR with G s in more detail. The focus was on the polar interactions, which generally play an important role for the specificity of protein-protein recognition [47].  . Sequence alignment of the histamine H 2 receptor (H2R) and the β 2 adrenoceptor (B2AR). The sequences are shaded according to the number of contacts (red > orange > yellow > green > blue) that the respective residues formed with the G s protein during the simulations (H 2 receptor: H-G s -cMD1, H-G s -cMD2) or in the crystal structure (β 2 adrenoceptor: PDB code 3SN6). The approximate length of the structural elements is indicated above the alignment. The ICL3 residues in lower case letters are not resolved in the β 2 -AR crystal structure. Figure 5 shows that both H 2 R and β 2 AR form numerous polar interactions with the C-terminal helix α5 of the G s -α-subunit. Both complexes contain a conserved salt bridge between D381 of the G s protein and R215 in H 2 R, or K232 in β 2 AR. H 2 R exhibits two additional salt bridges formed by R228 and R293 in H 2 R, which have no structural equivalent in the β 2 AR complex (Figure 5a). R293 recognizes two acidic residues of G s -namely D354 and E392 (Figure 5c). R228 also recognizes E392 and in addition forms a salt bridge with the charged carboxy-terminus of L394 (Figure 5c). Taken together R228 and R293 of H 2 R are part of a comprehensive network of polar interactions that involves D354, E392, and L394 of G s . These contacts may offer an explanation for the observation that H 2 R interacts with G s as a major signaling partner.
However, we want to emphasize that previous studies have shown that there is no simple GPCR sequence pattern explaining G protein coupling specificity [48,49] and that at least two additional structural factors may play an important role: Structural studies of GPCRs in complex with different G proteins indicate, that the position and dynamics of helix H6 significantly affect coupling selectivity [50][51][52].
(ii) As recently demonstrated for β 2 AR, transient interactions observed between the GPCR and GDP-bound G s protein may represent an intermediate on the way to the formation of the final complex and may contribute to coupling specificity [53]. Future studies of H 2 R in complex with different G proteins or investigations of different G protein binding modes will be required to assess the role of these structural factors for H 2 R coupling preferences.

Conformational Stability of the Active H 2 R
One aim of the present study was to assess the role of intra-and extracellular binding partners in conformational stability of the active H 2 R. In this context, we compared the dynamics of the histamine-H 2 R-G s ternary complex, the histamine-H 2 R binary complex, and apo H 2 R without any binding partner. In addition, a truncated ternary complex was generated that contains only the C-terminal G s -α5 helix (residues T369-L394) instead of the entire heterotrimeric G s protein to investigate whether the α5 helix can sufficiently stabilize the active H 2 R. To enhance conformational sampling and facilitate large structural rearrangements, Gaussian accelerated MD (GaMD) simulations [23] were performed in addition to cMD simulations on all systems investigated (see Methods section for the simulation details).
The most prominent motion detected in our simulations is the inward rotation of the intracellular end of helix 6 towards helix 3 (Figure 6a). This effect can be seen from the changes in the distance between the Cα atoms of R116 3.50 and T233 6.34 in simulations of the apo H 2 R (Figure 6b) and the H 2 R-histamine complex (Figure 6c). For the apo H 2 R, the motion may occur in less than 200 ns (see cMD2 and GaMD2 run), whereas in other simulations (cMD1, GaMD4) no approaching of the helices is observed on the timescales simulated (Figure 6b). This can be most likely attributed to the limited simulation time that does not allow a comprehensive sampling of slow motions, which are consequently only detected in a subset of the simulations. The simulations of the H 2 R-histamine complex (Figure 6c) qualitatively show the same behavior as the simulations of the apo H 2 R. In contrast, no decrease in H3-H6 distance is observed in the presence of G s or G s -α5 (Figure 6d,e). This data suggests that histamine alone in unable to stabilize the active H 2 R, whereas addition of the G s and G s -α5 can do so. The inward motion of the intracellular end of H6 is generally described as a hallmark of GPCR inactivation [13,54]. We investigated whether additional structural features of inactive H 2 R emerge during our simulations. One such feature is the formation of the "ionic lock" between the oppositely charged sidechains of E229 6.30 and R116 3.50 (Figure 7). Ionic lock formation is observed in simulations of apo H 2 R and histamine-bound H 2 R (Figure 7b,c), but not in complexes containing the G s or G s -α5 (Figure 7d,e). Comparison of the distances obtained from the individual simulations between Figure 6 and Figure 7 reveals that the inward motion of H6 generally correlates with ionic lock formation, i.e., the ionic lock appears to be the energetically most favorable sidechain arrangement when H3 and H6 are in close distance. This is remarkable, because the β 2 AR, which was used as a template for H 2 R modelling, does not display an ionic lock in the crystal structures of the inactive state (PDB code 5JQH [55]). However, in contrast to the crystal structure, long-timescale MD simulations of the β 2 AR show the ionic lock formation, which indicates the presence of a conformational equilibrium between conformations with the lock formed and the lock broken [56,57]. In summary, all analyses above indicate that the active conformation is retained in the ternary complex containing both histamine and the heterotrimeric G protein. In contrast, histamine alone is unable to stabilize this conformation. This is in line with previous work showing that the presence of the G protein is the key determinant for establishing an active conformation, whereas an agonist is insufficient to stabilize the active state [58,59]. In our study, the G s -α5 alone exerts a significant stabilization similar to that of the intact G protein. It is in line with the observation that mini G proteins or G protein mimicking nanobodies can stabilize the active conformation of GPCRs [33,59]. Microsecond cMD simulations have been able to capture beginning inactivation of the apo H 2 R and histamine-bound H 2 R in the absence the G s or G s -α5. GaMD simulations yield mostly similar results on this aspect, while with certain improved sampling of larger conformational space in the receptor residue distances plotted in Figures 6 and 7. Future GaMD simulations may be applied to investigate slower conformational changes underlying activation of the H 2 R and coupling of the receptor with different IBPs.
The data above raises the question about the role of agonists for the stabilization of the active state. A recent study of β 1 AR has shown that the presence of an agonist can shift the conformational equilibrium towards a pre-active state that facilitates G-protein binding [59]. In case of β 1 AR, the conversion between these states occurs on a millisecond to second timescale [59] and can therefore not fully be covered by MD simulations. Despite the limited timescale accessible to MD simulations, there is some evidence from our simulations that the presence of histamine affects the conformational equilibrium-namely the volume of the G protein binding pocket compared to the apo H 2 R (Figure 8).
In the presence of histamine, the volume of the pocket exhibits larger variations between the individual runs (Figure 8c), which might affect G protein binding. Reciprocally, G protein coupling was reported to increase agonist binding affinity in case of the adenosine A 2A receptor between 10-and 40-fold [60,61]. For H 2 R, competitive binding experiments also indicate a 20-fold change of K i in the presence of the G protein [62]. We have investigated this aspect in more detail by calculating the binding free energies for histamine in the presence and absence of G s from the cMD simulations (Table 1). In the presence of G s , histamine is bound approximately 1 kcal·mol −1 tighter, which reflects the same trend as observed in previous experiments. Interestingly, this effect is not observed, when only the G s -α5 is present instead of the intact G protein (Table 1). This finding is a first clue that G s -α5, although it stabilizes the cytosolic part of H 2 R (Figures 6-8), may be unable to enhance histamine binding affinity. Table 1. Binding energies between histamine and the H 2 receptor. Analysis were performed using cMD simulations of the binary histamine-H 2 R complex (H-cMD1/2), the ternary histamine-H 2 R-G s complex (H-G s -cMD1/2), and the ternary complex including only G s -α5 instead of G s (H-α5-cMD1/2). Averages and standard errors of mean (n = 10,000 frames from the second half of the simulations) were calculated for the binding energies using MM/GBSA.

System
Binding Taken together, the simulations of H 2 R in complex with different binding partners reveal that the presence of an intracellular partner (G s or G s -α5) appears crucial for maintaining an active conformation, whereas histamine alone is insufficient to stabilize this state. These mechanistic properties are in good agreement with those described previously for the related β 1 AR [59] and β 2 AR [58]. Our simulations also give first evidence for a weak allosteric coupling between the intraand extracellular binding site in H 2 R, which has also been described for other GPCRs [13,63,64]. However, more and longer MD simulations or experiments with atomic resolution (NMR, X-Ray, cryo-EM) will be needed to obtain a comprehensive picture of the underlying molecular mechanisms.

Molecular Modeling of the H 2 R-G s Complex
For model generation, the crystal structure of the human β 2 adrenoceptor in complex with a G s protein (PDB code 3SN6 [43]) was used as template. Both GPCRs exhibit a sequence identity of 31% (sequence similarity 66%). The homology model was generated with Modeller 9.16 [65] and comprises residues 15-304 of H 2 R. For the G s protein, the sequence from PDB entry 3SN6 was kept with only one minor modification: Residues 60-87 of the α-subunit, which are not resolved in the template crystal structure, were replaced by a heptameric GSGSGSG-linker in the modeled complex. The camelid antibody fragment present in the template crystal structure was also kept in the modeled complex. The resulting model exhibited a good sterochemistry with 97% of the residues in the most favored regions of the Ramachandran plot and no steric clashes >0.35 Å.

Overview of the Molecular Dynamics Simulations
Preparation of the H 2 R-G s complex including minimization, membrane embedding, and equilibration followed the protocol described for the H 1 R [66]. The parameters for histamine were also adapted from this work. An overview over all simulations performed is given in Table 2, which is subdivided in two sections: The first part contains information about the initial simulations that were used to derive the binding mode of histamine. The second part lists the simulations conducted for the final model to investigate H 2 R conformational stability.
Simulations to derive the binding mode of histamine were done with Gromacs [67]; the simulations conducted to investigate H 2 R conformational stability were done with AMBER [68]. The rationale for using two different simulation programs is that the PLUMED plugin required for metadynamics simulations works most efficiently with Gromacs, whereas the Gaussian-accelerated simulation method is not available for this program. The metadynamics protocol used for the determination of the histamine binding mode is described in detail in reference [66]. Briefly, Gromacs 2016.3 [67] was used in combination with the plumed 2.3.1 [69] plugin for all metadynamics simulations. The well-tempered metadynamics approach established previously by Saleh et al. [44], relies on the z component of the distances between the Cα atom of W247 6.48 and the histamine ammonium nitrogen as collective variable (CV). The lower and upper walls z low and z up were chosen as 0.3 nm and 4.9 nm as CV boundary conditions. To discourage an irrelevant exploration of the bulk solvent the same bell-shaped funnel as described in [66] was used. The initial metadynamics simulation to capture the histamine binding pathway were run with an initial bias height of 7 kJ·mol −1 , a gaussian width of 0.1 and bias factor 50. The multiple walker simulations were performed with bias factor 20 and an initial bias height of 5 kJ·mol −1 . Table 2. Overview of the simulations performed. The table lists all MD simulations conducted for the present study and the respective names by which they are referred to in the figures and manuscript text. It is stated whether a simulation was a conventional MD simulation (cMD), a Gaussian accelerated MD simulation (GaMD), or a metadynamics simulation. The number of runs and the respective simulation times are listed. Furthermore, the table describes the composition of the respective systems, i.e., the box dimensions and whether the histamine H 2 receptor (H 2 R), the ligand histamine (HSM), the G s protein ( ), or the G s -α5 helix ( ) was present. The symbol (×) denotes the absence of the respective component in the setup.

Investigation of H 2 R Conformational Stability
The investigation of H 2 R conformational stability in the presence of different intra-and extracellular binding partners ( Table 2) was performed with Amber [68] using the force fields ff14SB for the proteins, lipid14 for the DOPC molecules, and GAFF for histamine [70][71][72]. The initial structures were converted from the final structure of the cMD refined cluster5. Each system was energy minimized and equilibrated with the following protocol: Energy minimization consisted of three consecutive steps with restraints applied to different subsets of atoms (first to all atoms except for water molecules, second only to the Cα atoms, and third without any restraints). Each minimization stage was composed of 2500 steps with the steepest descent algorithm followed by 2500 steps of the conjugate gradient algorithm using a harmonic potential with a force constant of 10 kcal·mol −1 ·Å −2 for the atom restraints. The equilibration was also performed in three consecutive parts which were conducted with a time step of 2 fs: First, the system was heated from an initial temperature of 10 K to 310 K within 0.1 ns while all atoms except for water molecules were restrained with a force constant of 5 kcal·mol −1 ·Å −2 . Keeping pressure and temperature constant (NPT ensemble), another 0.4 ns equilibration was performed where only the Cα atoms were fixed. Finally, the whole system was equilibrated for 0.5 ns without any restraints. The time step was 2 fs for all systems. The temperature was kept constantly at 310 K by a Berendsen thermostat [73]. A reference z pressure of 1 bar and a reference surface tension of 1.1 nm·bar were applied using surface-tension coupling. The SHAKE algorithm [74] was used during equilibration and production runs for bonds with hydrogen atoms. Periodic boundary conditions were set for x, y and z direction. The same general simulation parameters and equilibration strategies were used for both cMD and GaMD simulations. The GaMD protocol was inspired by a similar study for the muscarinic M2 receptor published by Miao et al. [33]. For the statistical calculation of the acceleration parameters, each GaMD run was preceded by a short 10.4 ns cMD simulation. This step was followed by a 32 ns equilibration in which the calculated boost potentials were added. For all GaMD simulations, the dual-boost mode was chosen which means that both the total potential energy and the dihedral energy terms were boosted. The reference energy was set to the lower bound (E = V max ). Averages and standard deviations of the potential energies were updated every 400,000 steps (800 ps). For the standard deviations of the boost potentials applied to the total potential energy and the dihedral energy, an upper boundary of σ 0P = σ 0D = 6.0 kcal·mol −1 was defined. Table 3 gives an overview of the boost potentials (averages and standard deviations) that were actually added for the different simulation runs using these settings. Table 3. Boost potentials of the Gaussian accelerated molecular dynamics (GaMD) simulations. All GaMD simulations were run with the dual boost option which means that both the total potential energy and the dihedral energy were boosted. Averages ± standard deviations of the potentials applied are listed for the different simulations.

Structural Analysis
Post-processing and subsequent analysis were mainly performed using cpptraj from AmberTools 18 [68]. For the assessment of contacts, the nativecontacts command was used with a specified distance criteria of 5 Å. Interaction energies between ligand and receptor were calculated based on the MM/GBSA method using the script mm_pbsa.pl with default parameters [75][76][77].
The volume of the G protein binding pocket was calculated utilizing the tool POVME2 [78] according to a strategy published by Li et al. for the adenosine A 2A receptor [79]. The center of the binding pocket was defined as the geometric center of the last five Cα atoms of the α5 helix of G s α. In case of simulations without the G s protein, an overlay with a G s -bound H 2 receptor was performed to determine the respective coordinates. A 12 × 12 × 12 Å 3 rectangular box was defined as the maximum extension of the binding pocket. PISA [80] was used for interface analysis. Structure visualization was done with UCSF Chimera [81] and VMD [82] while plots were created with gnuplot [83]. For the visualization of sequence alignments, the T E Xshade package [84] was used.

Conclusions
The H 2 R represents an established drug target for the treatment of gastric diseases [8,9] and is also actively studied as a target for reducing the risk of heart failure [11]. Since there is yet no experimental structure available for H 2 R, we have generated a computational model of the ternary histamine-H 2 R-G s complex and assessed its conformational dynamics. Ligand binding and conformational changes of GPCRs mostly occur on timescales longer than microseconds [13,14] and can be covered only in part by cMD simulations. Therefore, metadynamics and GaMD simulations have been applied to enhance conformational sampling.
A multiple walker metadynamics simulation was used to identify the histamine binding mode, which was subsequently refined by cMD simulations according to an existing protocol [41]. The resulting histamine binding mode was in good agreement with data from site-directed mutagenesis [46] underscoring the usefulness of a cMD/metadynamics combination for the determination of GPCR ligand binding modes.
The analysis of the H 2 R-G s interface revealed a comprehensive network of polar interactions supporting the H 2 R preference for G s as a major signaling partner. However, there are additional structural factors that affect G protein coupling specificity, e.g. the exact position of helix H6 [50][51][52] or the existence of structural intermediates on the binding pathway [53]. Therefore, future studies of H 2 R in complex with different G proteins or investigations of different G protein binding modes will be required to assess the role of these structural factors for H 2 R coupling preferences.
We also compared dynamics of the histamine-H 2 R-G s ternary complex, the histamine-H 2 R binary complex and the apo H 2 R to assess the role of different interaction partners in the conformational stability of the active H 2 R. In addition, a truncated ternary complex that contains only the C-terminal G s -α5 helix instead of the entire G s protein was simulated to investigate whether the α5 helix could sufficiently stabilize the active state. To enhance conformational sampling and facilitate large structural rearrangements, 1-µs GaMD simulations were performed in addition to cMD simulations.
During the simulations of the apo H 2 R and the H 2 R-histamine complex, several structural features typical for inactive GPCRs emerged: An inward motion of the intracellular end of helix H6, the formation of the "ionic lock", and a contraction of the G protein binding pocket. In contrast, these motions were not observed in the presence of G s or G s -α5, suggesting that the presence of an intracellular partner is crucial for maintaining an active GPCR conformation, whereas histamine alone is insufficient to stabilize this state. This was in line with previous studies showing that the G protein is the key determinant for generating an active conformation, whereas an agonist alone is insufficient to stabilize the active state [58,59]. In addition, the simulations also provided first indications for a weak allosteric coupling between the intra-and extracellular binding sites, which has also been described for other GPCRs [13,63,64]. However, our investigations also showed that certain motions in H 2 R are so slow that they cannot be fully sampled even using 1-µs GaMD simulations. Therefore, future cMD and GaMD simulations with longer simulation times or modified simulation parameters will be required to investigate slower conformational changes underlying H 2 R activation.