Vibrational energy landscapes and energy flow in GPCRs: comparison between class A and class B GPCRs using all atom and coarse-grained models

Abstract

We compare vibrational energy landscapes and dynamics of the glucagon-like peptide 1 receptor (GLP-1R), a class B G-protein coupled receptor (GPCR), with corresponding properties of a class A GPCR studied previously. Energy flow in GLP-1R is computed by molecular dynamics (MD) simulations in active and inactive states using all atom (AA) and coarse-grained (CG) models. Based on the MD data, we construct and analyze the vibrational energy landscape of GLP-1R, focusing on the relative free energy of each residue and the minimum free energy barriers for energy transfer between them. We find that prolines and glycines, which contribute to GLP-1R plasticity and function, are bottlenecks to energy transport along the backbone, causing diversion of energy through alternative pathways via nearby noncovalent contacts. The probability distributions for the energy transfer time between numerous pairs of amino acids are computed, revealing pathways for energy transport that include noncovalent contacts connecting different helices. These distributions and mean energy transfer times are similar for the AA and CG models, validating the CG representation for future applications.

---

Dynamically coupled residues that can be quite distant from each other contribute to enzyme catalysis^1–3 and allostery,^4–7 motivating many studies to identify networks of residues that facilitate energy flow.^8,9 Investigations of energy transfer in proteins include experimental and computational studies probing energy dynamics in response to chemical reactions.^10–12 Some of these studies have addressed vibrational energy relaxation and flow in receptor proteins, notably rhodopsin in response to photoexcitation.^13–15 Rhodopsin is a G-protein coupled receptor (GPCR). Recent studies have also investigated energy flow in other GPCRs and connections to change in structure and structural dynamics upon activation or deactivation by binding of ligands.^16–19

Rhodopsin is a class A GPCR, the most abundant and best studied class, which also includes the adrenaline-binding receptor, β₂ adrenergic receptor (β₂AR). The structure of all GPCRs features 7 transmembrane (TM) helices, but members of different GPCR classes are distinguished by their motifs, regions with conserved residues, many of which undergo structural change during activation.²⁰ While the vibrational energy landscape and energy transport have been studied for β₂AR,¹⁹ a class A GPCR, much less is known about energy transport in GPCRs of other classes. In this report, we present results of a computational study of the vibrational energy landscape and energy transport in the glucagon-like peptide 1 receptor (GLP-1R), a class B GPCR, in active and inactive states. With the new results for GLP-1R, we provide a comparison between properties that mediate the flow of energy in class A and class B GPCRs.

There is a second aim of the present report. Protein dynamics occur over a wide range of time scales, including fast sub-picosecond to nanosecond times for vibrational dynamics and energy flow^8,9 with folding times of microseconds and beyond.²¹ While the fast times for vibrational energy transfer are favorable for computational study, noise in the dynamics often means that long computational studies are still needed to identify patterns in the flow of energy, even in small proteins.^22,23 For larger proteins or protein complexes, e.g., a GPCR complexed with the G-protein, coarse graining offers a means to reduce computational cost, but to date there have been no detailed studies comparing vibrational energy flow in sizable proteins using all atom (AA) and coarse-grained (CG) models. In this report, we compare vibrational energy dynamics of GLP-1R using AA and CG models.

GLP-1R (Fig. 1) regulates bone turnover, blood glucose control, and cardiovascular development,²⁴ and is a potential drug target of type 2 diabetes, neurodegeneration, and cardiovascular diseases.²⁵ Various endogenous agonists activate GLP-1R, some of which are either under development or approved for the treatment of type 2 diabetes.^25–27 Here, we present results of molecular dynamics (MD) simulations using AA and CG models for GLP-1R in an active state, where the agonist is a peptide ligand, and in an inactive state, where no ligand is present. In a previous computational study of GLP-1R using the same AA model, structural and dynamic properties of the protein mediating energy flow were identified,¹⁷ and the contribution of transmembrane water in stabilizing the active state was analyzed. Elucidation of the connections between the protein structure and dynamics to energy flow, a goal of many previous studies,^{8,9,22,23,28–43} has shown that measurement of energy transfer rates in a protein can yield information about underlying structural and dynamic properties^18,44–48 that contribute to the dynamics and rates of chemical reactions in proteins and allostery.^1–3,12,49

Fig. 1 GLP-1R is shown in (a) inactive and (b) active states with the peptide ligand in blue. The active state is also shown in (c) with the beta turn, discussed in text, colored cyan and prolines and glycines indicated, also discussed in the text. Views from the cytoplasmic region are shown in the (d) inactive and (e) active state. Note that TM6 is farther from TM3 in (e) than in (d). Prolines and glycines are helix destabilizers. (f) Kinks appear at these residues in transmembrane (TM) helices 6 and 7, which change angle between inactive (gray) and active (green) states. Their presence in the TM region contributes to GPCR flexibility and function.

We recently computed disconnectivity graphs^50–52 to represent the vibrational energy landscape governing energy flow in β₂AR in both active and inactive states,¹⁹ using data from MD simulations generated in earlier work on this protein.^16,18 Disconnectivity graphs provide a visualization of how free energy minima are organized in the energy landscape based on the lowest barriers between them.^50–52 Here we employ a construction representing the vibrational energy landscape of a protein using a measure of the free energy associated with each amino acid, related to its size. The minimum free energy barriers to energy transfer between residue pairs¹⁹ are obtained by converting the calculated rates using the Eyring–Polanyi relation.^53,54 Hence, the depth of the branches on the vertical scale of the disconnectivity graphs, associated with amino acid residues, encode their relative free energies, while the barriers return the calculated rates by construction.

The disconnectivity graphs computed for β₂AR provided new insight into the contributions of composition, structure and dynamics to energy flow. Noncovalent contacts between TM helices, such as hydrogen bonds and van der Waals interactions, can facilitate energy transfer between them. By construction, noncovalent contacts for which energy transfer rates are relatively large correspond to relatively low barriers between TM helices.

Shortcuts for energy flow in sequence via noncovalent contacts such as hydrogen bonds provide alternative pathways if there are bottlenecks to energy flow along the main chain. In the study of β₂AR, prolines and glycines were identified as bottlenecks to energy transfer along the protein backbone. Prolines and glycines destabilize helices,⁵⁵ and their presence in TM helices gives rise to kinks (Fig. 1), which contribute to flexibility and function.^56–59 The rate of energy transfer across a noncovalent contact depends both on the nature of the contact and on fluctuations in its length. Rates of energy transfer are typically faster across hydrogen bonds than across van der Waals contacts. Both are short range contacts. When either contact is present, the rate of energy transfer is fastest if fluctuations in the length of the contact are small,^18,60 corresponding to a lower barrier, which can be read from the disconnectivity graph.

Both structural and dynamic changes that occur upon binding of agonists contribute to allosteric regulation of GPCRs,^61–68 and these changes are related to differences in the flow of energy through the protein upon activation. Class A GPCRs, such as β₂AR, have been studied for some time, while structures for class B GPCRs such as GLP-1R have only become known more recently.^69–78 It is interesting to determine features found in the vibrational energy landscape of β₂AR shared by GPCRs of different classes and the study of GLP-1R allows us to make such a comparison, in addition to a comparison of results obtained using AA and CG models.

In the following section, we present the methods used in the computational study of vibrational energy flow in GLP-1R. We then compare energy transfer rates between residues of GLP-1R obtained using AA and CG models, disconnectivity graphs constructed for these systems, and first passage time distributions computed for various energy source-sink combinations across noncovalent contacts. We compare features of the vibrational energy landscape and energy flow computed for GLP-1R with previous results¹⁹ computed for β₂AR. We also compare mean first passage times computed for transport between pairs of residues (termed sources and sinks) where prolines or glycines are present or absent. This comparison provides another test for CG models in capturing how vibrational energy flows in the protein.

