Dynamic behavior of the active and inactive states of the adenosine A(2A) receptor.

The adenosine A2A receptor (A2AR) belongs to the superfamily of membrane proteins called the G-protein-coupled receptors (GPCRs) that form one of the largest superfamilies of drug targets. Deriving thermostable mutants has been one of the strategies used for crystallization of A2AR in both the agonist and antagonist bound conformational states. The crystal structures do not reveal differences in the activation mechanism of the mutant receptors compared to the wild type receptor, that have been observed experimentally. These differences stem from the dynamic behavior of the mutant receptors. Furthermore, it is not understood how the mutations confer thermostability. Since these details are difficult to obtain from experiments, we have used atomic level simulations to elucidate the dynamic behavior of the agonist and antagonist bound mutants as well the wild type A2AR. We found that significant enthalpic contribution leads to stabilization of both the inactive state (StaR2) and active-like state (GL31) thermostable mutants of A2AR. Stabilization resulting from mutations of bulky residues to alanine is due to the formation of interhelical hydrogen bonds and van der Waals packing that improves the transmembrane domain packing. The thermostable mutant GL31 shows less movement of the transmembrane helix TM6 with respect to TM3 than the wild type receptor. While restricted dynamics of GL31 is advantageous in its purification and crystallization, it could also be the reason why these mutants are not efficient in activating the G proteins. We observed that the calculated stress on each residue is higher in the wild type receptor compared to the thermostable mutants, and this stress is required for activation to occur. Thus, reduced dynamic behavior of the thermostable mutants leading to lowered activation of these receptors originates from reduced stress on each residue. Finally, accurate calculation of the change in free energy for single mutations shows good correlation with the change in the measured thermostability. These results provide insights into the effect of mutations that can be incorporated in deriving thermostable mutants for other GPCRs.

Introduction

G-protein-coupled receptors (GPCRs) belong to the superfamily of membrane proteins and play a crucial role in signal transduction. They also form the largest class of drug targets for many diseases.[1] GPCRs share a similar structural motif consisting of the seven helical transmembrane (TM) core connected by three extracellular and three intracellular loops.[2] GPCRs are dynamic proteins with several functionally important conformations ranging from the “inactive state” to the “fully active” state. The dynamics of GPCR conformations[3−5] pose a challenge for their purification and crystallization in various functional states, and site directed mutagenesis that confers thermostability to the receptor has been a successful strategy to stabilize these proteins in various conformations in detergents.[6−8] The adenosine A2A receptor (A2AR) is a class A GPCR and an emerging drug target for the treatment of Parkinson’s disease, inflammation, and cardiac ischemia.[9−12] A2AR has been thermostabilized and crystallized in two functional conformations.[13,14] A2AR-StaR2 is the inactive state mutant with eight point mutations,[13,15] while A2AR-GL31 with four point mutations is the thermostabilized mutant stabilized in an “active-like” state.[16] Hereafter, we refer to the inactive state mutant as just StaR2 and the mutant in the active-like state as GL31. The eight point mutations in StaR2 are A54L2.52, T88A3.36, R107A3.55, K122A4.43, L202A5.63, L235A6.37, V239A6.41, and S277A7.42. The four point mutations in GL31 are L48A2.46, A54L2.52, T65A2.63, and Q89A3.37. Here we have used the Ballesteros–Weinstein amino acid numbering system used for class A GPCRs.[17] The first number in the superscript is the TM helix in which the amino acid is present, and the second number is the position of this residue with respect to the most conserved residue in that helix which is numbered 50.

Analysis of the crystal structures of the inactive state[18] and the active state of the β2-adrenergic receptor with the G protein bound,[19] shows that the transmembrane helices TM5 and TM6 move significantly with respect to TM3 upon activation (when bound to both agonist and the G protein). The crystal structure of the agonist bound GL31 shows a similar type of movement of TM6 with respect to TM3 but not as far as observed in the active state of the β2-adrenergic receptor. Similar limited movement of TM6 has also been observed in the crystal structure of the agonist bound wild type A2AR.[20] Hence, we call the conformation of GL31 the “active-like” state henceforth in the paper.

