<figure><img src="image_1.png" /></figure>

HAL
open science

# Molecular Dynamics Simulations of the Allosteric Modulation of the Adenosine A2a Receptor by a Mini-G Protein

Pedro Renault, Maxime Louet, Jacky Marie, Gilles Labesse, Nicolas Floquet

► To cite this version:

Pedro Renault, Maxime Louet, Jacky Marie, Gilles Labesse, Nicolas Floquet. Molecular Dynamics Simulations of the Allosteric Modulation of the Adenosine A2a Receptor by a Mini-G Protein. Scientific Reports, 2019, 9 (1), 10.1038/s41598-019-41980-x. http://02322127

HAL Id: hal-02322127

https://hal.science/hal-02322127

Submitted on 2 Mar 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
SCIENTIFIC REPORTS

<figure><img src="image_3.png" /></figure>

OPEN

# Molecular Dynamics Simulations of the Allosteric Modulation of the Adenosine A2A Receptor by a Mini-G Protein

Pedro Renault ID<fnref n="1" />,², Maxime Louet<fnref n="1" />, Jacky Marie<fnref n="1" />, Gilles Labesse² & Nicolas Floquet ID³

Received: 31 December 2018

Accepted: 19 March 2019

Published online: 02 April 2019

Through their coupling to G proteins, G Protein-Coupled Receptors (GPCRs) trigger cellular responses to various signals. Some recent experiments have interestingly demonstrated that the G protein can also act on the receptor by favoring a closed conformation of its orthosteric site, even in the absence of a bound agonist. In this work, we explored such an allosteric modulation by performing extensive molecular dynamics simulations on the adenosine A2 receptor (A2AR) coupled to the Mini-Gs protein. In the presence of the Mini-Gs, we confirmed a restriction of the receptor's agonist binding site that can be explained by a modulation of the intrinsic network of contacts of the receptor. Of interest, we observed similar effects with the C-terminal helix of the Mini-Gs, showing that the observed effect on the binding pocket results from direct local contacts with the bound protein partner that cause a rewiring of the whole receptor's interaction network.

G protein-coupled receptors (GPCRs) constitute the largest family of human membrane receptors and are involved in numerous physiological processes. Among other partners, they can couple to G proteins to trigger signals that result in complex cellular responses<fnref n="1" />. It has been recognized that GPCRs are dynamical molecules that exist as an ensemble of conformations, whose equilibrium can be shifted by their interactions with other molecules²,³. Therefore, understanding GPCRs activation and designing drugs to modulate their function require a better comprehension of their dynamical behaviour.

The dynamics of GPCRs and the effect of ligands upon their conformational equilibria have been extensively investigated⁴⁻⁸. Computational approaches, including molecular dynamics simulations (conventional as well as enhanced sampling schemes), have contributed to provide a clearer picture of (1) the activation process and (2) the response of GPCRs to different types of ligands (agonists, antagonists or inverse agonists)⁹⁻¹³. The simulations have benefited enormously from the structural data accumulated in recent years, that culminated in the structural resolution of complexes formed by the GPCR bound to different heterotrimeric G proteins¹⁴⁻¹⁷. These structures revealed details of the interaction between receptor and G protein and established a basis for the comprehension of signal transmission from the agonist-bound receptor to the G protein. Accordingly, computational methods have been used to study the allosteric connections between the agonist binding site of the receptor and the nucleotide binding site of the G protein, usually focusing on the modulation of the G protein by the receptor¹⁸,¹⁹. However, in a study recently published by Devree et al., experiments with the β2aR have shown that the G protein can also allosterically modulate the receptor even in the absence of agonist²⁰. The authors captured the receptor coupled to the G protein but unbound to the agonist and demonstrated that the G protein hinders the association of ligands to the orthosteric site. Moreover, it has also been shown that the G protein slows the dissociation of ligands once they are bound to the site, resulting in the enhanced agonist affinity observed for many GPCRs in the presence of G protein²⁰. These observations led to the hypothesis that the G protein stabilizes a closed conformation of the receptor, with a restricted access to the agonist binding site²⁰. Based on their experiments and on available crystal structures of receptors coupled to G proteins, the authors of that study proposed a mechanism of signal transmission from the G protein to the orthosteric binding site of the receptor and suggested that this could

<footnote n="1">¹Institut des Biomolécules Max Mousseron (IBMM), CNRS UMR5247, Université de Montpellier, ENSCM, 34090 Montpellier, France. ²Centre de Biochimie Structurale, Université de Montpellier, CNRS, INSERM, 34090 Montpellier, France. Correspondence and requests for materials should be addressed to N.F. (email: nicolas.floquet@umontpellier.fr)</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 1.** The A2AR:Mini-Gs:NECA complex. (A) A2AR (blue) bound to Mini-Gs (green). The Cα's of residues of the agonist binding site are shown as purple spheres and the C-terminal helix (α5) of Mini-Gs is in black. The agonist NECA is in orange and is also represented separately. (B) Detail of the agonist site. Side chains are shown as sticks. The extracellular loops were omitted for clarity.</figcaption></figure>

be a general allosteric mechanism among the class A GPCRs²⁰. However, a dynamical view of this modulation is still necessary.

In this work, we proposed to explore the hypothesis of the stabilization of a closed conformation of the receptor by a G protein and to supply a dynamical view of this process.

In particular, we focused on the effect of the G protein upon the unbound receptor and performed molecular dynamics (MD) simulations to study the allosteric signaling from the intracellular partner to the agonist site. We intended to investigate if the G protein was able to stabilize a closed orthosteric site, even in the absence of a bound agonist. Our model system was the adenosine A2 receptor (A2AR), a prototypical class A GPCR, coupled to Mini-Gs (Fig. 1), an engineered protein that binds to A2AR through an interface similar to the one between Gs and β2AR²¹,²².

Although it is reduced to parts of the G protein alpha-subunit, experiments have shown that this mini-protein also causes an increase in agonist affinity, as observed with the full Gs protein²². Moreover, the smaller size of Mini-Gs in comparison to the full Gs represented a reduced computational cost and allowed longer simulation times.

