Molecular Basis of Class B GPCR Selectivity for the Neuropeptides PACAP and VIP

Introduction

Pituitary adenylate cyclase-activating peptide (PACAP, ADCYAP1, P18509) and vasoactive intestinal polypeptide (VIP, VIP, P01282) are two important neuropeptides for neural development, body calcium homeostasis, glucose metabolism, circadian rhythm, thermoregulation, inflammation, feeding behavior, pain modulation, stress and related endocrine response (Harmar et al., 2012; Bortolato et al., 2014; Culhane et al., 2015; Liao et al., 2019b; Liao et al., 2019c). Their important roles are mediated by activating three class B G protein-coupled receptors (GPCRs) including the PAC1 (ADCYAP1R1 and related splice variants), VPAC1 (VIPR1), and VPAC2 (VIPR2) receptors, which share over 60% sequence similarity within this GPCR subtype. Interestingly, PACAP and VIP display distinct selectivity for the PAC1 and VPAC1/2 receptors (Vaudry et al., 2009): PACAP (PACAP38 and the C-terminally truncated PACAP27) and VIP have comparable high affinity for the VPAC1 and VPAC2 receptors (Gottschall et al., 1990; Lam et al., 1990), whereas the PACAP peptides are 400–1,000 fold more potent than VIP as agonists for PAC1 receptor (Gottschall et al., 1990; Dautzenberg et al., 1999).

From recent progress in X-ray crystallography and cryogenic electron microscopy (cryo-EM) (Hollenstein et al., 2013; Yang et al., 2015; Jazayeri et al., 2017; Zhang et al., 2017a; Zhang et al., 2017b; Zhang et al., 2018; Duan et al., 2020; Kobayashi et al., 2020; Liang et al., 2020; Wang et al., 2020), several class B receptor structures have been established. We have now integrated the structural information with simulations, theory and available experimental findings to examine: 1) differences in intrinsic PAC1 and VPAC1/2 receptor dynamics; 2) differential PACAP and VIP interactions and energetic impacts at the receptor subtypes; and 3) the interactions and dynamical features of PACAP-induced activation of the PAC1null receptor. Class B GPCRs have been suggested to have receptor features and dynamics distinct from those described for class A receptors, presumably because of the presence of the extracellular domain (ECD) (Liao et al., 2017), the apparent absence of “toggle switches” (Li et al., 2013), and other signature features. The class B receptors are likely to adopt unique conformational changes in the heptahelical transmembrane (7TM) domain following neuropeptide binding for receptor activation. Different from some class B receptors (Hollenstein et al., 2014; Yang et al., 2015; Graaf et al., 2016; Jazayeri et al., 2016; Jazayeri et al., 2017; Song et al., 2017; Zhang et al., 2017a; Zhang et al., 2017b; Liang et al., 2018), the detailed and comparative examination of the PAC1 and VPAC1/2 receptor function and actions have not been well investigated. Accordingly, we have integrated long-timescale molecular dynamics (MD) simulations with other computational and theoretical techniques to study the conformational ensembles and transitions of two PAC1 receptor variants (namely PAC1null and PAC1s) (Liao et al., 2019c), VPAC1, and VPAC2 receptors, as well as their interactions with the neuropeptides PACAP and VIP. The PAC1null receptor has no intracellular loop 3 (ICL3) inserts (neither hip nor hop inserts) from alternative splicing, and represents one of the most prevalent receptor isoforms in the central nervous system (CNS). The PAC1s receptor is similar to the PAC1null variants, but has a 21-amino acid deletion (residues 89–109) in the ECD. For reasons unknown, the PAC1s receptor is not highly expressed in the CNS (or any tissue) and yet this variant was chosen for structural determination (Sun et al., 2007; Kumar et al., 2011; Wang et al., 2020). Distinct from structural determination, our simulations and analyses detail the receptor conformational transitions and peptide interactions that are biologically relevant to the receptor subtype. While the PAC1 receptor has been recognized as an emerging target for stress-related disorders (Ressler et al., 2011), our studies, for the first time, reveal the molecular basis underlying PACAP receptor selectivity which may be essential for the rational design of effective therapeutic agents.

Models and Methods

Model Preparation and MD Simulation Setup

The ligand-free and ligand-binding full-length receptor systems. Our homology modeling was carried out with the protein modeling program Prime (Schrödinger, Inc.) (Zhao et al., 2011; Li et al., 2011) which built the full-length VPAC1 and VPAC2 including the ECD, the 7TM, as well as the linker loop (see details in Supplementary Table S2). The PAC1s model was modified from our previous PAC1null receptor model (Liao et al., 2017) by deleting the 21-amino acid in ECD. For each receptor, we generated a number of initial models which were distinct in the ECD orientation to better sample different conformational states (Supplementary Figure S1). The Membrane Builder in CHARMM-GUI (Jo et al., 2008) was employed to build the protein/membrane complex systems, each of which contains a receptor with/without a ligand, a lipid bilayer of ∼ 225 POPC molecules, ∼ 27,000 TIP3P water molecules, counter ions, and 0.12 M NaCl, totaling ∼ 120,000 atoms in a periodic box ∼94 × 94 × 147ų.

The ligand-ECD complex systems. To understand the distinct binding affinity of PACAP/VIP for PAC1null, PAC1s, VPAC1 and VPAC2 receptors, we simulated the ECD of each receptor and PACAP/VIP in a solvent environment. PACAP/VIP was initially placed near the position based on known peptide-binding modes to class B GCPRs (Sun et al., 2007; Grace et al., 2010; Parthier et al., 2007; Pioszak and Xu, 2008) (Supplementary Figure S2). We used the PDB Reader of CHARMM-GUI (Jo et al., 2008) to prepare protein then solvated and neutralized with addition of 0.12 M NaCl in VMD (Humphrey et al., 1996) (see our simulation summary in Supplementary Table S1).