Methods

The inactive state structure was taken from the protein data bank (PDB) entry 6LN2⁷¹ and the active structure, which contains a peptide agonist, from PDB 5VAI.⁷⁴ We modeled all missing residues in both states, mainly in the extracellular domain (ECD) and TM helix 1 (TM1) using Modeller9.23,⁷⁹ which gives the same homology of 393 amino acid sequences in each state. Each system was initiated in a rectangular box of 170 1-palmitoyl-2-oleoyl-sn-glycero-3-phosphocholine (POPC) lipid molecules, with the protein in the center, solvated with 21 000 water molecules. The systems were neutralized using Na⁺ and Cl⁻ ions at 0.15 M NaCl concentration. We used the TIP3P,⁸⁰ Lipid17⁸¹ and AMBER ff14SB⁸² representations for water, lipids, and protein, respectively.

Details of the MD simulations using the AA model have been provided elsewhere.¹⁷ Briefly, we performed 10 000 steps of energy minimization using the steepest-descent method and another 10 000 steps with the conjugate gradient method. The NVT ensemble was heated starting from 0.1 K to 300 K for 1 ns and held at 300 K for an additional 1 ns using the Berendsen⁸³ thermostat. Positional constraints were applied to backbone atoms with a force constant of 1 kcal mol⁻¹ Å⁻². The SHAKE algorithm was applied to all bonds containing hydrogens. We performed equilibration using 2 fs time steps for 10 ns with position restraints and another 20 ns without restraints. The production simulations of 600 ns were performed for each system using 2 fs time steps in the NPT ensemble saving the trajectory files every 1 ns for subsequent microcanonical simulations. Energy currents, G, between residue–residue pairs of GLP-1R were computed from the trajectories of the 500 to 600 ns NPT simulations starting with 50 structures obtained every 2 ns. Each one was simulated for 150 ps in the NVE ensemble with an integration time step of 0.5 fs and an Ewald sum tolerance of 10⁻⁷ to reduce the energy drift. All energy current calculations were computed using the CURrent calculation for proteins (CURP) version 1.2.1⁸⁴ developed by Yamato and coworkers.

CG simulations were carried out for these systems following the protocol developed by Pantano and coworkers.^85,86 CG forcefield models preclude an accurate estimation of conformational flexibility if the secondary and tertiary structures are assumed to be fixed,^85,86 which could affect the study of energy transfer. To investigate this potential systematic error, we began by mapping results from the AA force field to those obtained with the CG model. To avoid the overlap of protein and lipid, lipid molecules within 3.5 Å of the protein were removed, following a standard protocol.^85,86 The system was solvated with an explicit CG water model, where again care was taken so that water molecules and other parts of the system do not overlap.^85,86 The unnecessary water near the hydrophobic region was removed.

The SIRAH force field was used,^85,86 with parameters generated and simulations performed using the AMBER software package. The energy minimization was performed in two steps. First, the side chain atoms were energy minimized with positional restraints for the backbone atoms with a force constant of 2.4 kcal mol⁻¹ Å⁻² for 10 000 steps, and then the whole system was energy minimized for another 10 000 steps without positional restraints. The system was heated for 500 ps from 1 K to 300 K and maintained at 300 K for an additional 500 ps. After heating, we performed equilibrium simulations for 2 ns with an integration time of 4 fs using the positional restraint for the protein backbone atoms with a force constant of 1.0 kcal mol⁻¹ Å⁻². Finally, we performed 1 μs simulations with an integration time step of 20 fs.

Upon analyzing the structures during the simulations, we observed that after 500 ns, TM6 of the active state turned towards TM3 in the coarse-grained simulations (see SI), yielding a structure resembling more the inactive state structure than the active state. The change in structure of the active state using the SIRAH CG model is not unique to that model. Separate simulations carried out using the MARTINI force field^87,88 for the active state produced a similar displacement of TM6 towards TM3 after about 500 ns. The CURP calculations were carried out for times where the active state structure was maintained and employed 50 structures, using the same protocol as described above for the AA model, considering the time range of 200 to 400 ns and taking trajectories at regular spacings.

We refer to the energy current, G, obtained using the AA and CG models as G_AA and G_CG, respectively. We need to introduce scaling factors to compare them; we expect these parameters to be transferable to different proteins. To scale G_CG, we systematically grouped residue pairs based on their sequence separation. Residue–residue separations of 1, 2, 3, and 4 were each treated as individual groups, independent of their G values. For longer-range interactions (≥5 residues apart in sequence), residue pairs were further classified into six bins based on their G_AA values in units of (kcal mol⁻¹)² ps⁻¹: G_AA > 100, 75–100, 50–75, 25–50, 1–25, and G_AA < 1. This classification yielded a total of ten distinct groups. For each group, we computed a scaling factor, A, listed in Table S1, defined as the ratio of the average G_AA value to G_CG value, i.e., . These scaling factors were then used to adjust the CG model values as G = G_CG × A. We expect the scaling factors to be applicable to other systems when comparing results using the same AA and CG models.

G is proportional to the rate of energy transfer, w, between residues. Rates of vibrational energy transfer, w, have also been computed from nonequilibrium simulations by Stock and coworkers.²² On average, the value of G, with units (kcal mol⁻¹)² ps⁻¹, can be converted to w using w = G/(187.93 (kcal mol⁻¹)²), as described in a previous study,¹⁹ and we use this conversion to produce rate constants from values of G for residue pairs of GLP-1R. The rate constants, moreover, satisfy detailed balance, i.e., the rate of energy transfer from residue i with f_i degrees of freedom to j with f_j degrees of freedom, w_ji, is w_ji = w_ij(f_j/f_i). We note that some of the water molecules in the transmembrane region of GPCRs stabilize the active state^17,66 and can also provide pathways for energy transport.^7,17 The rate constants for energy transfer between residue pairs were computed with the water present so that disconnectivity graphs computed with the rate constants account for the full system.

Disconnectivity graphs were constructed using the matrices of rate constants, as in ref. 89 and 90, to identify free energy minima for each amino acid and barriers between them. The relative free energies of the residues and the barriers between them are translated to reproduce the rates and the corresponding equilibrium distribution of the linear master equation. We compute first passage time (FPT) distributions within the same master equation representation for energy transfer between all residue pairs following the procedure detailed in ref. 89 and 91. The energy source is an amino acid with excess vibrational energy, due to collision or photoexcitation. The energy sink is a selected amino acid to which energy from the source flows. The first passage time is the time taken for energy to first reach the sink from the source, determined by the rate constants for energy transfer between all pairs of amino acids. Since there are an infinite number of pathways through which energy can flow between an energy source and an energy sink, we obtain a distribution of first passage times: some paths are much more likely than others. Dijkstra's algorithm⁹² is applied to identify shortest paths between an energy source and sink in terms of the transfer time. Here the metric for the path length is defined using edge weights −ln B_ij, where B_ij is the branching probability from residue j to i. This choice selects the path corresponding to the largest term in the overall rate if intervening minima are treated in steady state.⁹³

Results and discussion

In Fig. 2, we plot G for all residue pairs, referred to as energy exchange networks (EENs), for GLP-1R in active and inactive states computed using both AA and CG models. EENs for the inactive state are plotted in Fig. 2(a) and (b) and those for the active state are shown in Fig. 2(c) and (d). The results obtained with the AA and CG models are plotted in the left and right columns, respectively. Data plotted in red represent the largest values of G, G > 50 (kcal mol⁻¹)² ps⁻¹, in cyan the next largest values, 50 (kcal mol⁻¹)² ps⁻¹ > G > 25 (kcal mol⁻¹)² ps⁻¹, with smaller, though still significant, values in blue, 25 (kcal mol⁻¹)² ps⁻¹ > G > 1 (kcal mol⁻¹)² ps⁻¹.