The A2AR is a promising therapeutic target as its modulation by selective agonists and antagonists could help in the treatment of many different diseases including neurological disorders (e.g. Parkinson's disease), cancer and inflammation²³,²⁴. Also, the structure of A2AR has been obtained in different activation states (inactive²⁵⁻³⁰, active-intermediate³¹⁻³³ and active²²,³⁴), what makes it useful for the investigation of receptor activation. Therefore, this is a highly relevant model for the study of dynamics and allostery in GPCRs.

## Results

Starting from the crystal structure of the A2AR:Mini-Gs complex in the presence of the agonist NECA (PDB ID: 5G53; Fig. 1), we built four different systems: (1) A2AR; (2) A2AR:Mini-Gs; (3) A2AR bound to Mini-Gs C-terminal α5 helix only (A2AR:α5), since Mini-Gs contacts A2AR mainly through this helix and many experimental studies have also shown that it is able to mimic the full G protein, both functionally and structurally³⁵,³⁶; (4) A2AR:Mini-Gs:NECA, that served as a reference for the closed state of the orthosteric site. We emphasize that the initial conformation of the receptor was the same in all systems, as they were built by removing either the ligand, Mini-Gs or both. Two independent simulations of 0.5 µs were performed for each system (a total of 4 µs of simulation) to investigate the effect of the different partners on both the conformation of the agonist binding site and the global contact network of the receptor. The stability of all simulations performed in this study was attested by the backbone RMSD of the receptor, displayed in Fig. S1.

**Mini-Gs or the α5 helix restrict the conformational space of the agonist site.** We first looked at the conformational freedom of residues of the agonist binding site using a Principal Component Analysis (PCA) of their cartesian coordinates. Residues containing at least one atom within 4.5 Å of the agonist NECA in the crystal structure of the A2AR:Mini-Gs:NECA complex (PDB ID 5G53) were included as part of the binding
<figure><img src="image_5.png" /></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 2.** Dynamics and geometry of the agonist site. (A) PCA of the cartesian coordinates of the site. The agonist maximally restricted the conformation of the site (red), while the free receptor explored a larger conformational space (light blue); even in the absence of the agonist, Mini-Gs (green) or α5 (black) limited the conformational freedom of the site. The initial conformation of all trajectories (PDB ID 5G53) is represented by a blue square. For reference, an antagonist-bound conformation (PDB ID 5UIG) is also projected (yellow square). (B) Probability distribution of the volume of the agonist site. The y-axis indicates the kernel density estimation of the probability density function. The smallest volumes were observed when the agonist was bound (red). The free receptor (blue) explored conformations with the largest volumes. Even in the absence of agonist, Mini-Gs (green) or α5 (black) were able to maintain a reduced volume of the site, in comparison to the free receptor.</figcaption></figure>

site. These residues are listed in the Supporting Information (Text S1), their C'αs are shown as purple spheres in Fig. 1A and they are shown in detail in Fig. 1B. Figure 2A shows the projection of the conformations obtained for each studied system in a same space, defined by the first two PCs derived from a concatenated trajectory of the receptor, combining all simulated conditions. The initial conformation of A2AR in all simulations, which corresponds to the receptor bound to the agonist NECA and to Mini-Gs, is also projected onto the same PC subspace and represented in the figure (dark blue square). For reference, the projection of an antagonist-bound conformation is also shown (yellow square). In presence of NECA (red) we observed a maximal restriction of the conformation of the agonist site. This was explained by the rigidity of NECA during the simulations, which limited the flexibility of surrounding residues of the receptor: the average RMSD of NECA relative to the crystal structure was of 1.4 Å. In contrast, in the absence of any partner the A2AR agonist site explored a larger conformational space, visiting not only conformations sampled in the presence of intracellular partners or NECA but also others, apparently inaccessible to the bound forms. The most interesting results came from simulations of the receptor coupled only to the Mini-Gs: a restriction of the agonist site was also observed despite the absence of NECA (green). Even though the site was not as restricted as in the presence of the agonist, it could not reach all the conformations explored by the free receptor. Interestingly, a similar result was obtained when A2AR was bound to the α5 helix (black).

**Mini-Gs or the α5 helix stabilize closed conformations of the agonist site.** To investigate if the conformations explored by the distinct systems were associated to different volumes of the agonist site, we measured this volume and calculated its probability distribution in the different conditions. These distributions are illustrated in Fig. 2B. The presence of NECA determined a narrow peak centered at the smallest volume (red). The distribution for the free receptor (blue) was shifted towards higher volumes and was also wider, showing the increased structural variability already indicated by PCA. Coupling to Mini-Gs in the absence of NECA (green) led to an intermediate distribution, where volumes were decreased when compared to the free receptor. Remarkably, a similar result was observed when the α5 helix was bound to A2AR (black). Therefore, Mini-Gs or the α5 helix led not only to a conformational restriction of the agonist site, but also favored more closed conformations of the site in comparison to the uncoupled receptor, even in the absence of agonist.

To test whether the α5 helix could also behave similarly to Mini-Gs in the agonist-bound receptor, we have simulated A2AR bound to NECA and coupled to α5 (A2AR:α5:NECA). In this case, the conformational space of the orthosteric site was as restricted as in the A2AR:Mini-Gs:NECA complex, with volumes comparable to those observed with the full Mini-Gs (Fig. S2). These results suggest that the contacts between A2AR and the α5 helix are an essential ingredient of the allosteric stabilization of the closed forms of the agonist site, with or without the agonist bound to the receptor.

We have also extended the trajectories of A2AR, A2AR:Mini-Gs:NECA and A2AR:Mini-Gs to 2 µs (reaching a total of 12 µs of simulation time for these three systems) and we have found the same pattern of volume distributions (Fig. S3). PCA based on the coordinates of the agonist site performed in these extended trajectories (Fig. S4)
also corroborate the results after 500 ns, shown in Fig. 2. Therefore, a 500 ns simulation time was sufficient for observing the effect of Mini-Gs or α5 on the relaxation and geometry of the agonist site.

To assess the conformational changes responsible for the differences in the volume, we measured the pairwise distances between the residues of the agonist site in all simulated conditions. We evaluated the distances between the Co's and also between the side chains, and determined the Pearson correlation coefficient between every measured distance and the volume. All measured correlations are available as Supplementary Data and their distributions are depicted in histograms in Fig. S5. Backbone and side chain motions played important roles to modify the geometry of the agonist site. Distances between the Co's of residues in TM2 and TM3 were the most strongly correlated with the volume, particularly the distance between residues ALA 63 (2.61) and VAL 84 (3.32), with a Pearson correlation coefficient of 0.69 (Fig. S6; the numbers in parenthesis denote the Ballesteros-Weinstein numbering scheme³⁷). Higher volumes of the orthosteric site were also associated to a larger distance between the backbones of TM3 and TM7; most notably, the distance between the Co's of residues THR 88 (3.36) and HIS 278 (7.43) displayed a correlation coefficient of 0.64 with the volume of the site. HIS 278 was also important because of the conformational changes of its side chain, whose distance to the side chain of TRP 246 (6.48) was significantly associated to variations in the volume of the site (correlation coefficient of 0.63). These distances were indicative of the degree of openness of the site and we studied their distributions in the distinct simulated conditions. The results are illustrated in Fig. 3A-D, with displays of 2D density plots of the distance between the Co's of ALA 63 (2.61) and VAL 84 (3.32) relative to the perimeter of a triangle formed by the Co of THR 88 (3.36) and atoms on the side chains of TRP 246 (6.48) and HIS 278 (7.43) (see the Methods section for details).

In the free receptor (Fig. 3A) it is possible to distinguish a closed state of the site, characterized by a perimeter T88-W246-H278 of approximately 14–16 Å and the distance 63–84 in the range of approximately 6–7 Å; an open state, with the perimeter T88-W246-H278 of ~21–23 Å and the distance A63-V84 of ~13–15 Å; and a set of intermediate conformations, with the triangle formed by T88-W246-H278 assuming a large range of values between ~16 and ~22 Å and the distance A63-V84 of ~10 Å. In contrast, the A2AR:Mini-Gs:NECA system explores only the closed state of the site (Fig. 3B). Interestingly, for A2AR:Mini-Gs and A2AR:α5, the dynamical equilibrium was shifted towards the closed state and although some excursions into intermediate conformations were observed (Fig. 3C,D), the exploration of the open state was abolished in these cases. Therefore, even in the absence of the agonist the intracellular partners contributed to maintain the orthosteric site in the initial, NECA-bound conformation. This indicates a slower relaxation of the agonist site when Mini-Gs or the α5 helix were coupled to the receptor.

Figure 3E shows a superposition of the closed and open states and highlights the structural differences between them. The distance increase between TM2 and TM3 resulted from an outward movement of TM2, that was the single most important contribution to the expansion of the site. Moreover, in the open form of the free receptor TM2 also underwent a rotation (see Fig. S7), that moved ILE 60 (2.58) away from HIS 278 (7.43). By eliminating steric hindrances, the rearrangement of TM2 allowed a relaxation of the side chain of HIS 278 (7.43), that otherwise protruded into the site, contributing to restrict its volume. In agonist-bound forms, HIS 278 contacted NECA and, without the agonist, Mini-Gs or the α5 helix contributed to keep this residue in the NECA-bound conformation by maintaining the closed form of TM2. Our simulations also clearly pointed to the implication of TRP 246 (6.48) in the volume variations of the A2AR agonist site. This is a highly conserved residue known as the “turning switch” and its role in the activation of GPCRs has been extensively discussed in the literature³⁸,³⁹.

Additionally, we also calculated average distance matrices for the entire receptor, based on the distance between Co's, in all simulated conditions. To evaluate deviations from the closed conformation of the receptor, we have obtained difference distance matrices, by subtracting the matrix calculated in A2AR:Mini-Gs:NECA from the matrices of every other condition. The results are illustrated in Fig. S8 and clearly showed the contribution of the backbone motion of TM2 to open the receptor. We also paid special attention to a set of distances between 24 pairs of residues that form a consensus scaffold of non-covalent contacts in GPCRs, independent of the activation state⁴⁰. The frequencies of these contacts (distance of less than 4.5 Å between the closest heavy atoms of the residues) in the distinct systems are available in a Supplementary Table. In agreement with the literature⁴⁰, we observed the overall persistence of these contacts throughout the trajectories.

**Mini-Gs or the α5 helix modify the contact network of the receptor.** To investigate more closely how Mini-Gs was allosterically coupled to the othosteric site and led to the reduction of its volume even when the site was empty, we focused on the main structural differences between the closed state of A2AR:Mini-Gs and the open state of the free receptor (see the Methods section for precisions about the selection of these states). For these two groups of conformations, we evaluated the contacts in the entire receptor by means of Protein Structure Networks (PSNs)⁴¹. This method defines an interaction strength between the side chains of two residues, taking into account the number of atom pairs that are closer than a cutoff distance²⁸. Such interaction strengths were computed for all pairs of residues of the receptor, in the closed conformations of A2AR:Mini-Gs and in the open conformations of A2AR. This resulted in two matrices of average interaction strengths, one for each group. The difference between these two matrices highlighted the contribution of Mini-Gs to the structural changes between the two groups, illustrated in Fig. 4A. In this figure, pink cylinders represent interactions that prevail in the closed forms of A2AR:Mini-Gs and blue cylinders indicate predominant interactions in the open A2AR; the diameter of the cylinders is proportional to the difference of interaction strength between the two groups. We draw attention to the main features of this figure: (1) as expected, contacts located in the region of insertion of the α5 helix were strongly modified, but the effect of Mini-Gs extended beyond this region: it caused a global rearrangement of interactions, across the entire receptor; (2) some of the most prominent differences in interaction strength occurred around the agonist site and in the extracellular loop 2 (ECL 2) (as indicated by the thicker cylinders in these regions). We repeated the same analysis on closed conformations from the simulations of the A2AR:α5 complex with highly similar results (Fig. 4B), indicating that the interaction of the receptor with this helix is a key factor at the origin of the contact redistribution.
<figure><img src="image_7.png" /><figcaption>**Figure 3.** Conformations of the open and closed states of the agonist site. (A–D) 2D density plots showing the distribution of distances between the Co's of ALA 63 (2.61) and VAL 84 (3.32) relative to the perimeter of the triangle formed by the Co of THR 88 (3.36) and the side chains of TRP 246 (6.48) and HIS 278 (7.43). The free receptor explores closed, intermediate and open states. Mini-Gs or the α5 helix prevent the exploration of the open state and shift the equilibrium towards the agonist-bound conformation; (E) Superposition of representative conformations of the closed state of A2AR:Mini-Gs and the open state of A2AR. The side chains of TRP 246 and HIS 278 are shown in licorice representation and the Co's of ALA 63 (2.61), VAL 84 (3.32) and THR 88 (3.36) are shown as spheres.</figcaption></figure>

Notably, upon removal of Mini-Gs the side chain of ARG 102 (3.50) moved downward, loosening contacts with TYR 288 (7.53), which moved away from ILE 238 (6.40) (Fig. 4C). The latter came into contact with PHE 242 (6.44), which experienced a significant rotameric change and pointed down. This was accompanied by a relaxation of TRP 246 (6.46) (Fig. 4C). Rearrangements of these residues can also be observed in available X-Ray structures, when comparing A2AR coupled to Mini-Gs and uncoupled forms of the receptor. Though our sampling did not capture the complete conversion of the free receptor to its uncoupled conformation, even when the simulations were extended to 2 µs, we observed that ARG 102 (3.50) adopted conformations closer to those found in the inactive A2AR (Fig. S9).
<figure><img src="image_8.png" /><figcaption>**Figure 4.** Rearrangement of contacts in A2AR. (A) Differences between side chain interactions in A2AR and A2AR:Mini-Gs, revealed by analysis of Protein Structure Networks (PSNs). Pink cylinders: predominant contacts in A2AR:Mini-Gs; blue cylinders: predominant contacts in A2AR. The thickness of the cylinders is proportional to the difference in interaction strength. (B) Same analysis for the A2AR:α5 complex. (C) Conformational changes of side chains involved in the redistribution of contacts, that extended from the vicinity of the Mini-Gs binding site to the agonist site. Blue: A2AR; green: A2AR:Mini-Gs. The red arrows indicate the relaxation of W246, F242 and R102 upon removal of Mini-Gs.</figcaption></figure>

More surprisingly, Mini-Gs or the α5 helix also contributed to strengthen many interactions surrounding the agonist site, with reinforced connections between TM7 and TM2, TM4 and TM5 and TM6, 20 to 30 Å from the interaction site of the α5 helix. Residues TRP 246 (6.48) and HIS 278 (7.43), whose importance to the volume of the agonist site was shown (see Fig. 3), were also involved in some of these interactions: in the closed forms, TRP 246 (6.48) maintained contact with SER 277 (7.42) and HIS 278 (7.43) interacted more closely with ILE 60 (2.58), thus contributing to tighten the connections between TM6 and TM7 and between TM7 and TM2. We also observed that coupling to either partner was associated to a modification in the pattern of interactions of the ECL 2, with stronger contacts rigidifying part of this loop (Fig. S10).

As another evidence for the global and long-range effects of Mini-Gs or α5 helix binding, Fig. 5 shows 2D density plots that reveal the correlation between the conformation of residues next to the Mini-Gs binding site and the distance between TM2 and TM3, that in turn was determinant for the volume of the agonist site. The contacts between the side chains of TYR 288 (7.53) and ILE 238 (6.40) were enhanced in A2AR:Mini-Gs or A2AR:α5 relative to the free receptor and the closer interaction of these residues was associated to shorter distances between the Cα's of ALA 63 (2.61) and VAL 84 (3.32) in the closed forms of the complexes (Fig. 5A,B). In contrast, in the open state of A2AR a larger 238-288 distance was accompanied by a greater separation between TM2 and TM3 (that resulted in higher volumes of the agonist site; Fig. 5C). Interestingly, when the free receptor explored the closed state, the distance between TYR 288 (7.53) and ILE 238 (6.40) was also reduced (Fig. 5D). This result suggests that Mini-Gs or the α5 helix stabilize conformations of residues close to the intracellular domain that allosterically favor the closed conformation of the agonist site.
<figure><img src="image_9.png" /><figcaption>**Figure 5.** Coupling between the Mini-Gs and agonist binding sites. (A–D) 2D density plots of the distances between the Co's of residues A63 (2.61) and V84 (3.32), in the agonist site, and between the side chains of I238 (6.40) and Y288 (7.53), close to the region of insertion of the α5 helix. The closed state of the site is associated to a shorter I238-Y288 distance, even in the free receptor.</figcaption></figure>

To further explore this hypothesis we simulated a double mutant of A2AR, I238A/Y288A, coupled to the α5 helix and measured the volume distribution of the agonist site. In this case, α5 lost the ability to stabilize the closed conformation of the site. The mutated receptor visited open and intermediate states, with large values of the T88-W246-H278 perimeter, while the closed conformation was suppressed. In comparison to A2AR:α5, a pronounced outward motion of TM2 also contributed to open the receptor. These results are shown in Figs S11 and S12 and reinforced the suggestion that the local rearrangement of contacts in the region of interaction with α5 triggers global changes that reach the orthosteric site.

Finally, we also built other PSNs of the receptor, taking into account the whole trajectories and not only con-
formations from specific states. PSNs were then computed for each of the four different studied systems. After
being subjected to PCA, the matrices of average interaction strengths were clustered, confirming that the α5 helix
was able to mimic the global effect of the Mini-Gs with respect to the pattern of contacts (Fig. 6).

**Discussion**

In this work we employed molecular dynamics simulations to study the allosteric effects of the Mini-Gs protein or its C-terminal α5 helix upon the A2 adenosine receptor. Our simulations confirmed the hypothesis that the G-protein stabilizes a closed form of the receptor even in the absence of agonist and showed that Mini-Gs or α5 were able to shift the dynamical equilibrium towards a closed conformation of the orthosteric site, as suggested by previous experiments²⁰.

Considering that all our simulations started from the closed state of the receptor, Mini-Gs or α5 prevented or slowed down a complete relaxation of the agonist site towards an open conformation, even in the absence of
Clusterization Based on PCA of Interaction Matrices

<figure><img src="image_10.png" /><figcaption>**Figure 6.** Principal Component Analysis of matrices of average interaction strength of A2AR. Trajectories of A2AR coupled to Mini-Gs or to the α5 helix in the absence of NECA were grouped together, both being closer to the NECA bound form than to the free receptor. This indicates that α5 could maintain the pattern of contacts of the Mini-Gs coupled receptor.</figcaption></figure>

ligand. The fact that α5 was able to reproduce the global effect of Mini-Gs suggests that the reorganization of contacts in the G protein binding site is essential to trigger the allosteric signal. It is also consistent with the experimental observation that an antibody that mimics the G protein also leads to a closed conformation of the agonist site²⁰. The G protein, the α5 helix and the antibody share the same binding site of the receptor, possibly leading to similar local rearrangements of contacts that eventually affect the entire receptor.

In our simulations, this local rearrangement involved stronger interactions between I238 (6.40) and Y288 (7.53) in the presence of Mini-Gs or α5. Simulations with the double mutant I238A/Y288A confirmed the importance of this contact to the allosteric effect, since its disruption in A2AR:α5 abolished the closed state of the receptor. Y288 is part of the conserved NpxxY motif and, together with Y197, is possibly involved in an activation switch in class A GPCRs⁴². In the A3R, the Y288F mutation led to a 30-fold reduction in potency of a highly selective agonist, suggesting a role of this residue in the dynamical equilibrium of the receptor⁴³. The mutation I238A in A2AR led to an increased thermostability of agonist- and inverse-agonist-bound conformations⁴⁴, indicating that this residue can also be important for conformational equilibrium. Experiments with the double mutant I238A/Y288A would be valuable to confirm the pertinence of this theoretical model.

Our analyses showed that an important effect of Mini-Gs was the strengthening of contacts in ECL2, reducing its flexibility. The roles of ECL2 in ligand selectivity and stabilization of agonist-bound conformations in A2AR have been previously discussed⁴⁵. Our findings suggest that the allosteric effect of the G-protein on ECL2 could contribute to such stabilization. In particular, E169 is known to be involved in ligand binding⁴² and its interaction with the agonist could be favored by a reduced flexibility of ECL2.

Although it is tempting to discuss our observations in the framework of allosteric communication pathways, we avoided this description in view of results that point to the technical and conceptual difficulties in defining these paths. Hilser and co-workers have argued about the limitations of the concept of pathways in the context of an ensemble view of allostery⁴⁶,⁴⁷. More recently, simulations permitted a time-resolved study of communication in an allosteric protein, revealing that different trajectories had significant heterogeneity in the propagation of the allosteric signal and that nonlocal and simultaneous structural changes could occur in a single trajectory⁴⁸. Our sampling was not extensive enough for such a detailed study of the precise temporal order of conformational changes. This limited sampling was also reflected in the incomplete relaxation of the free receptor to its inactive conformation, even when the simulations were extended to 2 µs. Nevertheless, we captured rearrangements of residues such as W246, F242 and R102, this last one adopting conformations closer to those found in the inactive A2AR.

Moreover, our results suggest that the observed allosteric effect could result from a complex interplay between the entire network of contacts, that is globally rewired by Mini-Gs or the α5 helix. In spite of this complexity, it was possible to identify interactions that were stabilized in the different tested conditions, even if distinct allosteric pathways or mechanisms could contribute to such stabilization.

Our observations suggest that the closed state would be further stabilized by the binding of the agonist. Therefore, similarly to an agonist that shifts the equilibrium towards the active conformation, that is stabilized by the G protein, the G protein could select conformations more prone to interact favorably with the agonist, that would stabilize them when binding to the orthosteric site. This perspective could be particularly relevant in the study of receptor-ligand interactions in cases where the GPCR:G protein complex is pre-assembled, as recently shown for the ghrelin receptor⁴⁹. In these situations, docking experiments, for example, could benefit from an appropriate selection of more closed conformations of the receptor.
It has been suggested that a common allosteric mechanism could operate on class A GPCRs in order to close the agonist site upon coupling to the G protein<sup>20</sup>. This hypothesis cannot be excluded by our results, but they also point to the importance of considering the specificities of the distinct receptors. For example, in recent simulations, Miao and McCammon verified that the reduction in the volume of the orthosteric site of the M2 muscarinic receptor was associated to backbone motions of TM3, TM6 and TM7<sup>50</sup>. In our simulations of the A2AR, the volume was mainly affected by backbone motions of TM2 and, to a lesser extent, TM3 and TM7, with an important contribution of side chains, too. Even if the underlying mechanisms may differ, the theoretical and experimental results available so far point to the role of G protein in restricting the size of the orthosteric site of the receptor.

While this paper was under review, Lee *et al*.⁵¹ published an article about the stabilisation of the active state of A2AR by the G protein. They performed molecular dynamics simulations on four agonist-bound conformations of the receptor, including the complex with Mini-Gs. In contrast to our own study, they have not studied ligand-free forms of the receptor and carefully evaluated the interaction between the protein and the agonist in each case. They observed that coupling of Mini-Gs was associated to tighter contacts with the agonist and more favorable binding free energies⁵¹. In agreement with our findings, they point to the stabilisation of a closed conformation of orthosteric site by Mini-Gs⁵¹. They have also reported a stronger allosteric coupling between extra- and intracellular regions in the presence of Mini-Gs, mediated through a high level of correlated motions mainly in TM5⁵¹. In our investigation, we pointed to important side chain rearrangements in TM6 (Fig. 4C); however, we note that we analysed a distinct situation: instead of correlated motions in a system, our residue interaction analysis looked for differences in the conformations of side chains between A2AR and A2AR:Mini-Gs (or A2AR:α5). Despite the different approaches, both of our studies demonstrated an allosteric effect exerted by the G protein and the importance of dynamics in its evaluation.

## Methods

**Model preparation.** Models were based on the crystal structure of A2AR bound to Mini-Gs, PDB id: 5G53, pre-aligned in the Orientations of Proteins in Membranes database⁵². The chain B of this file was the basis for the A2AR construction and chain D was chosen to build the Mini-Gs model. Addition of missing residues and loop refinement were performed with Modeller⁵³. Missing residues in the receptor were located in the intracellular loop 3 (ICL3) region (residues 208–233); in Mini-Gs, residues 224–239 and a stretch of 15 aminoacids between residues 61 and 208 were also missing (note that the numbers refer to the numbering scheme in the original pdb file, 5G53).

Models for A2AR and Mini-Gs were built separately and reassembled through structural alignment to the original pdb file (in Pymol). The double mutant of the receptor (I238A/Y288A) was prepared using the Pymol Mutagenesis Wizard tool. GDP, initially bound to Mini-Gs, and NECA were included in the final model of the A2AR:Mini-Gs:NECA complex. This model was the starting point of the other simulated systems, constructed by removal of the corresponding components.

Models were inserted in 100% POPC membranes, solvated and ionized up to a concentration of 0.15 M NaCl using the Charmm-gui web server⁵⁴. Systems containing Mini-Gs were composed of ~105000 atoms (~70000 water atoms and ~27000 lipid atoms in ~200 POPC molecules); systems without Mini-Gs or with the receptor bound to the α5 helix were composed of ~79000 atoms (~50000 water atoms and ~23500 lipid atoms in ~175 POPC molecules).

**Molecular Dynamics Simulations.** Two completely independent replicas were launched for each system. Systems were equilibrated following the previously described CHARMM-GUI protocol⁵⁵,⁵⁶. Briefly, after energy minimization, successive steps of NVT and NPT (300 K, 1 bar) MD were performed, with progressive removal of position restraints applied to the membrane and protein atoms. Then, simulations were run for 600 ns (NPT ensemble, 300 K, 1 bar) and the first 100 ns were considered as equilibration and discarded; the last 500 ns of each replica were retained for analysis. Furthermore, for all systems except A2AR:α5 and the double mutant of A2AR coupled to α5, production simulations were extended until 2.0 µs.

Proteins, lipids and ions were described by the CHARMM36 force field⁵⁷⁻⁵⁹ and the parameters for GDP and NECA were obtained with the CGenFF force field⁶⁰, through the CHARMM-GUI interface. A time step of 2 fs was used in the production phase and PME (Particle Mesh Ewald)⁶¹ was employed for the treatment of long-range electrostatic interactions, with the application of a switch function between 1.0 and 1.2 nm. Temperature was kept at 300 K by the Nose-Hoover scheme⁶²,⁶³, using a time constant for coupling of 1 ps. Pressure was maintained at 1 bar by a semi-isotropic Parrinello-Rahman barostat⁶⁴, with coupling time constant of 5 ps and compressibility of 4.5 × 10⁻⁵ bar⁻¹. All bonds to hydrogen atoms were constrained by the LINCS algorithm⁶⁵. All simulations were performed with the Gromacs 5.1.4 software⁶⁶.

**Preparation for Analysis.** The coordinates of the receptor were extracted from all trajectories and concatenated. All conformations of the receptor in this “super trajectory”, containing snapshots simulated in all conditions, were superimposed according to a procedure that identifies a rigid core formed by the most invariant regions in terms of dynamics⁶⁷. Only the backbone of the residues was considered in this alignment procedure. The residues of the receptor identified as belonging to this most invariant core are listed in the Text S2. This core was used as the reference for the alignment of all snapshots and the initial model, based on the crystal structure (PDB id: 5G53).

**Principal Component Analysis of the Residues of the Agonist Binding Site.** The concatenated and aligned trajectory of the receptor was submitted to a Principal Component Analysis, consisting in the diagonalization of the covariance matrix, whose elements are:
$$C_{ij} = \langle (\mathbf{r}_i - \langle \mathbf{r}_i \rangle) \cdot (\mathbf{r}_j - \langle \mathbf{r}_j \rangle) \rangle$$

where i and j denote all pairs of the 3N Cartesian coordinates of the N atoms of the residues of the agonist site, including their side chains. The vector $\mathbf{r}_i$ indicates the instantaneous value of coordinate i and $\langle \mathbf{r}_i \rangle$ is the average value of this coordinate in the ensemble of conformations.

The coordinates of the selected residues in the snapshots corresponding to each individual trajectory were then projected on the space of the respective first two principal components. These analyses were performed with the Gromacs software<sup>66</sup>.

**Distances between residues.** The Mdtraj software<sup>68</sup> was used to measure the Euclidean distance between atoms. For a given pair of residues, two distances were measured in every trajectory frame: (i) between the Co's and (ii) between the side chains. In the latter case, it consisted in the distance between the closest heavy atoms in the side chains. To establish the perimeter of the triangle formed by residues 88, 246 and 278, we have additionally determined the distances between the CH2 atom in W246 and the CD2 atom in H278 (the closest side chain atoms in the initial conformation, extracted from pdb file 5G53) and also between the Co of T88 and each one of these atoms.

**Definition of open and closed states of the receptor.** The definition of states was based on the distance between the Co's of residues ALA 63 (2.61) and VAL 84 (3.32) and the perimeter of the triangle formed by the Co of THR 88 (3.36), the CH2 atom of TRP 246 (6.48) and the CD2 atom in HIS 278 (7.43). The following criteria were used to define the states:

Closed: perimeter T88-W246-H278 < 16 Å and distance A63-V84 < 8 Å

Open: perimeter T88-W246-H278 > 21 Å and distance A63-V84 > 13 Å

Otherwise, the states were considered as intermediate.

**Volume of the Agonist Site.** Direct estimations of the volume of the agonist site were obtained with the Epock software<sup>69</sup>. Before volume calculation, this program requires the definition of a region of interest (the "maximum encompassing region"). This region was defined as a sphere of 10 Å radius, centered in the geometric center of the residues of the agonist site, listed in the Supporting Information (Text S1). Please refer to the original publication for further details about the method<sup>69</sup>.

**Protein Structure Networks.** PSNs were obtained by the method introduced by Brinda and Vishveshwara<sup>41</sup>, using the Wordom software<sup>70</sup>. In this approach, an interaction strength between two amino acid side chains i and j is given by:

$$I_{ij} = (n_{ij}/\sqrt{N_i \times N_j}) \times 100,$$

where the result is expressed as a percentage. $n_{ij}$ is the number of atom pairs between the side chains that are closer than a cutoff, that was chosen as 4.5 Å; $N_i$ and $N_j$ are normalization factors that depend on the types of residues i and j (details in the original publication). Residues separated by 3 or less residues in the sequence were considered adjacent and skipped in the contact calculation. PSNs were calculated for all conformations in a given condition, resulting in average interaction strengths over the trajectories. Then, for each condition, one matrix with the interaction strengths between all residue pairs was obtained. These were the matrices used in subsequent analyses.

**PCA and clusterization of interaction strenght matrices.** Interaction strenght matrices calculated for the four simulated conditions (A2AR, A2AR:NECA, A2AR:Mini-Gs, A2AR:α5) were subjected to Principal Component Analysis with the Bio3d software<sup>67</sup>. The matrices were then clustered in PC space with the single linkage method.

## References

1. Katritch, V., Cherezov, V. & Stevens, R. C. Structure-function of the G protein-coupled receptor superfamily. *Annu. Rev. Pharmacol. Toxicol.* **53**, 531–556 (2013).

2. Marie, S. et al. Ligands and signaling proteins govern the conformational landscape explored by a G protein-coupled receptor. *Proc. Natl. Acad. Sci.* **109**, 8304–8309 (2012).

3. Manglik, A. & Kobilka, B. The role of protein dynamics in GPCR function: insights from the β2AR and rhodopsin. *Curr. Opin. Cell Biol.* **27**, 136–143 (2014).

4. De Lean, A., Stadel, J. M. & Lefkowitz, R. J. A ternary complex model explains the agonist-specific binding properties of the adenylate cyclase-coupled β-adrenergic receptor. *J. Biol. Chem.* **255**, 7108–7117 (1980).

5. West, G. M. et al. Ligand-dependent perturbation of the conformational ensemble for the GPCR β2 adrenergic receptor revealed by HDX. *Structure* **19**, 1424–1432 (2011).

6. Nygaard, R. et al. The dynamic process of β2-adrenergic receptor activation. *Cell* **152**, 532–542 (2013).

7. Kruse, A. C. et al. Activation and allosteric modulation of a muscarinic acetylcholine receptor. *Nature* **504**, 101–106 (2013).

8. Weis, W. I. & Kobilka, B. K. The Molecular basis of G protein-coupled receptor activation. *Annu. Rev. Biochem.* **87**, 897–919 (2018).

9. Louet, M. et al. Conformational restriction of G proteins coupled receptors (GPCRs) upon complexation to G proteins: a putative activation mode of GPCRs? *FEBS Lett.* **587**, 2656–2661 (2013).

10. Miao, Y. & McCammon, J. A. G protein coupled receptors: advances in simulation and drug discovery. *Curr. Opin. Struct. Biol.* **41**, 83–89 (2016).

11. Provasi, D., Artacho, M. C., Negri, A., Mobarec, J. C. & Filizola, M. Ligand-Induced modulation of the free-energy landscape of G protein-coupled receptors explored by adaptive biasing techniques. *PLoS Comput. Biol.* **7**, e1002193 (2011).

12. Ferruz, N. et al. Dopamine D3 receptor antagonist reveals a cryptic pocket in aminergic GPCRs. *Sci. Rep.* **8**, 1–10 (2018).

13. Li, J., Jonsson, A. L., Beuming, T., Shelley, J. C. & Voth, G. A. Ligand-dependent activation and deactivation of the human adenosine A2A receptor. *J. Am. Chem. Soc.* **135**, 8749–8759 (2013).
14. Rasmussen, S. G. F. et al. Crystal structure of the β2 adrenergic receptor-Gs protein complex. *Nature* **477**, 549–557 (2011).

15. Koehl, A. et al. Structure of the μ-opioid receptor-Gi protein complex. *Nature* **558**, 547–552 (2018).

16. Kang, Y. et al. Cryo-EM structure of human rhodopsin bound to an inhibitory G protein. *Nature* **558**, 553–558 (2018).

17. García-Nafria, J., Nehmé, R., Edwards, P. C. & Tate, C. G. Cryo-EM structure of the serotonin 5-HT1B receptor coupled to heterotrimeric Go. *Nature* **558**, 620–623 (2018).

18. Flock, T. et al. Universal allosteric mechanism for Gα activation by GPCRs. *Nature* **524**, 173–179 (2015).

19. Saleh, N., Saladino, G., Gervasio, F. L. & Clark, T. Investigating allosteric effects on the functional dynamics of β2-adrenergic ternary complexes with enhanced-sampling simulations. *Chem. Sci.* **8**, 4019–4026 (2017).

20. Devree, B. T. et al. Allosteric coupling from G protein to the agonist-binding pocket in GPCRs. *Nature* **535**, 182–186 (2016).

21. Carpenter, B. & Tate, C. G. Engineering a minimal G protein to facilitate crystallisation of G protein-coupled receptors in their active conformation. *Protein Eng. Des. Sel.* **29**, 583–593 (2016).

22. Carpenter, B., Nehmé, R., Warne, T., Leslie, A. G. W. & Tate, C. G. Structure of the adenosine A2A receptor bound to an engineered G protein. *Nature* **536**, 104–107 (2016).

23. de Lera Ruiz, M., Lim, Y. H. & Zheng, J. Adenosine A2A receptor as a drug discovery target. *J. Med. Chem.* **57**, 3623–3650 (2014).

24. Jacobson, K. A. & Gao, Z. G. Adenosine receptors as therapeutic targets. *Nature Rev. Drug Discov.* **5**, 247–264 (2006).

25. Jaakola, V. P. et al. The 2.6 angstrom crystal structure of a human A2A adenosine receptor bound to an antagonist. *Science* **322**, 1211–1217 (2008).

26. Dore, A. S. et al. Structure of the adenosine A(2A) receptor in complex with ZM241385 and the xanthines XAC and caffeine. *Structure* **19**, 1283–1293 (2011).

27. Liu, W. et al. Structural basis for allosteric regulation of GPCRs by sodium ions. *Science* **337**, 232–236 (2012).

28. Segala, E. et al. Controlling the dissociation of ligands from the adenosine A2A receptor through modulation of salt bridge strength. *J. Med. Chem.* **59**, 6470–6479 (2016).

29. Cheng, R. K. Y. et al. Structures of human A1 and A2A adenosine receptors with xanthines reveal determinants of selectivity. *Structure* **25**, 1275–1285.e4 (2017).

30. Sun, B. et al. Crystal structure of the adenosine A2A receptor bound to an antagonist reveals a potential allosteric pocket. *Proc. Natl. Acad. Sci.* **114**, 2066–2071 (2017).

31. Lebon, G. et al. Agonist-bound adenosine A2A receptor structures reveal common features of GPCR activation. *Nature* **474**, 521–525 (2011).

32. Xu, F. et al. Structure of an agonist-bound human A2A adenosine receptor. *Science* **332**, 322–327 (2011).

33. Lebon, G., Edwards, P. C., Leslie, A. G. & Tate, C. G. Molecular determinants of CGS21680 binding to the human adenosine A2A receptor. *Mol. Pharmacol.* **87**, 907–915 (2015).

34. Tate, C. G et al. Cryo-EM structure of the adenosine A2A receptor coupled to an engineered heterotrimeric G protein. *bioRxiv*, 267674, https://doi.org/10.1101/267674 (2018).

35. Rasenick, M. M., Watanabeg, M., Lazarevic, M. B., Hatta, S. & Hamm, H. E. Synthetic peptides as probes for G protein function. *J. Biol. Chem.* **269**, 21519–21525 (1994).

36. Semack, A., Sandhu, M., Malik, R. U., Vaidehi, N. & Sivaramakrishnan, S. Structural elements in the Gαs and Gαq C-termini that mediate selective G protein-coupled receptor (GPCR) signaling. *J. Biol. Chem.* **291**, 17929–17940 (2016).

37. Ballesteros, J. A. & Weinstein, H. [19] Integrated methods for the construction of three-dimensional models and computational probing of structure-function relations in G protein-coupled receptors. In *Methods in Neurosciences*, Vol. 25 (ed. Stuart, C. S.) 366–428 (Academic Press, 1995).

38. Preininger, A. M., Meiler, J. & Hamm, H. E. Conformational flexibility and structural dynamics in GPCR-mediated G protein activation: a perspective. *J. Mol. Biol.* **425**, 2288–2298 (2013).

39. Louet, M., Perahia, D., Martinez, J. & Floquet, N. A concerted mechanism for opening the GDP binding pocket and release of the nucleotide in hetero-trimeric G proteins. *J. Mol. Biol.* **411**, 298–312 (2011).

40. Venkatakrishnan, A. J. et al. Molecular signatures of G-protein-coupled receptors. *Nature* **494**, 185–194 (2013).

41. Brinda, K. V. & Vishveshwara, S. A network representation of protein structures: implications for protein stability. *Biophys. J.* **89**, 4159–4170 (2005).

42. Jespers, W. et al. Structural mapping of adenosine receptor mutations: ligand binding and signaling mechanisms. *Trends. Pharmacol. Sci.* **39**, 75–89 (2018).

43. Chen, A. et al. Constitutive activation of A3 adenosine receptors by site-directed mutagenesis. *Biochem. Biophys. Res. Commun.* **284**, 596–601 (2001).

44. Lebon, G., Bennett, K., Jazayeri, A. & Tate, C. G. Thermostabilisation of an agonist-bound conformation of the human adenosine A2A receptor. *J. Mol. Biol.* **409**, 298–310 (2011).

45. Seibt, B. F. et al. The second extracellular loop of GPCRs determines subtype-selectivity and controls efficacy as evidenced by loop exchange study at A2 adenosine receptors. *Biochem. Pharmacol.* **85**, 1317–1329 (2013).

46. Pan, H., Lee, J. C. & Hilser, V. J. Binding sites in Escherichia coli dihydrofolate reductase communicate by modulating the conformational ensemble. *Proc. Natl. Acad. Sci.* **97**, 12020–12025 (2000).

47. Motlagh, H. N., Wrabl, J. O., Li, J. & Hilser, V. J. The ensemble nature of allostery. *Nature* **508**, 331–339 (2014).

48. Buchenberg, S., Sittel, F. & Stock, G. Time-resolved observation of protein allosteric communication. *Proc. Natl. Acad. Sci.* **114**, E6804–E6811 (2017).

49. Damian, M. et al. Ghrelin receptor conformational dynamics regulate the transition from a preassembled to an active receptor: Gq complex. *Proc. Natl. Acad. Sci.* **112**, 1601–1606 (2015).

50. Miao, Y. & McCammon, J. A. Mechanism of the G protein mimetic nanobody binding to a muscarinic G protein-coupled receptor. *Proc. Natl. Acad. Sci.* **115**, 3036–3041.

51. Lee, S., Nivedha, A. K., Tate, C. G. & Vaidehi, N. Dynamic role of the G Protein in stabilizing the active state of the adenosine A2A receptor. *Structure* **27**, 1–10 (2019).

52. Lomize, M. A., Lomize, A. L., Pogozheva, I. D. & Mosberg, H. I. OPM: orientations of proteins in membranes database. *Bioinformatics* **22**, 623–625 (2006).

53. Sali, A. & Blundell, T. L. Comparative protein modelling by satisfaction of spatial restraints. *J. Mol. Biol.* **234**, 779–815 (1993).

54. Wu, E. L. et al. CHARMM-GUI membrane builder toward realistic biological membrane simulations. *J. Comput. Chem.* **35**, 1997–2004 (2014).

55. Jo, S., Kim, T. & Im, W. Automated builder and database of protein/membrane complexes for molecular dynamics simulations. *PLoS One* **2**, e880 (2007).

56. Jo, S., Lim, J. B., Klauda, J. B. & Im, W. CHARMM-GUI membrane builder for mixed bilayers and its application to yeast membranes. *Biophys. J.* **97**, 50–58 (2009).

57. Best, R. B. et al. 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. Theory Comput.* **8**, 3257–3273 (2012).

58. Klauda, J. et al. Update of the CHARMM all-atom additive force field for lipids: validation on six lipid types. *J. Phys. Chem. B* **114**, 7830–7843 (2010).

59. Lee, J. et al. CHARMM-GUI input generator for NAMD, GROMACS, AMBER, OpenMM, and CHARMM/OpenMM simulations using the CHARMM36 additive force field. *J. Chem. Theory Comput.* **12**, 405–413 (2016).
60. Vanommeslaeghe, K. et al. CHARMM general force field: a force field for drug-like molecules compatible with the CHARMM all-atom additive biological force fields. *J. Comput. Chem.* **31**, 671–690 (2010).

61. Darden, T., York, D. & Pedersen, L. Particle mesh Ewald: An N-log(N) method for Ewald sums in large systems. *J. Chem. Phys.* **98**, 10089–10092 (1993).

62. Nosé, S. A unified formulation of the constant temperature molecular dynamics methods. *J. Chem. Phys.* **81**, 511–519 (1984).

63. Hoover, W. G. Canonical dynamics: equilibrium phase-space distributions. *Phys. Rev. A* **31**, 1695–1697 (1985).

64. Parrinello, M. & Rahman, A. Polymorphic transitions in single crystals: a new molecular dynamics method. *J. Appl. Phys.* **52**, 7182–7190 (1981).

65. Hess, B., Bekker, H., Berendsen, H. J. C. & Fraaije, J. G. E. M. LINCS: a linear constraint solver for molecular simulations. *J. Comput. Chem.* **18**, 1463–1472 (1997).

66. Abraham, M. J. et al. GROMACS: High performance molecular simulations through multi-level parallelism from laptops to supercomputers. *SoftwareX* **1–2**, 19–25 (2015).

67. Grant, B. J., Rodrigues, A. P. C., ElSawy, K. M., McCammon, J. A. & Caves, L. S. D. Bio3d: an R package for the comparative analysis of protein structures. *Bioinformatics* **22**, 2695–2696 (2006).

68. McGibbon, R. T. et al. MD Traj: a modern open library for the analysis of molecular dynamics trajectories. *Biophys. J.* **109**, 1528–1532 (2015).

69. Laurent, B. et al. Epock: rapid analysis of protein pocket dynamics. *Bioinformatics* **31**, 1478–1480 (2015).

70. Seeber, M. et al. Wordom: a user-friendly program for the analysis of molecular structures, trajectories, and free energy surfaces. *J. Comput. Chem.* **32**, 1183–1194 (2011).

## Acknowledgements

We thank MESO@LR, CINES (DARI project n° A0020707530) and PRACE (project n° 2017174234) for calculation resources. This work has been funded by ANR (Agence nationale de la recherche) via the projects BIP:BIP (ANR 10-BINF-0003) and ALLOSIG (ANR-17-CE11-0011); we thank the CHEMISYST labex for financial support for PR and the University of Montpellier for PR.

## Author Contributions

P.R., M.L. and N.F. carried out the simulations and analyzed the results. P.R., M.L., J.M. and N.F. designed the experiments. G.L. and N.F. obtained funding support and supervised the overall project. All authors contributed to manuscript preparation.

## Additional Information

Supplementary information accompanies this paper at https://doi.org/10.1038/s41598-019-41980-x.

**Competing Interests:** The authors declare no competing interests.

**Publisher's note:** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

<figure><img src="image_11.png" /></figure>

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.

© The Author(s) 2019