Our simulations were performed with the CHARMM36-cmap force field (Best et al., 2012). Each system went through energy minimization, 285-ns multiple-step equilibrations, and MD simulations using the NAMD package (Phillips et al., 2005). Both equilibration and production runs were performed in the NPT ensemble (310K, 1 bar, Langevin thermostat and Nose-Hoover Langevin barostat) with a time step of 2 fs. All of the bond lengths to the hydrogen atoms in the protein or lipid molecules were constrained. The particle mesh Ewald (PME) technique was used for the electrostatic calculations. The van der Waals and short-range electrostatics were cut off at 12.0 Å with switch at 10.0 Å. Each system was carried out 1000 ns for 4–5 replicas.

Enhanced sampling for ligand-binding simulations. Adaptive tempering was employed to enhance conformational sampling under PACAP or PACAP6-38 insertion for ∼ 400 ns, in which the simulation temperature was dynamically updated as a continuously random variable in the range of 310–360 K with the Langevin equation (Zhang and Ma, 2010) implemented in NAMD (Supplementary Figure S3).

Markov state model

Markov state models (MSMs) of molecular kinetics were applied on the collection of MD trajectories, from which we identified the transition pathways and estimated the long-time statistical dynamics between the conformational states of interest (Pande et al., 2010; Prinz et al., 2011). The ligand-free full-length receptor MD trajectories were used for MSM construction. We used the MSMBuilder 3.8.0 program Bowman et al., 2009a; Bowman et al., 2009b to construct the transition matrices of MSM. Trajectories that contain stable conformational states (measured by C_α root-mean-square deviation (RMSD) in Supplementary Figure SM2) were prepared by saving the C_α coordinates of a GPCR protein with atomic indices. Thus, we have 4,519–14,281 conformations in each trajectory. We grouped the conformational ensembles into a set of clusters, called microstates, based on C_α RMSD using the k-centers algorithm with 33–55 cluster centers Bowman et al., 2009b.

With the microstate discretization, we used the maximum likelihood estimation to build the reversible transition matrices at the lag time series. The evolution of the transition probability between closed and open conformational states of a receptor protein was evaluated at a series of lag times τ, 2τ, …, nτ by performing the Chapman-Kolmogorov test (Prinz et al., 2011) which confirmed that our model at the chosen lag time is Markovian (Chodera et al., 2007; Noé et al., 2009; Prinz et al., 2011) (Supplementary Figures SM4,SM5). Finally, this was achieved at a lag time of 1.2–2.4 ns as judged by flat lines in the implied timescales plots (Supplementary Figure SM3) (see Ref. Liao et al., 2019a for our detailed implementation of the Chapman-Kolmogorov test).

Using the transition-path theory (TPT) (Weinan and Vanden-Eijnden, 2006; Berezhkovskii et al., 2009; Metzner et al., 2009; Noé et al., 2009), we calculated the transition path connecting between the stable conformational states in a transition matrix and estimated the time quantity to travel from one set of states to the others by the lag time divided by the minimum net flux in the transition path (Supplementary Table S3). Thus, the transition time refers to the time quantity to complete one transition between two stable conformational states. The division of macrostates were calculated from the eigenfunction structure using the Perron Cluster Cluster Analysis (PCCA) method (Noé et al., 2007). Supplementary Figure SM6 illustrates the 3D projections of the microstates and macrostate divisions in accordance to ECD orientation and position by backbone center-of-mass (COM) distances between ECD and extracellular loop 1 (ECL1) (d_ECD-ECL1) and extracellular loop 3 (ECL3) (d_ECD-ECL3), respectively. Detailed constructions and validation of MSM were described in the Supplementary Material.

Data Analysis

Binding free energy (ΔG_b) of PACAP/VIP peptide bound to the ECD and per-residue free energy decomposition were computed over the last 500 ns MD trajectories using the Molecular Mechanics Generalized Born Surface Area (MM-GBSA) method (Miller et al., 2012) by MMPBSA.py tool (version: 14.0) in the Amber package (Case, 2018). The second modified Bondi radii set were used in the GB method. ΔG_b was given as an average of five replicas in each system. Conformational analysis and distance analysis were performed with TCL scripts implemented in VMD 1.9.1 (Humphrey et al., 1996) and plotted by matplotlib. The average distances involving a group of residues were calculated by the distance between the backbone COMs of groups of residues. Water density map were calculated using MDAnalysis python package. The tilt angle of ECD is defined as the angle between the vector along the N-terminal helix and the Z-axis. For example, the vector along the N-terminal helix of PAC1 ECD is calculated by summing the C-O vectors along the helical residues 30–47.

Results and Discussion

Conformation Ensembles and Transitions of PAC1 and VPAC1/2 Receptors

Building on our previous study of PAC1null receptor (Liao et al., 2017), we first examined the conformations of full-length PAC1null, PAC1s, VPAC1, and VPAC2 receptors, as well as the transitions within the conformational ensembles, which are key in understanding the differential actions of these receptors beyond amino acid sequences. With MD simulations totaling 20 microseconds to sample different conformations, we constructed MSM transition pathways between conformational microstates for PAC1s and VPAC1/2 receptors to compare with those shown for the PAC1null receptor described previously (Liao et al., 2017) (Figure 1). Similar to the PAC1null receptor, the ligand-free PAC1s and VPAC1/2 receptors display a wide diversity of conformational states with respect to the ECD orientation (θ) and ECD-7TM backbone COM distance (d), which reflect in part the actions of the linker region that tethers the ECD to the 7TM. Free energy maps based on the conformation density are shown in Supplementary Figure S4, where the low energy regions correspond to the stable conformational states in Figure 1 are circled with the same color. The energy barriers between the stable conformational states vary from 0.5 to 3.5 kcal/mol. We define the closed (θ > 80 degree, d < 55 Å), semi-open (θ = ∼60 degree, d < 55 Å), and open states (θ = ∼40 degree, d < 55 Å) to distinguish these conformational states. ECD-unrestricted conformations refer to those with ECD not touching the 7TM (d > 55 Å), which include most intermediate states (colored in yellow in Figure 1), as well as some stable conformational states in systems of VPAC1 and VPAC2, e.g., the state in red in Figure 1C and the state in green in Figure 1D. Without binding to 7TM, they are capable to change to other stable states much faster. The conformational states in red and orange in Figure 1D that have the ECD-ECL1 backbone COM distance ∼10 Å shorter than ECD-ECL3 backbone COM distance are labeled as the on-side conformations.