Fig. 2 Energy exchange networks (EENs) for GLP-1R obtained using AA (left column) and CG models (right column). In the top row, panels (a) and (b) represent inactive states and in the bottom row, panels (c) and (d) represent active states. Red indicates residue pairs with G > 50 (kcal mol⁻¹)² ps⁻¹, cyan indicates 50 (kcal mol⁻¹)² ps⁻¹ > G > 25 (kcal mol⁻¹)² ps⁻¹, and blue indicates 25 (kcal mol⁻¹)² ps⁻¹ > G > 1 (kcal mol⁻¹)² ps⁻¹. There are 393 residues in each state and the active state has an additional peptide ligand corresponding to residues 394–424. Helices and the ECD are labeled at the edges of the panels.

Values of G are relatively large near the diagonal, where covalent interactions and hydrogen bonding along helices contribute to the EENs. Energy transfer due to interactions between helices appears near and perpendicular to the diagonal for neighboring TM helices, and away from the diagonal due to interactions between TM helices that are farther from each other in sequence. Significant differences between the latter contributions for the active and inactive states appear, due to the opening of the structure upon activation (Fig. 1), which rearranges the interactions, and therefore values of G, between TM3 and TM6, TM5 and TM6, TM6 and TM7, among others. Differences between the active and inactive states in the TM region obtained with the AA model are captured well by the CG representation.

There are differences between the AA and CG results within the extracellular domain (ECD), which corresponds to the first 110 residues. The ECD contains mainly loops and is generally more flexible than the helical TM region. Due to its greater structural variability, we expect to find larger differences in the ECD compared to the TM region.

Disconnectivity graphs for the inactive and active states of GLP-1R are plotted in Fig. 3. Each branch corresponds to an amino acid, and the height of each branch is proportional to the free energy of the amino acid, which is related to its size (the number of vibrational degrees of freedom). An energy scale bar of 1 kcal mol⁻¹ appears in Fig. 3(a). Deeper wells in the energy landscape, represented by longer branches, correspond to larger amino acids, while shallower ones appear for smaller residues. The branches, corresponding to particular amino acids, are connected via the lowest barriers where energy transfer is the fastest. The height of the barriers corresponds to the free energy of the transition state, which is a mapping of the calculated rate via the Eyring–Polanyi formulation.^53,54

Fig. 3 Disconnectivity graphs of GLP-1R computed using the AA model in (a) peptide-bound active and (b) inactive states, and the CG model for (c) peptide-bound active state and (d) inactive states. The transmembrane helices are colored and labeled in the figure. An energy scale bar of 1 kcal mol⁻¹ is indicated in (a).

If the fastest rate of energy transfer occurs between amino acids in sequence, a series of branches connected to nearest neighbors would be visible. However, because of noncovalent contacts such as hydrogen bonds, energy transfer between residues more distant in sequence may occur more rapidly than through the covalent bond network between neighbors. In this case, branches corresponding to amino acids far apart in sequence are connected in the disconnectivity graph. These shortcuts are highlighted by coloring the branches according to the helix in which each residue is located.

Because function of GPCRs is mediated largely by the TM helices, it is convenient to group the branches into helices and examine competition between energy flow along a helix and between helices. We focus on the patterns in energy flow along TM helices, which are labeled and indicated by the different colors, along with some of the specific residues. We also discuss helix 8, located in the cytoplasm and not a TM helix, which is in contact with some of the TM helices.

The disconnectivity graphs reveal a number of interesting features. First, we note that the lowest barriers separating minima are modest, of order 1 kcal mol⁻¹, as indicated by the energy scale bar in Fig. 3(a), consistent with the picosecond time scale for energy transfer between nearby amino acids. The relatively small barriers allow for facile flow of energy along the main chain of the protein. However, there is a noticeable variation in barrier heights, producing clusters of residues in the graph in which energy flow is particularly fast.

In both active and inactive states, computed using both AA and CG models, TM1 (red) consistently separates into two distinct clusters of branches, with Gly123 acting as a dividing residue. Each amino acid is represented by a branch, and groups of amino acids through which energy flows readily appear as a cluster of branches. We see that amino acids 110–122 of TM1 are clustered together, as are 123–139. Glycine acts as an energy flow bottleneck due to its small size, as seen in our previous analysis of β₂AR.¹⁹ When there are also noncovalent contacts present near the glycine, coupling different helices, the residues of one helix are seen to split into different clusters of branches, mediated by glycine.

TM2 was found to cluster into a largely single region in AA active, AA inactive, and CG inactive models, with three identifiable but closely grouped branches. Some interactions with residues outside the TM regions are apparent. In contrast, in the CG active model, these three groups were also separated and now mix, not only with residues in loops outside the TM region, but also with helix 8, which is in the intracellular region (Fig. 1). Since helix 8 can move in the cytoplasmic region, it may come into contact with one or several TM helices during a simulation, as observed for the active state using the CG model. TM2 notably lacks glycine and proline, both of which have been identified as energy transport bottlenecks.¹⁹

TM3 consistently forms two primary clusters of branches in all panels of Fig. 3, Gly188–Glu219 and Gly220–Leu227, with Gly220 serving as a clear boundary. The Gly188–Glu219 segment is intact with one cluster of branches in the active state. In the inactive state, this sequence of residues splits into smaller clusters in both the AA and CG models, because of interactions between this portion of TM3 and portions of TM4 and TM6.

Similar patterns appear for the remaining TM helices. TM4 is divided into three clusters in both states for both models, with Gly245 and Gly257 contributing to these separations. In TM5, fragmentation of clusters of branches occurs at Pro284 and Gly290 for both active and inactive states using both AA and CG models. TM6 displays clusters of branches separated by Gly333 and Pro330, which together form the kink in TM6 (Fig. 1). TM7 is split into two major clusters by Gly367, again at a kink in the helix (Fig. 1), with some difference in minor clusters of branches for the different models. Finally, helix 8 exhibits a range of clustering structures due to its intracellular location and flexible movement during the MD simulations, as noted above.

In a previous study,¹⁹ we constructed disconnectivity graphs for β₂AR, a class A GPCR, in active and inactive states using rate matrices obtained from MD simulations using the same AA force field as for GLP-1R. We now compare some of the results for β₂AR with those for GLP-1R. Glycine and proline residues act as bottlenecks to energy flow along the backbone in both systems, providing a route for energy transfer to other residues far apart in sequence via direct noncovalent interactions. Class A GPCRs, such as β₂AR, have several conserved motifs including Pro-Ile-Phe (PIF), Cys-Trp-xx-Pro (CWxP), Asn-Pro-xx-Tyr (NPxxY) in TM5, TM6, and TM7, where the proline residues are located. In addition, there are proline residues in intracellular loops 2 and 4, TM2 and TM4.

Class B GPCRs, such as GLP-1R, have two conserved motifs in which only one motif, Pro-x-x-Gly (PxxG), in TM6, consists of proline and glycine, which help to form the kink and shift the lower half of TM6 upon activation (Fig. 1). The HETx motif of class B GPCRs is spanned by several helices including TM7. Gly367 of TM7 gives rise to the formation of a sharp kink (Fig. 1) and hence the graph is split into different clusters.

For β₂AR, we identified a large barrier between TM1–3 and TM4–7, which appears in both the active and inactive states, caused by Pro110 in the intracellular loop connecting TM3 to TM4.¹⁹ No such barrier appears in this region for GLP-1R as there is no proline between TM3 and TM4.