While the crystal structures of the inactive and the active-like states of A2AR show conformational differences, they do not provide the answers for two important questions: (1) how do the mutations stabilize the receptor, and (2) what are the differences in the dynamics of the activation mechanism of the thermostable mutant compared to the wild type? We have used atomic level molecular dynamics (MD) simulations to answer these questions, since it is not straightforward to answer using experimental techniques. While some of the mutations in the thermostable A2AR confer stability, other mutations specifically stabilize the agonist bound state compared to the inactive state.[21] The role of the mutations and how certain mutations specifically stabilize the active-like state compared to the inactive state is not well understood. The insight provided by how certain mutations stabilize the active-like state compared to the inactive state would enormously benefit future design of thermostable mutants for other GPCRs. Additionally, the difference between the dynamics of the active-like state mutant compared to the wild type would provide insight into how the mutations in GL31 limit the receptor ability to activate G proteins although it is in the active-like conformation.

Methods

System Setup for MD Simulations

The starting conformations of A2AR for the MD simulations were taken from the crystal structures of GL31 and StaR2 (pdb ID: 2YDO for GL31 and 3PWH for StaR2). Hydrogens were added, the structures were solvated in the explicit palmitoyloleoylphosphatidylcholine (POPC) lipid and water, and the solvent was packed using the inflategro package in GROMACS.[22] Two initial conformations of the wild type were generated by mutating the residues in the crystal structure of the mutants back to the wild type using Maestro9.2.[23] We did not use the crystal structures of the wild type A2AR, since those structures were crystallized with T4L. However, we also performed simulations starting from the wild type A2AR crystallized with T4L.[24] We found no significant difference in the dynamics (Figure S1, Supporting Information) and hence pursued with the wild type receptors derived from the thermostable mutant structures. Residues within 5 Å of the sites of mutation were minimized using MacroModel with position restraints on all backbone atoms and all residues further than 5 Å from the site of mutation. The wild type generated from GL31 mutant (GL31struc-WTseq) is the “active-like” wild type conformation, and that generated from StaR2 (StaR2struc-WTseq) is the “inactive state” conformation of the wild type. MD simulations on A2AR in a POPC lipid bilayer with periodic boundary conditions were performed using the GROMACS package with the GROMOS96 force field[25] with SPC water molecules.[26] The SETTLE[27] and LINCS[28] algorithms were used for the bond and angle for water and all other bonds, allowing 2 fs of time step. For the analysis, the coordinates were saved every 2 ps. A cutoff distance of 12 Å for nonbond interactions was introduced, and the PME (particle mesh Ewald) method[29,30] was used for long-range vdW interactions. We performed MD simulations on eight systems, each 1 μs long. The eight systems are shown in Table S1 of the Supporting Information.

Each of the eight systems were equilibrated by performing 200 ps of MD at 310 K using a NVT ensemble followed by 5 ns of MD under NPT conditions with a pressure of 1 bar. The protein and ligand were kept in place during these equilibration steps using position restraints. After equilibration to the expected temperature and pressure, a total of 10 production simulations of up to 100 ns were performed for each initial conformation with different initial velocities using the NVT ensemble.

Trajectory Analysis

All trajectories obtained from the molecular dynamics simulations were analyzed using tools provided by GROMACS and Python script. PyMOL[31] and VMD[32] were used for the structural conformation analysis of the trajectories. For hydrogen bond analysis using 3.9 Å and 30° for the cutoff distance and angle, respectively, the g_hbond utility of GROMACS was used. Interhelical hydrogen bond interactions were derived from the stable hydrogen bond criteria having more than 50% occupancy (population) through 1 μs trajectories, which is normalized by setting the most densely populated point to 1.

Free Energy Simulations

The change in the free energy for single point mutations was calculated using the thermodynamic integration (TI) technique.[33,34] We used the Hamiltonian shown below which is a function of a coupling parameter λ that varies progressively from 0 to 1 to change the system interactions from the initial state (A) to the final state which is the single point mutation (B). The difference in free energy between the molecular systems A and B is then calculated usingThe simulations were done at a few discrete points λ along the pathway, and the integral was calculated numerically. Using the final snapshot of wild type receptor after 1 μs of MD trajectory as initial conformation, the free energy changes due to the mutations were calculated for both the wild type state (λ = 0) and the single mutant state (λ = 1) by integrating the average enthalpic contribution from each window. For the better convergence of the change in free energy, we used unequal window intervals: Δλ = 0.01 (λ = 0.00–0.10, λ = 0.90–1.00) and Δλ = 0.02 (λ = 0.10–0.90). Since the error estimation of the free energy generally is crucial, we monitored the standard deviation of the total calculations, expressing low fluctuation in the free energy change (±1.5 kcal/mol). To obtain the best estimates of the free energies, ΔG (WT → Mut), all the enthalpic contributions from each window from the entire 60 ns (12 ns of equilibration + 48 ns of production) were used. The free energy differences have been calculated using a united atom forcefield GROMOS, and this might affect the accuracy of the calculated free energies.

