Chenyi Liao
Jacob M. Remington
Victor May
Jianing Li- 1Department of Chemistry, University of Vermont, Burlington, VT, United States
- 2State Key Laboratory of Molecular Reaction Dynamics, Dalian Institute of Chemical Physics, Dalian, China
- 3Department of Neuroscience, University of Vermont, Burlington, VT, United States
The related neuropeptides PACAP and VIP, and their shared PAC1, VPAC1 and VPAC2 receptors, regulate a large array of physiological activities in the central and peripheral nervous systems. However, the lack of comparative and molecular mechanistic investigations hinder further understanding of their preferred binding selectivity and function. PACAP and VIP have comparable affinity at the VPAC1 and VPAC2 receptor, but PACAP is 400–1,000 fold more potent than VIP at the PAC1 receptor. A molecular understanding of the differing neuropeptide-receptor interactions and the details underlying the receptor transitions leading to receptor activation are much needed for the rational design of selective ligands. To these ends, we have combined structural information and advanced simulation techniques to study PACAP/VIP binding selectivity, full-length receptor conformation ensembles and transitions of the PACAP/VIP receptor variants and subtypes, and a few key interactions in the orthosteric-binding pocket. Our results reveal differential peptide-receptor interactions (at the atomistic detail) important for PAC1, VPAC1 and VPAC2 receptor ligand selectivity. Using microsecond-long molecular dynamics simulations and the Markov State Models, we have also identified diverse receptor conformational ensembles and microstate transition paths for each receptor, the potential mechanisms underlying receptor open and closed states, and the interactions and dynamics at the transmembrane orthosteric pocket for receptor activation. These analyses reveal important features in class B GPCR structure-dynamics-function relationships, which provide novel insights for structure-based drug discovery.
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Å3.
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) (dECD-ECL1) and extracellular loop 3 (ECL3) (dECD-ECL3), respectively. Detailed constructions and validation of MSM were described in the Supplementary Material.
Data Analysis
Binding free energy (ΔGb) 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. ΔGb 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. (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 (ΔGb) of PACAP or VIP to the ECD of PAC1null, PAC1s, VPAC1 and VPAC2 receptors, as well as the ΔGb contribution per residue. Our results (Figure 2A) show that PACAP binding to PAC1null receptor is the most favorable among receptor and peptide combinations. Specifically, ΔGb of PACAP is around 1.4 times of ΔGb 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 (Kd ≈ 0.5 nM) over VIP (Kd >500 nM) (Vaudry et al., 2009), which corresponds to that PACAP-binding free energy (−13.2 kcal/mol with Kd converted by ΔGb = RTlnKd, 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 (Kd ≈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 ΔGb values for comparison with our calculation estimations. Therefore, here we compare the relative values of ΔGb for the different systems as opposed to the absolute magnitude which differs from prior experiments.
FIGURE 2. Binding selectivity and recognition of PACAP and VIP on PAC1null, PAC1s, VPAC1 and VPAC2 receptors. (A) Binding free energy (ΔGb) of PACAP and VIP on PAC1null, PAC1s, VPAC1 and VPAC2 receptors respectively by averaging the ΔGb 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 V19xxY22LxxV/I26 sequence appears to contribute to the hydrophobic interactions with PAC1null/PAC1s and VPAC1/2 receptor residues I61PAC/L70VPAC1/I59VPAC2; the PACAP/VIP R14K15 residues contribute to hydrogen bonding with the ECD (Supplementary Figure S10). The β3-β4 loop in both PAC1null/PAC1s and VPAC2 receptor (environs L80FxxF84,PAC and V78FxxF82,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 F131xxxxF136,PAC1null, F110xxxxF115,PAC1s, Y118xxxxxL124,VPAC1, and V106xxxxY111,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–140PAC1null play major role in PACAP binding, such as M72, L80FxIF84xP86, VW90, I95, L106xLxxM111, F127, F131, F136, and Y139 (Figure 2B); only E142 and D145 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–87PAC1null such as E79, R82, N85 and D87. Lastly, PACAP interacts with residues 92–130PAC1null extensively, while VIP interacts with few. This also supports the finding that another short variant with deletion of residues 53–109PAC1null 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 PACAP1-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 D3PACAP-R199PAC1 interacting pair (Figure 2C). The depth of peptide N-terminus into the 7TM is not significantly changed, but H1PACAP 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 PACAP1-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), D3PACAP and PAC1null R1992.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 H1PACAP and the surrounding receptor K206/D207ECL1, E374/E380/E385ECL3, and K1541.40 decreases gradually and appears to form hydrogen bond (H-bond) networks with an average small distance of ∼5 Å. The hydrogen bonding between the D3PACAP sidechain and H1PACAP amino terminus weakens during adaptive tempering resulting in the ability for D3PACAP to form a new pairing with receptor R1992.60. As a result of D3PACAP-R1992.60 interactions, the average distance between R1992.60 and N2403.43, H3656.52 and Q3927.49 increases significantly (Figure 3B, middle panel), indicating weakened contacts between TM2, TM3, TM6, and TM7 (Figure 4D). At the receptor intracellular face, K180ICL1 transitions from a solvent orientation to face E247/R185, and decreasing K180ICL1-E2473.50/R185ICL1 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–13PACAP (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. 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. 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, E374ECL3 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 R1992.60 and N2403.43, H3656.52 and Q3927.49, and the distance between the R350ICL3-E406H8 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 D3PACAP-R1992.60 was also found in our simulations, although the N-terminus of PACAP underwent conformational rearrangements in our simulations. Together, the C-termligand-ECD and N-termligand-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 D3PACAP–R199PAC1null 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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
Keywords: G protein-coupled receptor, membrane protein, signaling, molecular dynamics, conformational transition, ligand selectivity
Citation: Liao C, Remington JM, May V and Li J (2021) Molecular Basis of Class B GPCR Selectivity for the Neuropeptides PACAP and VIP. Front. Mol. Biosci. 8:644644. doi: 10.3389/fmolb.2021.644644
Received: 21 December 2020; Accepted: 16 February 2021;
Published: 25 March 2021.
Edited by:
Joanna Trylska, University of Warsaw, PolandReviewed by:
Yinglong Miao, University of Kansas, United StatesShuguang Yuan, Institute of Synthetic Biology, Shenzhen Institutes of Advanced Technology, China
Copyright © 2021 Liao, Remington, May and Li. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Jianing Li, jianing.li@uvm.edu