FIGURE 1

FIGURE 1. (A) Transition pathways for (A) PAC1s, (B) PAC1null, (C) VPAC1 and (D) VPAC2 receptors. Stable conformational states are colored in blue, red, green, and orange respectively. Intermediate states are colored in yellow. Most intermediate states are ECD-unrestricted conformations. The VPAC1 conformation in red and the VPAC2 conformation in green also belong to the ECD-unrestricted conformational ensemble. PAC1null is remade from our recent study (Liao et al., 2017). The kinetic transition time between certain conformational states are labeled with units in nanoseconds.

For the PAC1null receptor, we previously showed that the 21-amino acid loop in ECD interacts with ECL3 to potentially sustain the open time of the receptor (Liao et al., 2017). Consistently in this work, with the deletion of the 21 amino acids in PAC1s, the receptor microstate transitions between the open conformational state and closed conformational states (Figure 1A) are faster by a factor of 2 (transition time up to 400 µs for PAC1s and up to 1,000 µs for PAC1null, Supplementary Table S3). The free energy map also exhibits ∼3 kcal/mol higher energy barrier between the open and closed states in PAC1null compared with it in PAC1s (Supplementary Figure S4). While the functional implications of these differences remain to be investigated, the open state of PAC1null, once formed, may be sustained longer than that for PAC1s to facilitate peptide trapping and PAC1null receptor activation. Overall, the conformational microstate transitions for the VPAC2 receptor appear the fastest (0.8–80 µs) compared to the VPAC1 (40–600 µs) or the PAC1 receptor variants (21–366 µs), in agreement with the free energy difference between states in Supplementary Figure S4. Therefore, despite significant primary amino acid sequence homology, the different PAC1 and VPAC receptors display unique intrinsic ensemble dynamics and conformational microstate transition times.

Several of the ligand-free receptor conformations resemble the peptide-bound receptor states identified in the structural solutions with respect to ECD orientation (θ) and ECD-7TM (d) or ECD-ECL backbone COM distances. For example, both the PAC1null and PAC1s receptors could adopt stable open conformations (θ = ∼40°, 45 < d < 55Å), which have been seen in the peptide-bound PAC1 (PDBID: 6M1I) and other Class B receptor structures, including the glucagon receptor (GCGR, PDBID: 5YQZ), parathyroid hormone 1 receptor (PTH1, PDBID: 6FJ3), and glucagon-like peptide-1 receptor (GLP-1, PDBID: 5VAI) with minimum RMSD of 4.2Å (Supplementary Figure S5). In addition, the “semi-open” conformations of the VPAC1 receptor are similar to the CGRP-bound calcitonin receptor-like receptor (CLR) structure (PDBID: 6E3Y, θ = ∼60°, d < 55Å) with RMSD of 5.4 Å. Moreover, the VPAC2 receptor “on-side” conformations with ECD in contact with ECL1 share similar features (θ < 80°, short ECD-ECL1 distance) with an antagonist-bound GCGR structure (PDBID: 5XEZ) with minimum RMSD of 4.6 Å. Beyond the relative orientation and position of ECD to 7TM, the 7TM conformations in the ligand-free receptor approximated those for the inactive GPCR structures (with lowest RMSD = 2.0–2.5Å in 7TM) than the peptide-bound active receptors (Supplementary Figures S6,S7), consistent with the actions of peptides to initiate transmembrane dynamics for receptor activation. In aggregate, these findings suggest that the PAC1 and VPAC1/2 receptors, even without peptide binding, have dynamic transitions between the open and closed conformational states, which prepare them for ligand binding and subsequent receptor activation/deactivation.

Distinct Affinity of PACAP and VIP to the Receptor ECDs

For Class B receptors, the ECD acts as a high-affinity peptide trap and subsequent peptide-ECD dynamics result in the presentation of the peptide N-terminus to the orthosteric binding site in the 7TM to initiate receptor activation. Our MSM analyses suggest the dynamic processes (peptide binding, large-scale ECD rotation, and peptide insertion) occur on a time scale greater than tens of microseconds. Thus, we applied a total of 40 microseconds MD simulations on the ligand-ECD systems; with the ECD-peptide complex trajectories we employed the MM-GBSA method (Miller et al., 2012) to estimate the binding free energy (ΔG_b) of PACAP or VIP to the ECD of PAC1null, PAC1s, VPAC1 and VPAC2 receptors, as well as the ΔG_b contribution per residue. Our results (Figure 2A) show that PACAP binding to PAC1null receptor is the most favorable among receptor and peptide combinations. Specifically, ΔG_b of PACAP is around 1.4 times of ΔG_b of VIP in binding PAC1null ECD; the difference decreases in PAC1s and VPAC1, and become comparable for PACAP and VIP binding VPAC2 ECD (Figure 2A). This provides an energetic basis for PAC1 receptor binding selectivity for PACAP (K_d ≈ 0.5 nM) over VIP (K_d >500 nM) (Vaudry et al., 2009), which corresponds to that PACAP-binding free energy (−13.2 kcal/mol with K_d converted by ΔG_b = RTlnK_d, T = 310 K) is about 1.5 times of VIP (−8.9 kcal/mol) in binding PAC1null, while VPAC1 and VPAC2 receptors appear to exhibit comparable binding affinity (K_d ≈1 nM) for PACAP and VIP (Vaudry et al., 2009). Notably, the experimental binding energies were obtained with full-length receptor systems, and the interactions between the ligand and 7TM may also affect the magnitude of the ΔG_b values for comparison with our calculation estimations. Therefore, here we compare the relative values of ΔG_b for the different systems as opposed to the absolute magnitude which differs from prior experiments.

FIGURE 2