Stress Calculation

The residue-based stress (forces) was calculated using both bond and nonbonded force on each residue coming from all residues within 3 Å of this residue except the ones that are directly bonded to the target residue. The force computation was performed using the GROMOS96 53a6 force field following the procedure described by Stacklies et al.[35] The average stress is the average residue-based stress over the entire MD trajectory for each system. The procedure is discussed in detail in Niesen et al.[36]

PCA of MD Trajectories

To understand the most important collective motion in the receptor within a few dominant modes, we performed the principal component analysis (PCA) over the entire 1 μs MD simulation trajectory. Only Cα atoms in transmembrane helices were included for analysis. Since the loops are highly flexible, they were omitted from the PCA analysis. We used the g_covar module of GROMACS to perform the PCA and to extract eigenvalues and eigenvectors. To investigate the crucial dominant motions, conformational changes along the two principal components (PC1 and PC2) were analyzed.

Results

The analysis presented in the results section was done by assembling all 10 trajectories for each system into an ensemble. Each of the 10 trajectories was run for 100 ns.

How Do the Mutations Stabilize the Active and Inactive States of A2AR

The crystal structures of both the antagonist-bound and the agonist-bound A2AR are available for both the wild type and thermostabilized receptors. To investigate the intrinsic dynamic behavior and the energetics of the mutants compared to the wild type A2AR, we performed MD simulation studies of 1 μs each in the following receptor systems: the crystal structures of inactive state StaR2 and active-like state GL31 thermostable mutants with and without their respective ligands and the wild type receptor with and without ligands. The wild type receptor structures in the inactive and active-like conformations were derived by mutating the crystal structure of the thermostable mutants to the wild type sequence. Details of all the simulation systems are shown in Table S1 (Supporting Information).

Stabilization of Both GL31 and StaR2 Mutants Comes from Enthalpic Contributions

Figure 1 shows the calculated enthalpy averaged over the MD trajectories of the thermostable mutants and the wild type receptor in both the inactive and active-like state conformations. While the StaR2 mutant (shown in red in Figure 1) is stable in its inactive state compared to the wild type in the inactive state, the GL31 mutant is stable in its active-like state compared to the wild type in the active-like conformation. This finding is in contrast to our previous calculations on the thermostable mutant of the inactive state of the β1-adrenergic receptor that showed very little enthalpic stabilization of the mutant compared to the wild type.[36] MD simulations starting from the active-like conformation of the StaR2 amino acid sequence showed a collapse of the active-like structure to the inactive state within 1 μs, implying the thermostabilizing mutations in StaR2 bias the conformation to the inactive state (see Figure S2, Supporting Information). These results show that the mutations in the thermostable mutants stabilize specific receptor conformations.

Figure 1

(A) Calculated enthalpy of the various thermostable mutant and wild type sequences in both the inactive (left) and active-like (right) conformations without ligands. The energies of the GL31 thermostable active-like mutant are shown in black, the inactive state mutant StaR2 in red, and the wild type in blue. (B) The free energy change upon mutations compared to the difference in the measured experimental stability of StaR2. (C) The measured Tm values compared to the calculated free energy difference for the wild type and thermostable mutants for GL31.

Entropic Contribution to Thermostability Is Less Significant than Enthalpy

We calculated the second order entropy as described in the Methods section, and observed that the entropic contribution to stability was insignificant (the TΔS factor is about 20 times less than the enthalpy) compared to the enthalpic contribution (see Figure S3, Supporting Information). It is interesting to note that the enthalpic contribution to stability in the β1-adrenergic receptor was not substantially different for the thermostable mutants compared to wild type.