In GLP-1R, there is a proline at the edge of TM1 separating a series of residues in the ECD from TM1, a separation that is not seen in the disconnectivity graphs for β₂AR, which does not have a proline in this location. In GLP-1R, there are glycines in the loops connecting TM1–2, TM4–5, and TM6–7. These glycines separate these TM helices, which is again different from β₂AR. For both β₂AR and GLP-1R, proline causes the fragmentation of clusters of branches belonging to TM helices. For β₂AR, this effect is seen for TM2, TM4, TM5, TM6, and TM7, and for GLP-1P, we find this pattern for TM5 and TM6.

Residues of the beta turn in the ECD of GLP-1R appear fragmented due to the prolines and glycines at each end of the turn, a feature absent for β₂AR, which lacks this beta turn in the ECD. Finally, for β₂AR, we see that the residues of helix 8 are well separated from other residues by a large barrier, due to the presence of Pro270 in the intracellular loop separating TM7 from helix 8. In GLP-1R, helix 8 is connected via a loop with TM7 without a barrier, due to the absence of either proline or glycine.

We turn now to first passage time (FPT) distributions for vibrational energy transfer, introduced in the Methods section, which have been computed for GLP-1R for a variety of energy sources and energy sinks. An energy source is an amino acid with excess energy, and an energy sink is a chosen amino acid to which the energy flows. Any pair of residues can be selected as the source and sink, and a separate first passage time distribution is calculated for each choice. For energy transport along the helices of β₂AR, we found that the FPT distributions look quite similar to one another.¹⁹ We therefore focus on FPT distributions for energy transport via noncovalent contacts. We consider energy sources and energy sinks on or near different helices that are linked by noncovalent interactions yielding the largest rates for energy transfer between helices. We plot results in Fig. 4 using the AA model as solid curves and the CG model as dashed curves for selected combinations. Noncovalent contacts between different helices for which the largest rates of energy transfer were obtained using the AA model have been chosen for the examples we discuss.

Fig. 4 FPT distributions plotted as P(ln t) vs. ln t, for energy source and sink combinations obtained from the all atom model (solid curves) and coarse-grained model (dotted curves) for active (blue) and inactive (red) states, respectively. Source-sink combinations are (a) Asn276–Phe341, (b) Cys198–Val259, (c) Val218–Ile280, and (d) Leu153–Ala210.

Each FPT distribution in Fig. 4 exhibits a large peak at around 100 ps. We observe this long-time peak in all the FPT distributions we have computed, including those for β₂AR.¹⁹ This long time peak corresponds to the random flow of energy among roughly 400 chemical groups and is unrelated to the structure of the disconnectivity graph.¹⁹ It apparently takes around 100 ps for energy to flow randomly to the sink, which far exceeds the time for excess vibrational energy to relax into the environment of the protein.^22,23 We therefore focus on peaks at shorter times, which sometimes appear as a shoulder at the short time side of the long-time peak. The peaks at shorter times correspond to the faster pathways for energy flow revealed by the disconnectivity graphs.

In Fig. 4(a), we plot the FPT distributions for the source-sink combination Asn276–Phe341, between which energy flows via the shortcut Arg282–Asp344 between TM5 and the loop just beyond TM6. The FPT distribution obtained with the CG model is similar to the FPT distribution obtained with the AA model. The peak at short times seen for the inactive state is clear for both the AA and CG data, and the shift to longer times for the active state due to a weaker interaction between Arg282 and Asp344 upon activation is captured well by both the AA and CG models.

A similar trend is found for the Cys198–Val259 source-sink combination plotted in Fig. 4(b). The slower energy transfer rate through the shortcut between TM3 and TM4 in the active state, as seen by a shift to longer times in the short-time peak of the FPT distribution compared to the inactive state, is captured by both the AA and CG models. The shortcut that transfers energy between TM3 and TM4 is TYR207–TRP246, which is intact in both the active and inactive states. The polar contact between this amino acid pair was found in previous work¹⁷ to be more flexible in the active state than in the inactive state, so that the rate of energy transfer across the contact is lower in the active state. The FPT distributions plotted in Fig. 4(b) clearly illustrate how the change in dynamics upon activation of a GPCR, in this case giving rise to a polar contact with greater fluctuations in its length in the active state, alters the energy flow between helices.

Rates of energy transfer via non-covalent contacts serve as sensitive probes of protein dynamics, and entropy associated with the dynamics, as confirmed by theoretical and computational studies of a variety of proteins.^44–47,60 Since both structural fluctuations and rates of vibrational energy transfer can be measured by 2-dimensional infrared (2D-IR) spectroscopy, connections between them can be analyzed by this experimental technique.⁹⁴

Fig. 4(c) and (d) represent cases where there is little shift in the short-time peak between the inactive and active states, consistent in both the AA and CG results. For the Val218–Ile280 source-sink combination plotted in Fig. 4(c), the short-time peak is mainly seen as a shoulder of the long-time peak at shorter times. Apparently, the noncovalent interaction that provides a shortcut for energy flow between TM3 and TM5 is not sufficiently strong in either active or inactive states to give rise to a more prominent feature at shorter times. The source-sink combination Leu153–Ala210 is plotted in Fig. 4(d). The shortcut is due to a hydrogen-bonded contact between Ser158–Asn212 of TM3 and TM4, which is found in both the active and inactive states. We observe only a small shift in the short-time peak of the FPT distribution when the state changes from active to inactive, for both the AA and CG models. However, the shift, while small, is in the opposite direction for the two models. There is a small shift to shorter times upon activation based on the AA data, whereas there is a small shift to longer times upon activation for the CG model.

The FPT distributions plotted in Fig. 4 depend mainly on a single noncovalent contact between sources and sinks on or near different TM helices. There may also be several noncovalent contacts between a source and sink. We plot in Fig. 5(a) the FPT distribution for the Gly50–His71 combination, which is connected by a beta turn with several hydrogen bonds providing alternative, competing pathways. The beta turn and the location of the source and sink are shown in Fig. 5(b) and highlighted in Fig. 1.

Fig. 5 (a) FPT distribution plotted as P(ln t) vs. ln t, for source-sink combination Gly50–His71 of the beta turn for active (blue) and inactive (red) states using AA (solid curve) and CG (dashed) models. (b) The section of the extracellular domain that consists of the beta turn of GLP-1R and nearby loop. (c) FPT distribution Gly50–His71 at the original temperature of 300 K using the CG model and at temperatures that are lower by the factor 0.5 (green), 0.4 (orange) and 0.25 (blue).

For the beta turn, we see that the short-time peak in the FPT distribution is more prominent than the peak at longer times. This result is consistent with the existence of multiple competing pathways for energy to flow between the source and sink. The position of the short-time peak is about the same for the active and inactive states, and for the AA and CG models. However, the short-time peak for the active state using the CG model arises from two separate, overlapping peaks. This superposition can be seen more clearly by lowering the temperature, as shown in Fig. 5(c). We note that the other short-time peaks may also mask features that overlap, but we could not separate them by lowering the temperature, suggesting that if they exist, the corresponding time scales are closer.

Mean first passage times (MFPTs) are calculated from the FPT distribution, p(t), as the integral of tp(t) to specific observation times, t_obs, following ref. 91. This approach enables us to account for the time scale accessible in a given experiment in our analysis of energy transfer. One important limit is the time for excess vibrational energy in the protein to relax into the environment, which is roughly 5 or 10 ps.^{22,23,95–97} We consider energy dissipation times of 5 ps to calculate values of the MFPT here, as we did in a previous study of MFPTs for source-sink combinations in β₂AR.¹⁹