FIGURE 2. Binding selectivity and recognition of PACAP and VIP on PAC1null, PAC1s, VPAC1 and VPAC2 receptors. (A) Binding free energy (ΔG_b) of PACAP and VIP on PAC1null, PAC1s, VPAC1 and VPAC2 receptors respectively by averaging the ΔG_b values from five replicas with bars representing the standard deviation. Representative binding conformations are displayed on the right, where residues to calculate free energy decomposition in (B) are colored in blue; only residues 1–27 are displayed for PACAP and VIP. (B) ECD per-residue free energy decomposition (average of five replicas) in PAC1 and VPAC1/2 binding PACAP (red column) and VIP (blue column). β3–β4 including the loop between is labeled. The 21 amino acids (residues 89–109) are labeled in PAC1null. (C) Close view of the N-terminal residues 1–13 of PACAP interacting with transmembrane orthosteric pocket of PAC1 receptor. Hydrogen bonds are shown in dashed lines. Comparison of PACAP with VIP at residues 1–13 by direct mutations is shown at the bottom right.

To examine the ligand-ECD complex conformations, we performed structural clustering on the last 200- ns trajectories and identified six complex conformational states based on ECD orientation and the bound peptide in a vertical pose for receptor activation (Figure 2A; a-e). State (a) is found in PACAP binding PAC1null, PAC1s and VPAC2 receptors, which was resolved in the PAC1 receptor cyro-EM structure (PDBID: 6M1I, Figure 2 state(f)). State (b) in which the ECD tile angle is decreased compared to state (a) is found prevalently in PACAP/VIP binding the VPAC1 receptor. While the PAC1s receptor adopts state (a) upon PACAP binding, it establishes state (c) if bound to VIP. Similarly, the VPAC2 receptor adopts state (a) in PACAP binding, but it alters to state (e) with VIP, which is comparable to the bound ligand bound CRF1 conformation (PDBID: 2L27).

According to previous theoretical and experimental findings (Vaudry et al., 2009; Wang et al., 2020), we examined the interacting residues between the receptor ECD (Figure 2B) and the C-terminal region of PACAP/VIP (amino acid residues 13–26). The PACAP and VIP C-terminal hydrophobic V¹⁹xxY²²LxxV/I²⁶ sequence appears to contribute to the hydrophobic interactions with PAC1null/PAC1s and VPAC1/2 receptor residues I61^PAC/L70^VPAC1/I59^VPAC2; the PACAP/VIP R¹⁴K¹⁵ residues contribute to hydrogen bonding with the ECD (Supplementary Figure S10). The β3-β4 loop in both PAC1null/PAC1s and VPAC2 receptor (environs L⁸⁰FxxF^84,PAC and V⁷⁸FxxF^82,VPAC2 respectively) interacts with PACAP/VIP (Figure 2, state (a)). In contrast, the ECD helix 1 rather than the β3-β4 loop in VPAC1 receptor interacts with the peptide (Figure 2, state (b)). In the different receptor ECD C-terminal region, segments F¹³¹xxxxF^136,PAC1null, F¹¹⁰xxxxF^115,PAC1s, Y¹¹⁸xxxxxL^124,VPAC1, and V¹⁰⁶xxxxY^111,VPAC2 also contribute to ligand-receptor hydrophobic interactions.

For the PAC1null receptor, there are a few notable differences in PACAP and VIP binding. Firstly, both states (a) and (d) in Figure 2 are identified in PAC1null-PACAP binding; state (a) is similar to the cyro-EM structure (Wang et al., 2020) and state (d) is close to the Nuclear Magnetic Resonance (NMR) model (Sun et al., 2007). PAC1null receptor VIP binding can also be in state (a). Secondly, PAC1null receptor hydrophobic residues 72–140^PAC1null play major role in PACAP binding, such as M⁷², L⁸⁰FxIF⁸⁴xP⁸⁶, VW⁹⁰, I⁹⁵, L¹⁰⁶xLxxM¹¹¹, F¹²⁷, F¹³¹, F¹³⁶, and Y¹³⁹ (Figure 2B); only E¹⁴² and D¹⁴⁵ at the ECD C-terminal region contributes significantly to hydrogen bonding with the ligand. Different from PACAP, VIP can also interact with charged/polar residues 78–87^PAC1null such as E⁷⁹, R⁸², N⁸⁵ and D⁸⁷. Lastly, PACAP interacts with residues 92–130^PAC1null extensively, while VIP interacts with few. This also supports the finding that another short variant with deletion of residues 53–109^PAC1null significantly reduces VIP binding affinity over PACAP (Dautzenberg et al., 1999). These remarkable interaction differences together contribute to the binding preference for PACAP. Compared with PAC1null, PAC1s missing the 21-amino acid loop loses the interaction formed by residues 86–111 of PAC1null which results in decrease in ligand-ECD binding affinity.

Integrating our finding on ligand-ECD interactions, PAC1 receptor cyro-EM structure (Wang et al., 2020) and biochemical data, the C-terminal region of PACAP/VIP plays a major role in binding with ECD, while the peptide N-terminus interacts with the receptor ECLs and 7TM activation site. Based on alignment with the cyro-EM structure (Wang et al., 2020), we simulated PACAP_1-13 insertion in PAC1 7TM, without PACAP C-terminal residues 13–26 associations to the ECD. Within a few hundred nanoseconds, PACAP residues 1–4 rearranges to form a D3^PACAP-R199^PAC1 interacting pair (Figure 2C). The depth of peptide N-terminus into the 7TM is not significantly changed, but H1^PACAP reorients to form hydrogen bond with ECL3, while PACAP residues D8, Y10, R12 and Y13 interact with ECLs (Figure 2C). These interactions result in outward conformational transition of TM6 at the intracellular face of receptor. However, in longer simulations, PACAP1-13 is also observed to shift away from its primary position toward the V-gap between TM5 and TM6, allowing TM6 to transition back to former inactive conformational state. Thus, these observations show the importance of PACAP C-terminus-ECD binding to constraint the peptide N-terminus within the 7TM orthosteric interaction site and sustain the activation process.

PACAP-Induced Activation of PAC1null Receptor

To further understand how PACAP interacts with the full-length receptor, we simulated the PAC1null receptor (a homology model of the inactive state (Siu et al., 2013; Liao et al., 2017)) bound to the PACAP_1-38 peptide, with the first 13 residues partially inserted into the receptor 7TM. We use the Wootten numbering (Wootten et al., 2013) in this section to better demonstrate the residue positions in PAC1null receptor. Given the timescale of PAC1R activation, we combined microsecond scale MD simulation and adaptive tempering (Zhang and Ma, 2010) to enhance conformational sampling under these conditions. The N-terminus of PACAP underwent conformational rearrangement in the orthosteric site with a clear TM6 outward movement of ∼ 4 Å with the TM6-TM3 distance close to the relaxed GCGR and GLP-1 receptor structures (Supplementary Figures S11,S12), indicating an “active” state.