Calculated free energy changes for single point mutations in GL31 and StaR2 correlate with experimentally measured thermostability: To understand the thermodynamic consequences of single point mutations, we have used the alchemical free energy simulations - thermodynamic integration (TI) method (detailed in the Methods section) to calculate the free energy change upon single point mutations and compare them to the experimental stabilities. The positions of mutations in GL31 are L48A2.46, A54L2.52, T65A2.63, and Q89A3.37, and in StaR2, they are A54L2.52, T88A3.36, R107A3.55, K122A4.43, L202A5.63, L235A6.37, V239A6.41, and S277A7.42. The results for the convergence tests of each of these simulations are described in the Supporting Information (Figure S4). Figure 1B shows the quantitative correlation of the calculated change in free energy to the experimental measured stabilities measured by a single-point binding assay using [3H]-ZM241385.[13] There is significant correlation between the calculated free energy difference and the measured stabilities except for the S277A mutation. T88A with the highest measured thermostability has a 28.2 kcal/mol difference in free energy. The structural basis for this stability by T88A is discussed in the next section. Figure 1C shows the correlation of the calculated change in free energies for single mutants and experimental stability ΔTm. Here Tm is the temperature at which 50% of the solubilized receptor can still bind the radiolabeled NECA (agonist) after heating for 30 min.[16] L48A mutation shows the highest thermostabilization (with ΔTm = +14 °C) and also has the largest change in free energy upon mutation. We observe a good correlation between the calculated free energies and the experimentally measured stability. Details of the convergence of the free energy calculations for point mutations are shown in Figure S4 (Supporting Information).

Interhelical Hydrogen Bond and van der Waals Packing Increase the Stability of the Thermostable Mutants

Energetic contributions from interhelical hydrogen bonds and hydrophobic interactions are important in determining the stability and the packing of the TM core in GPCRs.[37−41] Therefore, to analyze the structural basis of the enthalpic stabilization of the thermostable mutants shown in the previous section, we calculated the difference in the number of interhelical hydrogen bond and van der Waals (vdW) interactions that are stably formed during the dynamics of the thermostable mutants GL31 and StaR2 and their respective wild type conformations.

Interhelical Hydrogen Bond Interactions

We calculated all possible interhelical hydrogen bonds (backbone–backbone, backbone–side chain, and side chain–side chain) formed by A2AR in both GL31 and StaR2 mutants and the wild type that show a significant population (see the Methods section) over the course of the MD simulations.

Figure 2 shows the total number of interhelical hydrogen bond and vdW interactions where each helix is represented as a circle. The total number of interhelical contacts for each system is shown below each figure. The numbers of interhelical hydrogen bonds and vdW contacts are higher in both thermostable mutants GL31 and StaR2 compared to their respective wild type conformations. This makes both StaR2 and GL31 more stable enthalpically than the wild type in the respective states. Both GL31 and the corresponding wild type in the active-like conformation show a greater number of interhelical hydrogen bonds than their corresponding inactive state structures. This is especially interesting, since the number of interhelical hydrogen bonds is higher in the intracellular region of TM5–TM6 and TM6–TM7 where the G protein couples with the receptor (8 hydrogen bonds in GL31 and 4 hydrogen bonds in StaR2). Interestingly, the centrally located TM3 helix shows a greater number of interhelical hydrogen bonds with other helices in the inactive StaR2 (10 H-bonds) mutant showing tighter helical core packing in the inactive state. In GL31, we observed better packing between TM5 and TM6 compared to StaR2, which is known to move as a unit with respect to TM3 upon activation.

Figure 2

Interhelical interaction networks for active-like and inactive mutants (GL31 and StaR2) and wild type receptor without ligands. Seven transmembrane α-helices are shown as circles. The black lines with different thicknesses show the interhelical interactions between each pair of transmembrane helices, and the number of such interactions is shown on the lines. Parts A–D show the interhelical hydrogen bond interaction, while parts E–H indicate the interhelical van der Waals (vdW) interactions.

Interhelical van der Waals Interaction

Mutation of hydrophobic residues in the interhelical interface has shown that these contacts contribute significantly to the thermal stabilization of receptors.[42,43] Figure 2E–H shows that both of the thermostable mutants have a greater number of favorable interhelical vdW interactions compared to the wild type. In contrast to the interhelical hydrogen bonds, the inactive state mutant StaR2 and the corresponding wild type in the inactive state show more interhelical vdW interactions than the active-like state. The vdW packing centered around TM3 and also between TM3 and TM5 is stronger in the inactive state mutant than in the active-like GL31 (StaR2, 15; GL31, 11). In the GL31 active-like state, TM5 interacts more strongly with TM6 and weakly with TM3. This could facilitate the outward motion of TM6 away from TM3 in the active-like conformation. In summary, the number of interhelical interactions is higher in the thermostable mutants compared to the wild type receptor, thus accounting for the enthalpic stabilization of the receptor mutants.

Single Point Mutations Have a Crucial Structural Role in Stabilizing the Active-Like State of A2AR