MFPTs are generally expected to increase with the distance between residues in a helix. We therefore plot MFPT as a function of the center of mass distance between source and sink residues, as shown in Fig. 6. Data for both the active and the inactive states are combined in Fig. 6, as we do not expect systematic differences for energy transfer along any of the helices that depends on the state of the system, as concluded in the previous study of β₂AR.¹⁹ We plot source-sink combinations obtained using (a) the AA model and (b) the CG model.

Fig. 6 Mean first passage time (MFPT) at observation times of 5 ps vs. center of mass distance between the source-sink combinations on the same helix. Results using (a) AA and (b) CG models. MFPTs for combinations between which there is no glycine or proline are in red and via glycine and/or proline are in gray. Linear fits to the gray (red) appear as gray (red) lines. The inset shows linear fits for distances > 20Å.

The source-sink combinations in Fig. 6(a) and (b) are grouped in two sets. If there is a glycine or proline along the sequence of amino acids between the combinations, we plot the MFPT in gray, and if there is no glycine or proline we plot the MFPT in red. In Fig. 6(a), containing data obtained with the AA model, the linear fit to the gray data appears as a gray line with a slope of 0.081 ps Å⁻¹. The linear fit to the red data appears as a red line with a slope of 0.073 ps Å⁻¹. We see that the speed of energy transfer along the helix, estimated here as the ratio of the distance between the centers of mass to the fitted MFPT, is faster, 13.7 Å ps⁻¹, when no proline or glycine is present than if one or both are present, 12.3 Å ps⁻¹. The impact of the longer MFPT when glycine or proline is present is magnified at greater distances, as shown in the inset.

Turning to the MFPT results obtained with the CG model, plotted in Fig. 6(b), we see that the linear fit to the data for source-sink combinations where there is proline or glycine in between, plotted in gray, yields a slope of 0.087 ps Å⁻¹, and the linear fit to the red data, representing source-sink combinations without proline or glycine, yields a slope of 0.074 ps Å⁻¹. These slopes are similar to those obtained using the AA model. For the speed of energy transfer along the helix, we again find it is faster, 13.5 Å ps⁻¹, when no proline or glycine is present than if one or both are present, 11.5 Å ps⁻¹. Both these speeds are close to the corresponding values found using the AA model. As for the results based on the AA data, we find the impact of the presence of proline or glycine on the MFPT to be more noticeable at greater distances between the source and sink, as shown in the inset.

We now compare the MFPT values as a function of distance along helices of GLP-1R and β₂AR. For β₂AR, a rate of 12.7 Å ps⁻¹ was obtained when no glycine or proline is present and is otherwise 10.6 Å ps⁻¹.¹⁹ The corresponding values for GLP-1R are 13.7 Å ps⁻¹ when no proline or glycine is present and 12.3 Å ps⁻¹ otherwise. These values are quite similar, and the difference between the rate when no proline or glycine is present and when at least one of them is present is essentially the same. We note that all these rates are below the speed of sound in a protein, which has been computed⁴⁰ and measured¹⁰ to be about 20 Å ps⁻¹.

Conclusions

Vibrational energy dynamics of GLP-1R and properties of the protein that control this process were studied computationally using all atom (AA) and coarse grained (CG) models. In a previous analysis of GLP-1R,¹⁷ information about energy transport in this protein in active and inactive states was obtained using an AA model. In this report, we computed disconnectivity graphs representing the vibrational energy landscape using both AA and CG models and compared vibrational energy dynamics for them. We recently analyzed the corresponding landscape for another GPCR, β₂AR,¹⁹ revealing the impact of certain residues and some structural properties on energy transport that had not been seen before. Specifically, prolines and glycines were found to impede energy flow along the backbone and redirect energy transfer between helices via nearby noncovalent contacts. We find the same effect in GLP-1R, even though the proline and glycine residues are generally located in different regions of the protein compared to β₂AR. Hence, our conclusions about the role of proline and glycine in vibrational energy transfer appear to have some generality. By impeding energy flow along the backbone, prolines and glycines in GPCRs, in combination with nearby noncovalent contacts such as hydrogen bonds, facilitate energy flow between helices, helping to create channels of dynamically coupled residues that can be quite remote in sequence. We suggest that this observation may have important consequences in understanding how these proteins perform their cellular functions.

Both energy transport along the backbone and across noncovalent contacts contribute to the energy transfer dynamics of GLP-1R. When the protein changes its structure due to activation by a ligand, some of the noncovalent contacts in the inactive state are broken, while others form. Even if the noncovalent contacts remain intact upon activation, changes in fluctuations in the length of the contact impact the rate of energy transfer across them. FPT distributions for the energy transfer time were computed for energy source-energy sink combinations of residues in different TM helices, where specific noncovalent contacts facilitate energy flow from the source to the sink. Differences in the FPT distributions for the active and inactive states were found when the noncovalent contact providing a pathway was either broken in the active or inactive state, or there were differences in the dynamics of the contact in the two states, in both cases affecting the rate of energy transfer.

The FPT distributions for both inactive and active states were found to be similar using data obtained from the AA and CG models. For energy transport along TM helices, we computed MFPTs between numerous source-sink combinations to examine the impact of prolines or glycines on energy transport along the helix. For GLP-1R, we found approximately the same speeds, as determined by the ratio of the distance between the residues to the MFPT, when using the AA model and the CG model. When no glycine or proline is present between the energy source and sink, the speed is about 2 Å ps⁻¹ faster than when at least one such residue is present, a consistent result for the AA and CG models, and for β₂AR in previous work using an AA model.¹⁹

The consistent description of energy flow found for GLP-1R, using the AA and CG models, is encouraging for the use of the more computationally efficient CG models in studying energy flow in proteins. Coarse graining has obvious advantages in the study of larger systems, for instance, if the GPCR is complexed with the G-protein, as well as other protein complexes. The present results suggest that appropriate coarse-graining can preserve the dynamics of interest. It will be particularly interesting to explore the role of proline and glycine in structure–function relations for related systems. We suspect that the evolutionary role of these residues, which is well known in terms of structure, may also be exploited in important functional properties, including allosteric interactions.

Conflicts of interest

There are no conflicts to declare.

Data availability

Supplementary information is available. See DOI: https://doi.org/10.1039/d5cp02918j.

The rate constant matrices used to generate disconnectivity graphs and FTP distributions are available on GitHub, https://github.com/PathickHalder/GLP1R.git. MFPT data plotted in Fig. 6 are listed in the SI.

Acknowledgements

Support from NSF grant CHE-2245240 is gratefully acknowledged.

References