A detailed examination suggests that the formation of this state results from several key interactions. As shown in Figure 3B (top panel), D3^PACAP and PAC1null R199^2.60 near the receptor extracellular face appears to migrate closer together in the first 700 ns, but becomes entrapped in a separate conformation for another 1500 ns. During adaptive tempering, the average distance between H1^PACAP and the surrounding receptor K206/D207^ECL1, E374/E380/E385^ECL3, and K154^1.40 decreases gradually and appears to form hydrogen bond (H-bond) networks with an average small distance of ∼5 Å. The hydrogen bonding between the D3^PACAP sidechain and H1^PACAP amino terminus weakens during adaptive tempering resulting in the ability for D3^PACAP to form a new pairing with receptor R199^2.60. As a result of D3^PACAP-R199^2.60 interactions, the average distance between R199^2.60 and N240^3.43, H365^6.52 and Q392^7.49 increases significantly (Figure 3B, middle panel), indicating weakened contacts between TM2, TM3, TM6, and TM7 (Figure 4D). At the receptor intracellular face, K180^ICL1 transitions from a solvent orientation to face E247/R185, and decreasing K180^ICL1-E247^3.50/R185^ICL1 distance (Figure 3B, bottom panel) converts the PAC1null receptor from an apparent “inactive” to “active” conformation that allows G protein interactions. Direct mutagenesis on residues 4–13^PACAP (i.e., N-terminal G4A, I5V, S9N, S11T, Y13L found in VIP sequence) in the “active” conformation caused the revocation of the spreading H-bond network around H1 and ECLs and led to inward movement of TM6 (Supplementary Figure S12) and stronger TM2, TM3, TM6 and TM7 interactions (Figure 4D).

FIGURE 3

FIGURE 3. Representative conformation of PACAP-activated PAC1null receptor. (A) Structural alignment of PACAP-activated (red) and ligand-free (cyan) conformations. Ligand-free conformation is from our recent study(Liao et al., 2017) with only the 7TM shown. PACAP is shown in green cartoon and ECD is shown in surface representation. (B) Time evolution of average distance or pair distance of key charged/polar residues locating at the extracellular side, middle transmembrane, and intracellular side of the PAC1null receptor, which are displayed in (C) accordingly. Adaptive tempering period is in light orange strip background.

FIGURE 4

FIGURE 4. Time evolution of average distance or pair distance of key charged/polar residues locating at extracellular side, middle transmembrane, and intracellular side of PAC1null binding PACAP6-38, which are illustrated in (B) accordingly. Adaptive tempering period is in light orange strip background. (B) Final snapshot of PACAP6-38 bound PAC1null (cyan) superposed on PACAP-activated PAC1null structure (red), in which an inward shifting of TM6 was observed. (C) Wheel plot of major polar and hydrophobic interactions in the N-terminus of PACAP bound PAC1null (in red), N-terminus of PACAP6-38 bound PAC1null (in blue), and residue interactions within ligand-free PAC1null (in green). (D) Stack histogram of number of contacts between the intracellular half of TM6 (residues 344–360) and TM2, TM3, and TM7 in PAC1null receptor with: PACAP, PACAP deleted (ΔPACAP), PACAP6-38 (P6-38), and PACAP mutant (G4A, I5V, S9N, S11T, Y13L). Collections of the last 120 ns of each replica were used for the contact calculations for the PACAP, ΔPACAP, P6-38, and PACAP mutant systems. G1-G4 represent the four major ligand-free conformational states in Figure 1B. We used the last 960 ns of each 2- μs ligand-free system (Liao et al., 2017) for the contact calculations.

To investigate the effects of the first few residues at the N-terminal, we removed PACAP residues 1-5 from the PACAP-bound PAC1null complex, to simulate the effects of the PAC1 receptor antagonist PACAP6-38. Starting from an active conformation, we employed 45 ns adaptive tempering to accelerate the process, as conventional MD simulations show that longer times may be needed to relax an “active” conformational state (Supplementary Figure 12B). As shown in Figures 4A,B a 4∼6 Å inward movement of TM6 was observed, indicating the closing of the intracellular G protein binding site into an inactive state. Compared to the “active” conformation, E374^ECL3 is released and shifted to face the solvent. A return of enhanced TM2, TM3, TM6 and TM7 interactions are subsequently observed (Figure 4D), resulting in decreased average distance between R199^2.60 and N240^3.43, H365^6.52 and Q392^7.49, and the distance between the R350^ICL3-E406^H8 pair (Figure 4A).

Key interactions between the N-terminus of PACAP or PACAP6-38 and PAC1null receptor is summarized in Figure 4C comparing with residue interactions in the ligand-free conformations from our previous study (Liao et al., 2017). Compared with the recent PAC1 receptor cryo-EM structure (Wang et al., 2020), key contact like D3^PACAP-R199^2.60 was also found in our simulations, although the N-terminus of PACAP underwent conformational rearrangements in our simulations. Together, the C-term^ligand-ECD and N-term^ligand-7TM interactions present the unique structure features of class B GPCR in activation: ligand-ECD binding initially constraints the peptide close to 7TM; under optimized binding state, interactions between residues 8–15 of PACAP/VIP and ECLs help dock the N-terminus of the peptide into the pocket and initiate activation (Schäfer et al., 1999).

Conclusion