In this section, we have analyzed the structural basis of the thermostability due to single point mutations that specifically favor stabilization of the active-like state GL31 compared to the inactive state. The four mutations in GL31 are clustered on TM2 and TM3: L48A2.46, A54L2.52, T65A2.63, and Q89A3.37. L48A2.46 and Q89A3.37 are mutations that selectively stabilize the agonist-bound state of A2AR.[14]

L48A2.46 shows a marked 14 °C increase in thermostability specifically to the agonist-bound receptor.[21] We observed that the mutation of L48A led to the formation of the interhelical hydrogen bond between the backbone amide nitrogen of L48A and the side chain of S943.42, as shown in Figure 3, and this hydrogen bond is not observed in the crystal structure of GL31.

Figure 3

Interhelical interactions of the agonist specific mutations L48A in active-like GL31. (A) Population of the interhelical hydrogen bond between the amide nitrogen of L48A mutant and the side chain of S94 in the GL31 (black curve) and its wild type (red curve). Three inset structures represent the three maxima in the population density: (I) GL31, (II) one possible conformation of the wild type receptor within the hydrogen bond, and (III) the second populated conformation of the wild type where this hydrogen bond is broken. (B) The inter-relationship between the formation of the hydrogen bond and movement of the intracellular region of TM6 away from TM3. (C) The population distribution of the receptor conformations (GL31 in black and wild type in red) that have a TM3–TM6 distance greater than 7.9 Å with the hydrogen bond distance between L48 and S94.

This interhelical hydrogen bond is not present in the wild type receptor due to the steric hindrance from the side chain of L482.46 located between the TM2 and TM3 helices (Figure 3A, insets II and III). As seen in Figure 3A, this hydrogen bond is well populated in the GL31 mutant, while it is weakly populated in the wild type.

We further examined if the formation of this hydrogen bond is correlated to the stabilization of the active-like state of the receptor. We observed that formation of the interhelical hydrogen bond between the backbone of L48A2.46 and the side chain of S943.42 (these two residues are located in the lower half of their respective helices) correlates with the increase in distance between the intracellular half of TM3 and TM6, as shown in Figure 3B. The distance between the backbone atoms of the last pair of residues in the intracellular half of TM3 and TM6 (R102–A231) is shown to increase above 7.8 Å (which is the distance in the crystal structure of GL31) when the probability of a hydrogen bond between TM2 and TM3 is high. Thus, when the intracellular part of TM3 is pulled toward the backbone of TM2, TM6 is free to move away from TM3 and hence this mutation facilitates the stabilization of the active-like state. This is illustrated in Figure 3C, which shows the population distribution of the conformations that have the TM3–TM6 distance above 7.8 Å as a function of the hydrogen bond distance between L48A2.46 and S943.42. Figure 3C shows that when the hydrogen bond is formed between L48A2.46 and S943.42 the population of receptor conformations of GL31 (shown in black bars) showing larger movement of TM6 is higher than wild type (red bars).

The Q89A3.37 mutation shows preferential stabilization of the active-like state by 6 °C but also shows destabilization of the inactive state by 8 °C.[21] Mutation of the large Q893.37 residue to a smaller Ala residue also results in formation of a backbone side chain hydrogen bond between its neighbor V863.34 of TM3 and S1324.53 on TM4, as seen from Figure S5 in the Supporting Information.

This specific hydrogen bond called Cα—H···O=C type, which is weaker than the N—H···O hydrogen bond, has been observed in transmembrane helical proteins.[44−48] In addition, the neighboring amphiphilic residue T883.36 shows interhelical vdW interaction with W2466.48 of TM6 in GL31 which is insignificant in the wild type. This mutation also leads to rearrangement in the intracellular part of TM3 and TM6 with reduced vdW interactions between I1063.54 and the aliphatic chain of K2276.29. This weakened interaction could release TM6 to move into the active-like state (bottom part of Figure S5, Supporting Information). The population density of the hydrogen bond V86–S132 and the vdW interaction T88–W246 in both the mutant and wild type are shown in Figure S6 of the Supporting Information.

Structural Basis for Stability of Mutations in the Inactive State StaR2