1. J. P. Klinman, S. M. Miller and N. G. J. Richards, A Foundational Shift in Models for Enzyme Function, J. Am. Chem. Soc., 2025, 147, 14884–14904.
2. R. Kozlowski, J. Zhao and R. B. Dyer, Acceleration of catalysis in dihydrofolate reductase by transient, site-specific photothermal excitation, Proc. Natl. Acad. Sci. U. S. A., 2021, 118, e2014592118.
3. J. P. T. Zaragoza, A. R. Offenbacher, S. Hu and J. P. Klinman, Temporal and spatial resolution of distal protein motions that activate hydrogen tunneling in soybean lipoxygenase, Proc. Natl. Acad. Sci. U. S. A., 2023, 120, e2211630120.
4. M. Post, B. Lickert, G. Diez, S. Wolf and G. Stock, Cooperative protein allosteric transition mediated by a fluctuating transmission network, J. Mol. Biol., 2022, 434, 167679.
5. G. Stock and P. Hamm, A Nonequilibrium Approach to Allosteric Communication, Philos. Trans. R. Soc., B, 2018, 373, 20170187.
6. Q. Cui and M. Karplus, Allostery and cooperativity revisited, Protein Sci., 2008, 17, 1295–1307.
7. D. M. Leitner, C. Hyeon and K. M. Reid, Water-mediated biomolecular dynamics and allostery, J. Chem. Phys., 2020, 152, 240901.
8. D. M. Leitner and J. E. Straub, Proteins: Energy, Heat and Signal Flow, CRC Press, Taylor & Francis Group, Boca Raton, FL, 2010.
9. D. M. Leitner, Energy flow in proteins, Annu. Rev. Phys. Chem., 2008, 59, 233–259.
10. M. Levatino, G. Schiro, H. T. Lemke, G. Cottone, J. M. Glownia, D. Zhu, M. Chollet, H. Ihee, A. Cupane and M. Cammarata, Ultrafast myoglobin structural dynamics observed with an X-ray free-electron laser, Nat. Commun., 2015, 6, 6772.
11. T. R. M. Barends, L. Foucar, A. Ardevol, K. Nass, A. Aquila, S. Botha, R. B. Doak, K. Falahati, E. Hartmann, M. Hilpert, M. Heinz, M. C. Hoffmann, J. Köflinger, J. E. Koglin, G. Kovacsova, M. Liang, D. Milanthianaki, G. J. Williams, I. Burghardt, G. Hummer, S. Boutet and I. Schlichting, Direct observation of ultrafast collective motions in CO myoglobin upon ligand dissociation, Science, 2015, 350, 445–450.
12. L. U. L. Brinkmann and J. S. Hub, Ultrafast anisotropic protein quake propagation after CO photodissociation in myoglobin, Proc. Natl. Acad. Sci. U. S. A., 2016, 113, 10565–10570.
13. R. W. Schoenlein, L. A. Peteanu, R. A. Mathies and C. V. Shank, The first step in vision: femtosecond isomerization of rhodopsin, Science, 1994, 254, 412–415.
14. J. E. Kim, M. J. Tauber and R. A. Mathies, Analysis of the mode-specific excited-state energy distribution and wavelength-dependent photoreaction quantum yield in rhodopsin, Biophys. J., 2003, 84, 846–867.
15. P. Kukura, D. W. McCamant, S. Yoon, D. B. Wandschneider and R. A. Mathies, Structural observation of the primary isomerization in vision with femtosecond-stimulated Raman, Science, 2005, 310, 1006–1009.
16. H. Poudel and D. M. Leitner, Activation-Induced Reorganization of Energy Transport Networks in the B2 Adrenergic Receptor, J. Phys. Chem. B, 2021, 125, 6522–6531.
17. H. Poudel and D. M. Leitner, Energy transport in class B GPCRs: Role of protein-water dynamics and activation, J. Phys. Chem. B, 2022, 42, 8362–8373.
18. H. Poudel and D. M. Leitner, Locating dynamic contributions to allostery via determining rates of vibrational energy transfer, J. Chem. Phys., 2023, 158, 015101.
19. H. Poudel, D. J. Wales and D. M. Leitner, Vibrational energy landscapes and energy flow in GPCRs, J. Phys. Chem. B, 2024, 128, 7568–7576.
20. W. I. Weis and B. K. Kobilka, The molecular basis of G-protein coupled receptor activation, Annu. Rev. Biochem., 2018, 87, 897–919.
21. M. B. Prigozhin and M. Gruebele, Microsecond folding experiments and simulations: A match is made, Phys. Chem. Chem. Phys., 2013, 15, 3372–3388.
22. S. Buchenberg, D. M. Leitner and G. Stock, Scaling rules for vibrational energy transport in proteins, J. Phys. Chem. Lett., 2016, 7, 25–30.
23. D. M. Leitner, S. Buchenberg, P. Brettel and G. Stock, Vibrational energy flow in the villin headpiece subdomain: Master equation simulations, J. Chem. Phys., 2015, 142, 075101.
24. A. Bortolato, A. S. Doré, K. Hollenstein, B. G. Tehan, J. S. Mason and F. H. Marshall, Structure of Class B GPCRs: New Horizons for Drug Discovery, Br. J. Pharmacol., 2014, 171, 3132–3145.
25. C. de Graaf, D. Donnelly, D. Wootten, J. Lau, P. M. Sexton, L. J. Miller, J. M. Ahn, J. Liao, M. M. Fletcher, D. Yang, A. J. H. Brown, C. Zhou, J. Deng and M.-W. Wang, Glucagon-like Peptide-1 and Its Class B G Protein-Coupled Receptors: A Long March to Therapeutic Successes, Pharmacol. Rev., 2016, 68, 954–1013.
26. M. A. Nauck and J. J. Meier, Are All GLP-1 Agonists Equal in the Treatment of Type 2 Diabetes?, Eur. J. Endocrinol., 2019, 181, R211–R234.
27. F. Malik and Z. Li, Non-Peptide Agonists and Positive Allosteric Modulators of Glucagon-like Peptide-1 Receptors: Alternative Approaches for Treatment of Type 2 Diabetes, Br. J. Pharmacol., 2021, 179, 511–525.
28. Y. Zhang and J. E. Straub, Directed Energy Funneling in Proteins: From Structure to Function, in Proteins: Energy, Heat and Signal Flow, ed. D. M. Leitner and J. E. Straub, Taylor & Francis Group, CRC Press, Boca Raton, 2009, pp. 199–228.
29. E. H. G. Backus, P. H. Nguyen, V. Botan, A. Moretto, M. Crisma, C. Toniolo, O. Zerbe, G. Stock and P. Hamm, Structural flexibility of a helical peptide regulates vibrational energy transport properties, J. Phys. Chem. B, 2008, 112, 15487–15492.
30. D. M. Leitner and T. Yamato, Mapping energy transport networks in proteins, in Rev. Comp. Chem., ed. A. L. Parrill and K. B. Lipkowitz, John Wiley & Sons, Inc., 2018, vol. 31, pp. 63–114.
31. D. M. Leitner and T. Yamato, Recent developments in the computational study of protein structural and vibrational energy dynamics, Biophys. Rev., 2020, 12, 317–322.
32. K. Ota and T. Yamato, Energy exchange network model demonstrates protein allosteric transition: An application to an oxygen sensor protein, J. Phys. Chem. B, 2019, 123, 768–775.
33. T. Wang, T. Yamato and W. Sugiura, Thermal energy transport through nonbonded native contacts in protein, J. Phys. Chem. B, 2024, 128, 8641–8650.
34. E. Deniz, L. Valiño-Borau, J. G. Löffler, K. B. Eber, A. Gulzar, S. Wolf, P. M. Durkin, R. Kaml, N. Budisa, G. Stock and J. Bredenbeck, Through bonds or contacts? Mapping protein vibrational energy transfer using non-canonical amino acids, Nat. Commun., 2021, 12, 3284.
35. N. Fujii, M. Mizuno, H. Ishikawa and Y. Mizutani, Observing vibrational energy flow in a protein with the spatial resolution of a single amino acid residue, J. Phys. Chem. Lett., 2014, 5, 3269–3273.
36. M. Kondoh, M. Mizuno and Y. Mizutani, Importance of atomic contacts in vibrational energy flow in proteins, J. Phys. Chem. Lett., 2016, 7, 1950–1954.
37. M. Mizuno and Y. Mizutani, Role of atomic contacts in vibrational energy transfer in myoglobin, Biophys. Rev., 2020, 12, 511–518.
38. S. Yamashita, M. Mizuno, D. P. Tran, H. Dokainish, A. Kitao and Y. Mizutani, Vibrational Energy Transfer from Heme through Atomic Contacts in Proteins, J. Phys. Chem. B, 2018, 122, 5877–5884.
39. A. Buchmann, L. Caminiti, L. L. Stiewe, S. Jung, D. M. Leitner and M. Havenith, Internal vibrational energy redistribution precedes energy dissipation into the solvent upon photoexcitation of heme proteins, Phys. Chem. Chem. Phys., 2025, 27, 9470–9477.
40. X. Yu and D. M. Leitner, Vibrational energy transfer and heat conduction in a protein, J. Phys. Chem. B, 2003, 107, 1698–1707.
41. X. Yu and D. M. Leitner, Anomalous diffusion of vibrational energy in proteins, J. Chem. Phys., 2003, 119, 12673–12679.
42. M. B. Enright and D. M. Leitner, Mass fractal dimension and the compactness of proteins, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2005, 71, 011912.
43. J. K. Agbo, Y. Xu, P. Zhang, J. E. Straub and D. M. Leitner, Vibrational energy flow across heme-cytochrome c and cytochrome c-water interfaces, Theor. Chem. Acc., 2014, 133, 1504.
44. K. M. Reid, T. Yamato and D. M. Leitner, Scaling of rates of vibrational energy transfer in proteins with equilibrium dynamics and entropy, J. Phys. Chem. B, 2018, 122, 9331–9339.
45. K. M. Reid, T. Yamato and D. M. Leitner, Variation of Energy Transfer Rates across Protein–Water Contacts with Equilibrium Structural Fluctuations of a Homodimeric Hemoglobin, J. Phys. Chem. B, 2020, 124, 1148–1159.
46. K. M. Reid, X. Yu and D. M. Leitner, Change in vibrational entropy with change in protein volume estimated with mode Grüneisen parameters, J. Chem. Phys., 2021, 154, 055012.
47. H. Poudel, K. M. Reid, T. Yamato and D. M. Leitner, Energy transfer across nonpolar and polar contacts in proteins: role of contact fluctuations, J. Phys. Chem. B, 2020, 124, 9852–9861.
48. H. Poudel and D. M. Leitner, What measurements of energy transfer across non-covalent contacts in proteins reveal about the energy landscape, J. Innovative Mater. Extreme Cond., 2024, 5, 45–56.
49. S. Gao, X. T. Wu, W. Zhang, T. Richardson, S. L. Barrow, C. A. Thompson-Kucera, A. T. Iarone and J. P. Klinman, Temporal resolution of activity-related solvation dynamics in the TIM barrel enzyme murine adenosine deaminase, ACS Catal., 2024, 14, 4554–4567.
50. O. M. Becker and M. Karplus, The topology of multidimensional potential energy surfaces: Theory and application to peptide structure and kinetics, J. Chem. Phys., 1997, 106, 1495–1517.
51. D. J. Wales, M. A. Miller and T. R. Walsh, Archetypal energy landscpaes, Nature, 1999, 394, 758–760.
52. J. J. Joseph, K. Röder, D. Chakraborty, R. G. Mantell and D. J. Wales, Exploring biomolecular energy landscapes, Chem. Commun., 2017, 53, 6974–6988.
53. H. Eyring, The activated complex in chemical reactions, J. Chem. Phys., 1935, 3, 107–115.
54. M. G. Evans and M. Polanyi, Some applications of the transition state method to the calculation of reaction velocities, especially in solution, Trans. Faraday Soc., 1935, 31, 875–894.
55. C. N. Pace and J. M. Scholtz, A helix propensity scale based on experimental studies of peptides and proteins, Biophys. J., 1998, 75, 422–427.
56. D. M. Rosenbaum, S. G. Rasmussen and B. K. Kobilka, The structure and function of G-protein-coupled receptors, Nature, 2009, 459, 356–363.
57. E. B. VanArnam, H. A. Lester and D. A. Dougherty, Dissecting the Functions of Conserved Prolines within Transmembrane Helices of the D2 Dopamine Receptor, ACS Chem. Biol., 2011, 6, 1063–1068.
58. P. Mazna, L. Grycova, A. Balik, H. Zemkova, E. Friedlova, V. Obsilova, T. Obsit and J. Teisinger, The Role of Proline Residues in the Structure and Function of Human MT2 Melatonin Receptor, J. Pineal Res., 2008, 45, 361–372.
59. X. Gan, Z. Ma, N. Deng, J. Wang, J. Ding and L. Li, Involvement of the C-Terminal Proline-Rich Motif of G Protein-Coupled Receptor Kinases in Recognition of Activated Rhodopsin, J. Biol. Chem., 2004, 279, 49741–49746.
60. D. M. Leitner, H. D. Pandey and K. M. Reid, Energy Transport across Interfaces in Biomolecular Systems, J. Phys. Chem. B, 2019, 123, 9507–9524.
61. J. S. Burg, J. R. Ingram, A. J. Venkatakrishnan, K. M. Jude, A. Dukkipati, E. N. Feinberg, A. Angelini, D. Waghray, R. O. Dror, H. L. Ploegh and K. C. Garcia, Structural basis for chemokine recognition and activation of a viral G protein-coupled receptor, Science, 2015, 347, 1113–1117.
62. B. Trzaskowski, D. Latek, S. Yuan, U. Ghoshdastider, A. Debinski and S. Filipak, Action of molecular switches in GPCRs – Thoertical and experimental studies, Curr. Med. Chem., 2012, 19, 1090–1109.
63. N. R. Latorraca, A. J. Venkatakrishnan and R. O. Dror, GPCR dynamics: Structures in motion, Chem. Rev., 2017, 117, 139–155.
64. A. J. Venkatakrishnan, A. K. Ma, R. Fonseca, N. R. Latorraca, B. Kelly, R. M. Betz, C. Asawa, B. K. Kobilka and R. O. Dror, Diverse GPCRs Exhibit Conserved Water Networks for Stabilization and Activation, Proc. Natl. Acad. Sci. U. S. A., 2019, 116, 3288–3293.
65. Y. Lee, S. Choi and C. Hyeon, Mapping the intramolecular signal transduction of G-protein coupled receptors, Proteins: Struct., Funct., Bioinf., 2014, 82, 727–743.
66. Y. Lee, S. Kim, S. Choi and C. Hyeon, Ultraslow water-mediated transmembrane interactions regulate the activation of A2A adenosine receptor, Biophys. J., 2016, 111, 1180–1191.
67. S. K. Huang, A. Pandey, D. P. Tran, N. L. Villanueva, A. Kitao, R. K. Sunahara, A. Sljoka and R. S. Prosser, Delineating the conformational landscape of the adenosine A2A receptor during G protein coupling, Cell, 2021, 184, 1884–1894.
68. S. Hashem, A. Dougha and P. Tuffery, Ligand-induced biased activation of GPCRs: Recent advances and new directions from in silico approaches, Molecules, 2025, 30, 1047.
69. A. J. Venkatakrishnan, X. Deupi, G. Lebon, C. G. Tate, G. F. Schertler and M. M. Babu, Molecular signatures of G protein coupled receptors, Nature, 2013, 494, 185–194.
70. K. Hollenstein, C. deGraaf, A. Bortolato, M.-W. Wang, F. H. Marshall and R. C. Stevens, Insights into the structure of class B GPCRs, Trends Pharmacol. Sci., 2014, 35, 12–22.
71. F. Wu, L. Yang, K. Hang, M. Laursen, L. Wu, G. W. Han, Q. Ren, N. K. Roed, G. Lin, M. A. Hansen, H. Jiang, M.-W. Wang, S. Reedtz-Runge, G. Song and R. C. Stevens, Full-length human GLP-1 receptor structure without orthosteric ligands, Nat. Commun., 2020, 11, 1272.
72. B. Krumm and B. L. Roth, A structural understanding of Class B GPCR selectivity and activation revealed, Structure, 2020, 28, 277–279.
73. Y.-L. Liang, M. Khoshouei, M. Radjainia, Y. Zhang, A. Glukhova, J. Tarrasch, D. M. Thal, S. G. B. Furness, G. Christopoulos, T. Coudrat, R. Danev, W. Baumeister, L. J. Miller, A. Christopoulos, B. K. Kobilka, D. Wootten, G. Skiniotis and P. M. Sexton, Phase-plate cryo-EM structure of a class B GPCR–G-protein complex, Nature, 2017, 546, 118–123.
74. Y. Zhang, B. Sun, D. Feng, H. Hu, M. Chu, Q. Qu, J. T. Tarrasch, S. Li, T. S. Kobilka, B. K. Kobilka and G. Skiniotis, Cryo-EM structure of the activated GLP-1 receptor in complex with a G protein, Nature, 2017, 546, 248–253.
75. G. Song, D. Yang, Y. Wang, C. DeGraaf, Q. Zhou, S. Jiang, K. Liu, X. Cai, A. Dai, G. Lin, D. Liu, F. Wu, Y. Wu, S. Zhao, L. Ye, G. W. Han, J. Lau, B. Wu, M. A. Hanson, Z. J. Liu, M. W. Wang and R. C. Stevens, Human GLP-1 Receptor Transmembrane Domain Structure in Complex with Allosteric Modulators, Nature, 2017, 546, 312–315.
76. H. Zhang, A. Qiao, D. Yang, L. Yang, A. Dai, C. DeGraaf, S. Reedtz-Runge, V. Dharmarajan, H. Zhang, G. W. Han, T. D. Grant, R. G. Sierra, U. Weierstall, G. Nelson, W. Liu, Y. Wu, L. Ma, X. Cai, G. Lin, X. Wu, Z. Geng, Y. Dong, G. Song, P. R. Griffin, J. Lau, V. Cherezov, H. Yang, M. A. Hanson, R. C. Stevens, Q. Zhao, H. Jiang, M. W. Wang and B. Wu, Structure of the Full-Length Glucagon Class B G-Protein-Coupled Receptor, Nature, 2017, 546, 259–264.
77. P. Zhao, Y.-L. Liang, M. J. Belousoff, G. Deganutti, M. M. Fletcher, F. S. Willard, M. G. Bell, M. E. Christe, K. W. Sloop, A. Inoue, T. T. Truong, L. Clydesdale, S. G. B. Furness, A. Christopoulos, M.-W. Wang, L. J. Miller, C. A. Reynolds, R. Danev, P. M. Sexton and D. Wootten, Activation of the GLP-1 Receptor by a Non-Peptidic Agonist, Nature, 2020, 577, 432–436.
78. H. Ma, W. Huang, X. Wang, L. Zhao, Y. Jiang, F. Liu, W. Guo, X. Sun, W. Zhong, D. Yuan and H. E. Xu, Structural insights into the activation of GLP-1R by a small molecule agonist, Cell Res., 2020, 30, 1140–1142.
79. A. Fiser and A. Sali, ModLoop: Automated Modeling of Loops in Protein Structures, Bioinformatics, 2003, 19, 2500–2501.
80. P. Mark and L. Nilsson, Structure and Dynamics of the TIP3P, SPC, and SPC/E Water Models at 298 K, J. Phys. Chem. A, 2001, 105, 9954–9960.
81. C. J. Dickson, B. D. Madej, Å. A. Skjevik, R. M. Betz, K. Teigen, I. R. Gould and R. C. Walker, Lipid14: The Amber Lipid Force Field, J. Chem. Theory Comput., 2014, 10, 865–879.
82. J. A. Maier, C. Martinez, K. Kasavajhala, L. Wickstrom, K. E. Hauser and C. Simmerling, Ff14SB: Improving the Accuracy of Protein Side Chain and Backbone Parameters from Ff99SB, J. Chem. Theory Comput., 2015, 11, 3696–3713.
83. H. J. C. Berendsen, J. P. M. Postma, W. F. van Gunsteren, A. DiNola and J. R. Haak, Molecular Dynamics with Coupling to an External Bath, J. Chem. Phys., 1984, 81, 3684–3690.
84. T. Ishikura, Y. Iwata, T. Hatano and T. Yamato, Energy Exchange Network of Inter-Residue Interactions Within a Thermally Fluctuating Protein: A Computational Study, J. Comput. Chem., 2015, 36, 1709–1718.
85. M. R. Machado, E. E. Barrera, F. Klein, M. Sóñora, S. Silva and S. Pantano, The SIRAH 2.0 Force Field: Altius, Fortius, Citius, J. Chem. Theory Comput., 2019, 15, 2719–2733.
86. L. Darré, M. R. Machado, A. F. Brandner, H. C. González, S. Ferreira and S. Pantano, SIRAH: A Structurally Unbiased Coarse-Grained Force Field for Proteins with Aqueous Solvation and Long-Range Electrostatics, J. Chem. Theory Comput., 2015, 11, 723–739.
87. L. Monticelli, S. K. Kandasamy, X. Periole, R. G. Larson, D. P. Tieleman and S.-J. Marrink, The MARTINI Coarse-Grained Force-Field: Extension to Proteins, J. Chem. Theory Comput., 2008, 4, 819–834.
88. P. C. T. Souza, R. Alessandri, J. Barnoud, S. Thallmair, I. Faustino, F. Grünewald, I. Patmanidis, H. Abdizadeh, B. M. H. Bruininks, T. A. Wassenaar, P. C. Kroon, J. Melcr, V. Nieto, V. Corradi, H. M. Khan, J. Domański, M. Javanainen, H. Martinez-Seara, N. Reuter, R. B. Best, I. Vattulainen, L. Monticelli, X. Periole, D. P. Tieleman, A. H. de Vries and S. J. Marrink, Martini 3: A General Purpose Force Field for Coarse-Grained Molecular Dynamics, Nat. Methods, 2021, 18, 382–388.
89. E. J. Woods and D. J. Wales, Analysis and interpretation of first passage time distributions featuring rare events, Phys. Chem. Chem. Phys., 2024, 26, 1640–1657.
90. E. J. Woods, D. Kannan, D. J. Sharpe, T. D. Swinburne and D. J. Wales, Analysing ill-conditioned Markov chains, Philos. Trans. R. Soc., A, 2023, 381, 20220245.
91. D. J. Wales, Dynamical signatures of multifunnel landscapes, J. Phys. Chem. Lett., 2022, 13, 6349–6358.
92. E. W. Dijkstra, A note on two problems in connection with graphs, Numer. Math., 1959, 1, 269–271.
93. D. J. Wales, Discrete path sampling, Mol. Phys., 2002, 100, 3285–3306.
94. R. Maia and C. R. Baiz, Membrane Interactions Drive Hydrogen-Bond Dynamics in the pHLow Insertion Peptide (pHLIP), J. Phys. Chem. Lett., 2025, 16, 9583–9589.
95. P. H. Nguyen, S. M. Park and G. Stock, Nonequilibrium molecular dynamics simulation of the energy transfer through a peptide helix, J. Chem. Phys., 2010, 132, 025102.
96. P. H. Nguyen, P. Hamm and G. Stock, Nonequilibrium molecular dynamics simulation of photoinduced energy flow in peptides: Theory meets experiment, in Proteins: Energy, Heat and Signal Flow, ed. D. M. Leitner and J. E. Straub, CRC Press, Taylor & Francis Group, Boca Raton, 2010, pp. 149–168.
97. V. Botan, E. H. G. Backus, R. Pfister, A. Moretto, M. Crisma, C. Toniolo, P. H. Nguyen, G. Stock and P. Hamm, Energy transport in peptide helices, Proc. Natl. Acad. Sci. U. S. A., 2007, 104, 12749–12754.

---

Footnote

† These authors contributed equally.

---