In summary, we have used advanced simulation and modeling techniques combined with prior experimental knowledge to examine the differential binding selectivity, conformation ensembles and microstate dynamic transitions among the PACAP/VIP receptor variants and subtypes, and a few key interactions in the orthosteric-binding pocket. The PAC1null, PAC1s, VPAC1 and VPAC2 receptors demonstrate unique conformation and transition states. Expectedly, from residue interactions and free energy calculations, PACAP C-terminal binding to PAC1null receptor ECD is favored over VIP. The PAC1s receptor variant harboring a 21-amino acid ECD deletion has been enigmatic because after identification from cloning, the variant appears absent or expressed at very low levels in all tissues examined to date; the PAC1null receptor by contrast is dominant. Our analyses show that the ECD 21-amino acid loop in the PAC1null receptor allows several functional advantages including potential enhancement of the receptor open state to facilitate peptide entrapment and increasing peptide binding affinity. The simulations reveal the receptor conformational changes upon ECD-bound peptide dynamics and peptide N-terminal presentation to the 7TM orthosteric activation site. Notably, the interaction of PACAP residues 8–15 to the ECLs is observed to help dock the peptide N-terminal into the orthosteric pocket and facilitate D3^PACAP–R199^PAC1null interactions that result in 7TM transitions and TM6 dynamics to initiate activation. Interestingly, the ECD not only functions for high affinity binding but also appears to constraint the orientation of the peptide N-terminus in the orthosteric site to maintain the activation state. These observations may be common to other Class B GPCRs and provide mechanistic insights for the development of therapeutics.

Data Availability Statement

The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

Author Contributions

CL and JR performed simulations and data analysis. All authors discussed and approved the contents of the manuscript and contributed to its final version by reading and editing.

Funding

The work was supported by grants R01-GM129431 (JL) from the National Institutes of Health.

Conflict of Interest

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

Acknowledgments

Computational resources were provided by VACC, PSC-ANTON, and XSEDE.

Supplementary Material

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

References

Berezhkovskii, A., Hummer, G., and Szabo, A. (2009). Reactive flux and folding pathways in network models of coarse-grained protein dynamics. J. Chem. Phys. 130 (20), 205102. doi:10.1063/1.3139063

PubMed Abstract | CrossRef Full Text | Google Scholar

Best, R. B., Zhu, X., Shim, J., Lopes, P. E., Mittal, J., Feig, M., et al. (2012). Optimization of the additive CHARMM all-atom protein force field targeting improved sampling of the backbone φ, ψ and side-chain χ(1) and χ(2) dihedral angles. J. Chem. Theor. Comput. 8 (9), 3257–3273. doi:10.1021/ct300400x

CrossRef Full Text | Google Scholar

Bortolato, A., Doré, A. S., Hollenstein, K., Tehan, B. G., Mason, J. S., and Marshall, F. H. (2014). Structure of class B GPCRs: new horizons for drug discovery. Br. J. Pharmacol. 171 (13), 3132–3145. doi:10.1111/bph.12689

PubMed Abstract | CrossRef Full Text | Google Scholar

Bowman, G. R., Beauchamp, K. A., Boxer, G., and Pande, V. S. (2009a). Progress and challenges in the automated construction of Markov state models for full protein systems. J. Chem. Phys. 131 (12), 124101. doi:10.1063/1.3216567

PubMed Abstract | CrossRef Full Text | Google Scholar

Bowman, G. R., Huang, X., and Pande, V. S. (2009b). Using generalized ensemble simulations and markov state models to identify conformational states. Methods 49 (2), 197–201. doi:10.1016/j.ymeth.2009.04.013

CrossRef Full Text | Google Scholar

Case, D. A. (2018) Amber 2018. (San Francisco: University of California).

Chodera, J. D., Singhal, N., Pande, V. S., Dill, K. A., and Swope, W. C. (2007). Automatic discovery of metastable states for the construction of Markov models of macromolecular conformational dynamics. J. Chem. Phys. 126, 155101. doi:10.1063/1.2714538

PubMed Abstract | CrossRef Full Text | Google Scholar

Culhane, K. J., Liu, Y., Cai, Y., and Yan, E. C. (2015). Transmembrane signal transduction by peptide hormones via family B G protein-coupled receptors. Front. Pharmacol. 6, 264. doi:10.3389/fphar.2015.00264

PubMed Abstract | CrossRef Full Text | Google Scholar

Dautzenberg, F. M., Mevenkamp, G., Wille, S., and Hauger, R. L. (1999). N-terminal splice variants of the type I PACAP receptor: isolation, characterization and ligand binding/selectivity determinants. J. Neuroendocrinol 11, 941–949. doi:10.1046/j.1365-2826.1999.00411.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Duan, J., Shen, D. D., Zhou, X. E., Bi, P., Liu, Q. F., Tan, Y. X., et al. (2020). Cryo-EM structure of an activated VIP1 receptor-G protein complex revealed by a NanoBiT tethering strategy. Nat. Commun. 11, 4121. doi:10.1038/s41467-020-17933-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Gottschall, P. E., Tatsuno, I., Miyata, A., and Arimura, A. (1990). Characterization and distribution of binding sites for the hypothalamic peptide, pituitary adenylate cyclase-activating polypeptide. Endocrinol. 127, 272–277. doi:10.1210/endo-127-1-272

CrossRef Full Text | Google Scholar

Graaf, C., Donnelly, D., Wootten, D., Lau, J., Sexton, P. M., Miller, L. J., et al. (2016). Glucagon-like peptide-1 and its class B G protein-coupled receptors: a long march to therapeutic successes. Pharmacol. Rev. 68, 954–1013. doi:10.1124/pr.115.011395

PubMed Abstract | CrossRef Full Text | Google Scholar

Grace, C. R., Perrin, M. H., Gulyas, J., Rivier, J. E., Vale, W. W., and Riek, R. (2010). NMR structure of the first extracellular domain of corticotropin-releasing factor receptor 1 (ECD1-CRF-R1) complexed with a high affinity agonist. J. Biol. Chem. 285, 38580–38589. doi:10.1074/jbc.M110.121897

PubMed Abstract | CrossRef Full Text | Google Scholar

Harmar, A. J., Fahrenkrug, J., Gozes, I., Laburthe, M., May, V., Pisegna, J. R., et al. (2012). Pharmacology and functions of receptors for vasoactive intestinal peptide and pituitary adenylate cyclase-activating polypeptide: IUPHAR Review 1. Br. J. Pharmacol. 166, 4–17. doi:10.1111/j.1476-5381.2012.01871.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Hollenstein, K., de Graaf, C., Bortolato, A., Wang, M. W., Marshall, F. H., and Stevens, R. C. (2014). Insights into the structure of class B GPCRs. Trends Pharmacol. Sci. 35, 12–22. doi:10.1016/j.tips.2013.11.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Hollenstein, K., Kean, J., Bortolato, A., Cheng, R. K., Doré, A. S., Jazayeri, A., et al. (2013). Structure of class B GPCR corticotropin-releasing factor receptor. Nature 499, 438–443. doi:10.1038/nature12357

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