There are eight mutations in StaR2 spread over TM2 to TM7: A54L2.52, T88A3.36, R107A3.55, K122A4.43, L202A5.63, L235A6.37, V239A6.41, and S277A7.42. Here, A54L2.52 is common to both inactive StaR2 and active-like GL31 mutants.[15] The mutations R107A3.55 and L202A5.63 are proximal and improve the vdW packing between I2005.61, I1063.54, and A2045.65 on the intracellular regions of TM3 and TM5 (Figure 4A). L2025.63 is close to the conserved residue Y1975.58 (TM5), which in turn makes contact with the backbone oxygen atom of L235A6.37 (Figure 4B). This hydrogen bond is not possible in the wild type due to the long side chain of L202 and L235 being in the way of this interhelical hydrogen bond. Mutation of L202A5.63 and V239A6.41 also improves the interhelical vdW packing between TM3, TM5, and TM6. Thus, while bulkier hydrophobic side chains provide good interhelical contact, mutation of such residues to Ala provides a tighter packing of the helical backbone by facilitating interhelical backbone–side chain hydrogen bonds and improved van der Waals packing of the neighboring residues. K1224.43 is in the intracellular region of TM4 facing the lipid. Mutation of this residue to Ala favors better packing of the side chain of the neighboring I1244.45 with F442.42 (brown arrow, Figure S7, Supporting Information).

Figure 4

Interhelical packing comparison between inactive StaR2 (right) and wild type (left) in the TM3, TM5, and TM6 transmembrane helices. (A) Enhanced vdW packing between residue pairs (shown in sticks) that are near the mutation positions R107A3.55 and L202A5.63 (shown in spheres). (B) The interhelical hydrogen bond between conserved Y1975.58 on TM5 and backbone of L235A6.37 mutant on TM6 shown in dotted red line. Residues that are involved in interhelical interaction are shown as cyan sticks, while the vdW packing is shown as surface and the hydrogen bond as red dotted line.

The T88A3.36 mutation leads to a significant change in measured stability. Two main clusters of hydrophobic interactions near T88A3.36 and S277A7.42 mutants facilitate strong interhelical hydrophobic interactions between TM2, TM3, TM6, and TM7. T88A3.36 and its neighboring residues V843.32 and L853.33 interact with F622.60, F2426.44, and W2466.48 (orange arrows, Figure 5). S277A7.42 and its spatial neighbor T88A3.36 tighten the packing between W2466.48 and L2496.51 on TM6 and L2767.41 on TM7 that is close to S277A7.42 mutant (cyan arrows). The distance distributions of the residue pairs that show enhanced vdW interactions are shown for the wild type and the mutants in Figure S8 of the Supporting Information.

Figure 5

Enhanced interhelical vdW packing of residue pairs (shown in sticks) due to mutations T88A3.36 and S277A7.42 mutants in the inactive mutant StaR2 (right) and its wild type (left). The residue pairs with enhanced vdW interactions close to T88A3.36 and S277A7.42 mutants are highlighted in the orange and cyan surface and double ended arrows, respectively, and the positions of mutations are shown as spheres.

Difference in Activation Mechanism of the Thermostable Mutants and Wild Type Receptor

The crystal structures of the inactive and active-like thermostable mutants of A2AR show movement of TM5 and TM6 away from TM3 upon binding of the agonist, adenosine. We calculated the dynamic range of the movement of TM5 and TM6 away from TM3 for the thermostable mutants and the wild type receptors.

A salt bridge is formed between the side chains of Arg1023.50 (that is highly conserved in class A GPCRs) and Glu2286.30. This salt bridge, also called the “ionic lock” shown in Figure 6A, has been observed in the inactive state structures of rhodopsin,[49] but it appears to be more dynamic in other class A GPCRs and may be broken even when an antagonist is bound.[50] Figure 6 shows the population density of the “ionic lock” distance that characterizes the extent of movement of TM6 away from TM3. The ionic lock is formed in the inactive StaR2 crystal structure with the antagonist ZM241385 bound and not in the crystal structure of the active-like GL31 structure with the agonist adenosine bound. In our MD simulations, the ionic lock is dynamic in the inactive states of both StaR2 and the wild type A2AR. There is a higher population of the ensemble forming a strong ionic lock both in the inactive state StaR2 mutant and the wild type in the inactive state (red lines in Figure 6B) than in GL31 mutant and the wild type receptor in the active-like state (black lines in Figure 6B). However, the Cα distances between the residues that make the ionic lock are distributed uniformly about the distances in the crystal structure of both StaR2 and GL31 (shown in Figure 6C).

Figure 6

The salt-bridge interaction known as the “ionic lock” between Arg1023.50 on TM3 and Glu2286.30 on TM6 in the inactive mutant state. (A) Superposition of X-ray crystal structures of active-like GL31 (cyan) and inactive StaR2 (orange). Active-like GL31 shows the absence of the ionic lock owing to side chain rotation of E2286.30 and outward movement of TM6. (B) Population density of the ionic lock between the cationic guanidinium group (NH*) of R1023.50 and the anionic carboxylate (OE*) of E2286.30 in GL31 dynamics with adenosine bound (solid black lines), StaR2 dynamics with antagonist ZM241385 bound (solid red lines), and wild type receptors in the inactive (with antagonist ZM241385 bound, dashed red lines) and active-like states (with adenosine bound, dashed black lines). (C) The same color scheme as in part B showing the distribution of the Cα–Cα distance of R1023.50 on TM3 and E2286.30 on TM6. In part B and C, red and black vertical dotted lines show the Cα–Cα distances observed in the respective crystal structures.

Wild Type A2AR Is More Dynamic than the Thermostable Active-Like or the Inactive State Mutant

To analyze the extent of structural dynamics shown by the wild type A2AR compared to the thermostable mutants, we calculated the number of microstates in the most populated conformational state from the MD simulation trajectories. This was done by clustering analysis based on root-mean-square deviation in coordinates (see the Methods section for details). Figure S9A of the Supporting Information shows the population densities of various conformational states sampled by the thermostable mutants and the wild type receptors. Figure S9B (Supporting Information) shows the number of microstates within the most populated conformation in each system. It is seen that the number of microstates for the wild type receptor is higher than both GL31 and StaR2 thermostable mutant receptors, showing more flexibility than the thermostable mutants.

Residues with High Stress Are Required for Activation of the Receptor

The distribution of net internal force on each residue, known as stress, provides insight into the regions of high stress in the wild type A2AR and the thermostable mutants, as we had previously shown for the thermostable mutant of the β1-adrenergic receptor.[36] We hypothesize that high stress at the location of functionally important residues leads to large scale conformational changes required to activate the receptor. To examine this, we have calculated the stress on each residue using both the bonded and nonbonded components,[35] averaged over the MD simulations for the thermostable mutants and wild type A2AR (Figure S10, Supporting Information). Both the inactive state and the active-like state of the wild type A2AR show higher stress than their respective thermostable mutants StaR2 and GL31. Figure 7 shows the positions of high stress calculated for the wild type A2A receptor in the active-like state. The residues shown in pink sticks in Figure 7 are positions of high stress. We also calculated the root-mean-square deviations (RMSDs) in coordinates of the corresponding residues in the fully active G-protein-coupled state of the β2-adrenergic receptor (pdb ID: 3SN6) and the inactive state of the β2-adrenergic receptor (pdb ID: 2RH1).

Figure 7

Residue positions of high stress in the wild type A2A receptor show maximum movement upon activation. Blue to red color coding is the difference between the coordinates of the inactive state to the active state receptor.

The RMSD value ranging from low to high, blue to red in Figure 7, shows the extent of movement of each residue upon activation. The residues such as T2797.44 and N2807.45 on TM6 and TM7 show large fluctuations upon activation, as well as the high stress. Thus, we observe that residue positions of high stress are required for movement of helices during activation and reducing this stress could stabilize the receptor but could make the receptor inefficient in G protein coupling or activation.

Effect of Ligand Binding on the Stability of GL31 and StaR2

Tate and co-workers observed that the stabilization of the thermostable mutants by ligand binding (either adenosine or ZM241385) is required for crystallization of the mutant A2A receptors.[20] To understand the effect of ligand binding on the thermostability, we calculated the enthalpic contributions due to ligand binding for active-like GL31 and inactive StaR2 structures and wild type A2AR in both the inactive and active-like conformational states (Figure S11, Supporting Information). We observed that the binding of the agonist adenosine stabilizes the active-like state mutant GL31 more than the stabilization resulting from the binding of the antagonist ZM241385 to the inactive state StaR2 mutant. This is consistent with the experimental observation that agonist binding was required to purify the GL31 mutant, whereas the StaR2 mutant could be purified in the absence of antagonist. Also, ligand binding leads to greater stability increase for the thermostable mutants compared to the wild type. Adenosine binding stabilizes the GL31 mutant better than the wild type by 19 kcal/mol. Similar stabilization of the thermostable mutant StaR2 (16 kcal/mol) was observed relative to the wild type receptor by the antagonist ZM241385 binding.

The Dynamics of the Ligand in the Thermostable Mutants

Crystal structures show partial overlap of the binding sites of the agonist adenosine and the antagonist ZM241385 with some features that are distinct to each ligand. In ZM241385, the furan group (O25) forms a hydrogen bond with N2536.55 in TM6, but the ribose ring moiety in the adenosine agonist forms hydrogen bonds with S2777.42 and H2787.43 in TM7 (Figure S12A, Supporting Information).