Jazayeri, A., Doré, A. S., Lamb, D., Krishnamurthy, H., Southall, S. M., Baig, A. H., et al. (2016). Extra-helical binding site of a glucagon receptor antagonist. Nature 533, 274–277. doi:10.1038/nature17414

PubMed Abstract | CrossRef Full Text | Google Scholar

Jazayeri, A., Rappas, M., Brown, A. J. H., Kean, J., Errey, J. C., Robertson, N. J., et al. (2017). Crystal structure of the GLP-1 receptor bound to a peptide agonist. Nature 546 (7657), 254–258. doi:10.1038/nature22800

PubMed Abstract | CrossRef Full Text | Google Scholar

Jo, S., Kim, T., Iyer, V. G., and Im, W. (2008). CHARMM-GUI: a web-based graphical user interface for CHARMM. J. Comput. Chem. 29 (11), 1859–1865. doi:10.1002/jcc.20945

PubMed Abstract | CrossRef Full Text | Google Scholar

Kobayashi, K., Shihoya, W., Nishizawa, T., Kadji, F. M. N., Aoki, J., Inoue, A., et al. (2020). Cryo-EM structure of the human PAC1 receptor coupled to an engineered heterotrimeric G protein. Nat. Struct. Mol. Biol. 27 (3), 274–280. doi:10.1038/s41594-020-0386-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Kumar, S., Pioszak, A., Zhang, C., Swaminathan, K., and Xu, H. E. (2011). Crystal structure of the PAC1R extracellular domain unifies a consensus fold for hormone recognition by class B G-protein coupled receptors. PLoS ONE 6 (5), e19682. doi:10.1371/journal.pone.0019682

PubMed Abstract | CrossRef Full Text | Google Scholar

Lam, H. C., Takahashi, K., Ghatei, M. A., Kanse, S. M., Polak, J. M., and Bloom, S. R. (1990). Binding sites of a novel neuropeptide pituitary-adenylate-cyclase-activating polypeptide in the rat brain and lung. Eur. J. Biochem. 193 (3), 725–729. doi:10.1111/j.1432-1033.1990.tb19392.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, J., Abel, R., Zhu, K., Cao, Y., Zhao, S., and Friesner, R. A. (2011). The VSGB 2.0 model: a next generation energy model for high resolution protein structure modeling. Proteins 79 (10), 2794–2812. doi:10.1002/prot.23106

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, J., Jonsson, A. L., Beuming, T., Shelley, J. C., and Voth, G. A. (2013). Ligand-dependent activation and deactivation of the human adenosine A(2A) receptor. J. Am. Chem. Soc. 135 (23), 8749–8759. doi:10.1021/ja404391q

PubMed Abstract | CrossRef Full Text | Google Scholar

Liang, Y. L., Belousoff, M. J., Zhao, P., Koole, C., Fletcher, M. M., Truong, T. T., et al. (2020). Toward a structural understanding of class B GPCR peptide binding and activation. Mol. Cell 77 (3) 656–668. doi:10.1016/j.molcel.2020.01.012

PubMed Abstract | CrossRef Full Text | Google Scholar

Liang, Y. L., Khoshouei, M., Glukhova, A., Furness, S. G. B., Zhao, P., Clydesdale, L., et al. (2018). Phase-plate cryo-EM structure of a biased agonist-bound human GLP-1 receptor-Gs complex. Nature 555 (7694), 121–125. doi:10.1038/nature25773

PubMed Abstract | CrossRef Full Text | Google Scholar

Liao, C., May, V., and Li, J. (2019a). Assessment of conformational state transitions of class B GPCRs using molecular dynamics. Methods Mol. Biol. 1947, 3–19. doi:10.1007/978-1-4939-9121-1_1

CrossRef Full Text | Google Scholar

Liao, C., May, V., and Li, J. (2019b). PAC1 receptors: shapeshifters in motion. J. Mol. Neurosci. 68 (3), 331–339. doi:10.1007/s12031-018-1132-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Liao, C., de Molliens, M. P., Schneebeli, S. T., Brewer, M., Song, G., Chatenet, D., et al. (2019c). Targeting the PAC1 receptor for neurological and metabolic disorders. Curr. Top. Med. Chem. 19 (16), 1399–1417. doi:10.2174/1568026619666190709092647

CrossRef Full Text | Google Scholar

Liao, C., Zhao, X., Brewer, M., May, V., and Li, J. (2017). Conformational transitions of the pituitary adenylate cyclase-activating polypeptide receptor, a human class B GPCR. Sci. Rep. 7 (1), 5427. doi:10.1038/s41598-017-05815-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Metzner, P., Schütte, C., and Vanden-Eijnden, E. (2009). Transition path theory for Markov jump processes. Multiscale Model. Simul. 7 (3), 1192–1219. doi:10.1137/070699500

CrossRef Full Text | Google Scholar

Miller, B. R., McGee, T. D., Swails, J. M., Homeyer, N., Gohlke, H., and Roitberg, A. E. (2012), MMPBSA.py: an efficient program for end-state free energy calculations. J. Chem. Theor. Comput. 8 (9), 3314–3321. doi:10.1021/ct300418h

CrossRef Full Text | Google Scholar

Noé, F., Horenko, I., Schütte, C., and Smith, J. C. (2007). Hierarchical analysis of conformational dynamics in biomolecules: transition networks of metastable states. J. Chem. Phys. 126 (15), 155102. doi:10.1063/1.2714539

PubMed Abstract | CrossRef Full Text | Google Scholar

Noé, F., Schütte, C., Vanden-Eijnde, E., Reich, L., and Weikl, T. R. (2009). Constructing the equilibrium ensemble of folding pathways from short off-equilibrium simulations. Proc. Natl. Acad. Sci. United States 106 (45), 19011–19016. doi:10.1073/pnas.0905466106