We have examined the dynamics of the ligand receptor interaction during the MD simulations. The populations of the hydrogen bonds between adenosine ligand and S2777.42 and H2787.43 in GL31 were not stable during the dynamics. Instead, T883.36 on TM3 and N1815.42 on TM5 showed direct hydrogen bonds with adenosine (Figure S12B, Supporting Information). However, the ligand receptor interactions were well maintained during the dynamics of the StaR2 mutant. The Asn2536.55 within the binding pocket remained connected to the ZM241385 ligand directly, as observed in the crystal structure of the StaR2 mutant. Analysis of the number of water molecules in the ligand binding pocket showed that StaR2 has a few more waters than active-like GL31 (average number of 10 in GL31 and 11 in StaR2 (Figure S12C, Supporting Information)). While GL31 and the wild type A2AR in the active-like conformation have similar water densities around the ligand, StaR2 mutant has less water in the ligand binding site than the wild type in the inactive state (average waters 11 and 13 in StaR2 and wild type in the inactive state). Both adenosine and ZM241385 retain direct contact with the residues in their respective binding sites in the mutants compared to the wild type, showing that the ligand stabilizes the mutants more than the wild type.

Discussion and Conclusion

Using long time scale molecular dynamics simulations and the thermodynamic integration method for calculating the free energy change upon mutation, we have provided insights into the distinct structural and energetic characteristics that contribute to the stability of the active-like (GL31) and inactive state (StaR2) thermostable mutants of A2AR. We have shown that the wild type receptor is less stable even when ligand is bound compared to both of the thermostable mutants GL31 and StaR2. While the StaR2 mutant is stable in the inactive state, GL31 is more stable in its active-like state than its inactive state. MD simulations starting from the amino acid sequence of StaR2 in the active-like conformation collapse to the inactive state, showing that the StaR2 sequence is optimized to stabilize the inactive receptor conformation. Enthalpic contributions to the free energy stabilize the thermostable mutants of A2AR. The entropic contributions to stabilization are low, as reflected by the entropy calculated using mutual information. Additionally, the number of microstates present in the ensemble of the thermostable mutants is less than that in wild type. This is in contrast to the thermostable mutant m23 of the β1-adrenergic receptor, which is stabilized by increasing the side chain entropy.[36]

We observed the stress calculated on each residue is less in the thermostable mutants than in the corresponding conformations of the wild type receptor. There are residue positions with high stress (net force) in both the wild type A2AR and β1-adrenergic receptor receptors. This stress may be essential to cause the receptor movement upon activation. However, in both the β1-adrenergic receptor and A2AR systems, the stress on each residue is reduced upon making the thermostable mutations. We observed that the high points of stress in the receptors are very often the most conserved positions such as Pro5.50, Pro6.50, and Pro7.50. These proline residues play an important part in activation of the receptor by enabling the movement of the helices. Reducing the stress at these high stress points by mutating neighboring residues could also reduce the potential for the receptor to get activated which is possibly why the thermostable mutants show markedly reduced G protein activation.[36] We observed that both GL31 and StaR2 are less dynamic than their corresponding wild type conformations. The movement of TM6 with respect to TM3 in GL31 is more restrictive than in the wild type, as represented in Figure 8. While this restricted movement is advantageous in purification and crystallization of GL31 mutant receptor, it could be the reason why these mutants are not efficient in activating the G protein.[16]

Figure 8

Schematic representation of the extent of dynamic motion observed in the MD simulations for the wild type (left), and GL31 mutant (right). TM3, TM5, and TM6 helices are shown in blue.

The alchemical free energy changes calculated for single point mutants correlate well with the measured stabilities. Most of the thermostabilizing mutations in GL31 and StaR2 are mutating large residues such as Phe that show well-packed side chain conformations, to Ala. We observed that, in large to small side chain mutations, the loss of side chain packing is compensated by the main chain of the TM helices forming a hydrogen bond with the side chain, or main chain of residues in neighboring helices, thus providing a stronger packing interaction. Mutations also lead to rearrangement in the side chain conformations of nearby residues that improve hydrophobic packing with neighboring helices. The two mutations that specifically stabilize the active-like state (L48A2.46 and Q89A3.37) show a correlation of an interhelical hydrogen bond formed as a result of these mutations to the opening of TM3 and TM6 in the mutant that was not observed in the dynamics of the wild type.