CrossRef Full Text | Google Scholar

Pande, V. S., Beauchamp, K., and Bowman, G. R. (2010). Everything you wanted to know about Markov State Models but were afraid to ask. Methods 52 (1), 99–105. doi:10.1016/j.ymeth.2010.06.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Parthier, C., Kleinschmidt, M., Neumann, P., Rudolph, R., Manhart, S., Schlenzig, D., et al. (2007). Crystal structure of the incretin-bound extracellular domain of a G protein-coupled receptor. Proc. Natl. Acad. Sci. United States 104 (35), 13942–13947. doi:10.1073/pnas.0706404104

CrossRef Full Text | Google Scholar

Phillips, J. C., Braun, R., Wang, W., Gumbart, J., Tajkhorshid, E., Villa, E., et al. (2005). Scalable molecular dynamics with NAMD. J. Comput. Chem. 26 (16), 1781–1802. doi:10.1002/jcc.20289

PubMed Abstract | CrossRef Full Text | Google Scholar

Pioszak, A. A., and Xu, H. E. (2008). Molecular recognition of parathyroid hormone by its G protein-coupled receptor. Proc. Natl. Acad. Sci. United States 105 (13), 5034–5039. doi:10.1073/pnas.0801027105

PubMed Abstract | CrossRef Full Text | Google Scholar

Prinz, J. H., Wu, H., Sarich, M., Keller, B., Senne, M., Held, M., et al. (2011). Markov models of molecular kinetics: generation and validation. J. Chem. Phys. 134 (17), 174105. doi:10.1063/1.3565032

PubMed Abstract | CrossRef Full Text | Google Scholar

Ressler, K. J., Mercer, K. B., Bradley, B., Jovanovic, T., Mahan, A., Kerley, K., et al. (2011). Post-traumatic stress disorder is associated with PACAP and the PAC1 receptor. Nature 470 (7335), 492–497. doi:10.1038/nature09856

PubMed Abstract | CrossRef Full Text | Google Scholar

Schäfer, H., Zheng, J., Morys-Wortmann, C., Fölsch, U. R., and Schmidt, W. E. (1999). Structural motifs of pituitary adenylate cyclase-activating polypeptide (PACAP) defining PAC1-receptor selectivity. Regul. Pept. 79 (2-3), 83–92. doi:10.1016/s0167-0115(98)00147-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Siu, F. Y., He, M., de Graaf, C., Han, G. W., Yang, D., Zhang, Z., et al. (2013). Structure of the human glucagon class B G-protein-coupled receptor. Nature 499 (7459), 444–449. doi:10.1038/nature12393

PubMed Abstract | CrossRef Full Text | Google Scholar

Song, G., Yang, D., Wang, Y., de Graaf, C., Zhou, Q., Jiang, S., et al. (2017). Human GLP-1 receptor transmembrane domain structure in complex with allosteric modulators. Nature 546, 312. doi:10.1038/nature22378

PubMed Abstract | CrossRef Full Text | Google Scholar

Sun, C., Song, D., Davis-Taber, R. A., Barrett, L. W., Scott, V. E., Richardson, P. L., et al. (2007). Solution structure and mutational analysis of pituitary adenylate cyclase-activating polypeptide binding to the extracellular domain of PAC1-Rs. Proc. Natl. Acad. Sci. United States 104 (19), 7875–7880. doi:10.1073/pnas.0611397104

CrossRef Full Text | Google Scholar

Vaudry, D., Falluel-Morel, A., Bourgault, S., Basille, M., Burel, D., Wurtz, O., et al. (2009). Pituitary adenylate cyclase-activating polypeptide and its receptors: 20 Years after the discovery. Pharmacol. Rev. 61 (3) 283–357. doi:10.1124/pr.109.001370

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, J., Song, X., Zhang, D., Chen, X., Li, X., Sun, Y., et al. (2020). Cryo-EM structures of PAC1 receptor reveal ligand binding mechanism. Cell Res. 30 (5), 436–445. doi:10.1038/s41422-020-0280-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Weinan, E., and Vanden-Eijnden, E. (2006). Towards a theory of transition paths. J. Stat. Phys. 123, 503–523. doi:10.1007/s10955-005-9003-9

CrossRef Full Text | Google Scholar

Wootten, D., Simms, J., Miller, L. J., Christopoulos, A., and Sexton, P. M. (2013). Polar transmembrane interactions drive formation of ligand-specific and signal pathway-biased family B G protein-coupled receptor conformations. Proc. Natl. Acad. Sci. United States 110 (13), 5211–5216. doi:10.1073/pnas.1221585110

CrossRef Full Text | Google Scholar

Yang, L., Yang, D., de Graaf, C., Moeller, A., West, G. M., Dharmarajan, V., et al. (2015). Conformational states of the full-length glucagon receptor. Nat. Commun. 6, 7859. doi:10.1038/ncomms8859

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, C., and Ma, J. (2010). Enhanced sampling and applications in protein folding in explicit solvent. J. Chem. Phys. 132 (24), 244101. doi:10.1063/1.3435332

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, H., Qiao, A., Yang, D., Yang, L., Dai, A., de Graaf, C., et al. (2017a). Structure of the full-length glucagon class B G-protein-coupled receptor. Nature 546 (7657), 259–264. doi:10.1038/nature22363

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, Y., Sun, B., Feng, D., Hu, H., Chu, M., Qu, Q., et al. (2017b). Cryo-EM structure of the activated GLP-1 receptor in complex with a G protein. Nature 546 (7657), 248–253. doi:10.1038/nature22394

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, H., Qiao, A., Yang, L., Van Eps, N., Frederiksen, K. S., Yang, D., et al. (2018). Structure of the glucagon receptor in complex with a glucagon analogue. Nature 553 (7686), 106–110. doi:10.1038/nature25153

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhao, S., Zhu, K., Li, J., and Friesner, R. A. (2011). Progress in super long loop prediction. Proteins 79 (10), 2920–2935. doi:10.1002/prot.23129

PubMed Abstract | CrossRef Full Text | Google Scholar