Advertisement
Browse Subject Areas
?

Click through the PLOS taxonomy to find articles in your field.

For more information about PLOS Subject Areas, click here.

  • Loading metrics

Open Access

Peer-reviewed

Research Article

Molecular Basis of Ligand Dissociation in β-Adrenergic Receptors

  • Angel González,

    Affiliations Laboratori de Medicina Computacional, Unitat de Bioestadística, Facultat de Medicina, Universitat Autònoma de Barcelona, Bellaterra, Catalunya, Spain, Departamento de Ciencias Biológicas, Facultad de Ciencias Biológicas, Universidad Andrés Bello, Santiago, Chile

    ⨯
  • Tomas Perez-Acle,

    Affiliations Computational Biology Lab, Center for Mathematical Modeling, Facultad de Ciencias Físicas y Matemáticas, Universidad de Chile, Santiago, Chile, Centro Interdisciplinario de Neurociencias de Valparaíso, Playa Ancha, Valparaíso, Chile, Fundación Ciencia para la Vida, Ñuñoa, Santiago, Chile

    ⨯
  • Leonardo Pardo,

    Affiliation Laboratori de Medicina Computacional, Unitat de Bioestadística, Facultat de Medicina, Universitat Autònoma de Barcelona, Bellaterra, Catalunya, Spain

    ⨯
  • Xavier Deupi

    xavier.deupi@psi.ch

    Affiliation Condensed Matter Theory Group and Laboratory of Biomolecular Research, Paul Scherrer Institut, Villigen PSI, Switzerland

    ⨯

Abstract

The important and diverse biological functions of β-adrenergic receptors (βARs) have promoted the search for compounds to stimulate or inhibit their activity. In this regard, unraveling the molecular basis of ligand binding/unbinding events is essential to understand the pharmacological properties of these G protein-coupled receptors. In this study, we use the steered molecular dynamics simulation method to describe, in atomic detail, the unbinding process of two inverse agonists, which have been recently co-crystallized with β1 and β2ARs subtypes, along four different channels. Our results indicate that this type of compounds likely accesses the orthosteric binding site of βARs from the extracellular water environment. Importantly, reconstruction of forces and energies from the simulations of the dissociation process suggests, for the first time, the presence of secondary binding sites located in the extracellular loops 2 and 3 and transmembrane helix 7, where ligands are transiently retained by electrostatic and Van der Waals interactions. Comparison of the residues that form these new transient allosteric binding sites in both βARs subtypes reveals the importance of non-conserved electrostatic interactions as well as conserved aromatic contacts in the early steps of the binding process.

Citation: González A, Perez-Acle T, Pardo L, Deupi X (2011) Molecular Basis of Ligand Dissociation in β-Adrenergic Receptors. PLoS ONE 6(9): e23815. https://doi.org/10.1371/journal.pone.0023815

Editor: Hendrik W. van Veen, University of Cambridge, United Kingdom

Received: June 3, 2011; Accepted: July 25, 2011; Published: September 7, 2011

Copyright: © 2011 González et al. This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

Funding: This work was supported by the Swiss National Science Foundation (SNF) grant 31003A_132815 and the ETH Zürich within the framework of the National Center for Competence in Research Structural Biology Program (to XD); by grants from Ministerio de Ciencia e Innovación (MICINN) SAF2010-22198-C02-02 and Instituto de Salud Carlos III (ISCIII) RD07/0067/0008 (to LP); and by the Proyectos de Financiamiento Basal PFB03 and PFB16 grants and the Millenium Institute Centro Interdisciplinario de Neurociencias de Valparaíso (CINV) (to TPA). AG is recipient of a Comisión Nacional de Investigación Científica y Tecnólogica de Chile (CONICYT) doctoral fellowship. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Competing interests: The authors have declared that no competing interests exist.

Introduction

G protein-coupled receptors (GPCRs) represent one of the largest protein families in mammals [1] and constitute 2%–3% of the human proteome [2]. GPCRs transduce sensory signals of external origin, such as photons, odors or pheromones, and endogenous signals including biogenic amines, (neuro)peptides, proteases, glycoprotein hormones and ions, into the cell. Thus, these receptors are essential in cell physiology, and their malfunction is commonly translated into pathological outcomes [3]. As a result, GPCRs constitute one of the most important pharmaceutical targets, as around 40% of prescribed drugs act through this family of proteins [4]. These receptors feature a common fold of seven transmembrane helices (TMs 1 to 7) connected by three extracellular (ECLs 1 to 3) and three intracellular (ICLs 1 to 3) loops [5], with an extracellular N-terminus and an intracellular C-terminus. Extracellular regions are very diverse in structure and amino acid composition, and in many GPCRs, as glycoprotein hormone and peptide receptors in family A or most receptors in families B and C, they are directly involved in ligand binding [6]. While smaller ligands bind in a pocket relatively buried within the TM bundle, they must also interact with the extracellular regions in order to reach the binding site. Understanding the molecular basis of ligand-receptor interactions in the extracellular domains is of great importance, as they are implicated in many aspects of receptor function, as ligand binding [7] and specificity [8], allosterism [9] or receptor activation [10], [11]. Importantly, recent NMR data show ligand-specific conformational changes in the extracellular surface of the β2-adrenergic receptor (β2AR) [12].

While there is a vast amount of pharmacological, functional and pathophysiological information about GPCRs deposited in specialized databases (e.g. IUPHAR-DB, at http://www.iuphar-db.org), structural data of GPCRs is still scarce. Presently, only the structures of eight Class A GPCRs (bovine and squid rhodopsins, human β2-adrenergic, turkey β1-adrenergic, human A2A adenosine (reviewed in [13], [14], [15]), human CXCR4 chemokine [16], human dopamine D3 [17] and human histamine H1 [18] receptors) are known. The availability of the structure of the β1AR [19] and β2AR [20] represents a unique opportunity to investigate the similarities and/or differences in the ligand entry process between these closely related subtypes. While these receptors have slightly different pharmacological properties [21], they present a strong similarity in sequence and structure, particularly in the TM bundle and orthosteric binding pockets [19]. Thus, it is plausible to argue that extracellular regions can have an impact on the different pharmacological properties between subtypes. Previous theoretical studies, using random acceleration molecular dynamics simulations, have suggested that ligands access the orthosteric binding site of the β2AR mainly through an opening at the extracellular surface [22]. Conversely, ligand docking calculations in opsin located the paths for access/egress between transmembrane helices [23]. This difference is due to both the different architecture of the extracellular regions and the different chemical nature of their respective ligands. While the β2AR binding pocket is relatively exposed to the solvent, ECL 2 and the N-terminal of opsin cover the binding pocket, which form a “plug” that prevents the access of the ligand from the extracellular environment.

In this work, we have conducted a comparative analysis of the process of ligand dissociation in β1 and β2ARs using the steered molecular dynamics (SMD) simulation method [24]. SMD has been very successful in the study of dissociation reactions of several small-molecules/protein complexes through application of external forces on nanosecond time scales [25], [26], [27], [28], and is particularly useful to describe the interactions occurring in the binding/unbinding of ligands [25]. Our results suggest that both receptors have two putative ligand entry channels located at the extracellular region, discarding the entry channels located between the transmembrane segments that lead to the lipid environment. By monitoring the forces and energies of the ligand-dissociation along these extracellular channels in both βAR structures, we have identified for the first time two secondary binding pockets in the extracellular region of the receptors. In addition, we discuss the importance for the ligand exit/entry process of non-conserved charged residues and conserved aromatic interactions shared by the two entry channels.

Results

Ligand entry/exit channels in β1 and β2 adrenergic receptors

Using the skeleton search algorithms implemented in the CAVER program [29], we explored routes that connect the buried orthosteric binding pocket to the extracellular surface in the structures of the human β1AR and β2AR. Figure 1 displays two entry channels identified in each receptor, located between TMs 3, 5, 6 and 7 (C1) and TMs 1, 2, 3 and 7 (C2). These channels are separated from each other by charged residues in ECLs 2 and 3; D217 and D356 in β1AR (Figure 1a) and D192 and K305, forming a salt bridge, in β2AR (Figure 1b). These residues, in combination with other neighboring polar/charged amino acids, confer a negative electrostatic potential to both channels, which suggests the existence of an electronegative funnel to attract positively charged ligands into the orthosteric binding site of beta adrenoceptors [30]. On the other hand, the entrance/exit channels for retinal in rhodopsin have been proposed to occur through the lipid bilayer, via two openings located between TMs 1 and 7, and TMs 5 and 6, respectively [23]. While CAVER does not detect these alternative channels in the βAR structures, in order to further assess their possible relevance, we identified these two channels on the structure of the ligand-free apoprotein opsin (PDB entry 3CAP [31]) and mapped them onto the βAR structures by coordinate superimposition (C3 and C4 in Figure 1c).

thumbnail
Download:
Figure 1. Extracellular molecular surfaces of the β1AR (panel a) and β2AR (panel b), embedded in a lipid bilayer (in yellow).

The electrostatic potential was calculated using the program APBS with nonlinear Poisson-Boltzmann equation and contoured at ±10 kT/e (negatively and positively charged surface areas in red and blue, respectively). The accessible channels (C1 and C2) identified by CAVER are depicted as green wires. D217/D356 in β1AR (panel a) and the salt bridge D192/K305 in β2AR (panel b) are represented by circled − and + symbols. Panel c displays the extraction vectors along the four channels (C1 to C4) at the end of the equilibration run. The β2AR ribbon structure is colored as follows: TM1 (grey), TM2 (yellow), TM3 (red), TM4 (black), TM5 (green), TM6 (blue), TM7 (cyan), and helix 8 (red). Carazolol is shown in white sticks. Pictures were prepared using PyMOL (http://www.pymol.org/).

https://doi.org/10.1371/journal.pone.0023815.g001

Channel route preferences for ligand dissociation

To study the process of ligand release from β1AR and β2AR orthosteric binding pockets, we performed SMD simulations of the antagonist-receptor complexes embedded in a model lipid bilayer (see Methods). Ten nanoseconds of equilibration were performed to obtain constant values of energy, cell volume and lipid density. The root mean square deviation (rmsd) of the protein backbone atoms from the initial coordinates during equilibration stabilizes rapidly to a value in the vicinity of 2.0 Å (Figure S1). Following this equilibration period, steered forces were applied to both ligands along the four calculated channels (C1 to C4 in Figure 1c). Figure S2 displays representative force profiles of the pulling experiments of cyanopindolol (Figure 2a) and carazolol (Figure 2b) along extracellular C1 (black) and C2 (blue) and lipid C3 (red) and C4 (yellow) channels. The initial force peaks to remove ligands from the orthosteric binding site via extracellular C1 or C2 channels were on average ∼600 pN, a typical value in ligand diffusion SMD experiments [26], [32]. On the contrary, pulling the ligands through the proposed rhodopsin channels (C3 and C4 in Figure 1c), required forces two-fold larger than for the extracellular routes. These results strongly suggest that, in βARs, the transit of molecules through the lipidic phase, via TMs 1 and 7 or TMs 5 and 6, is not favored compared to the extracellular routes. Consequently, the C3 and C4 channels were not included in the rest of the analysis.

thumbnail
Download:
Figure 2. PMF and force profiles of ligand extraction along the extracellular channel C1.

Panel a, dissociation of cyanopindolol from β1AR. Panel b, dissociation of carazolol from β2AR. The star symbols correspond to the snapshots depicted in Figure 5. The statistical error in the PMF data is shown in bars. Inset figures display representative force profiles of the repeated trajectories. The force simulation data is shown in grey and smoothed to a black line. Horizontal bars denote regions with positive slope in the force profile.

https://doi.org/10.1371/journal.pone.0023815.g002

Residues implicated in ligand-receptor interactions during dissociation

Figures 2 and 3 display the potential of mean force (PMF) and representative force profiles (insets) of the pulling experiments of cyanopindolol (Figures 2a and 3a) and carazolol (Figures 2b and 3b) along the extracellular C1 and C2 channels. In all cases, small fluctuations were observed in receptor structures during ligand extraction, which were in similar ranges to the rmsd values of the equilibration runs (data not shown). These results indicate that selected velocities, force constants, and extraction vectors were adequate to achieve smooth ligand releases. Thus, no steric clashes occur between molecules and receptors during dissociation. Horizontal bars in the insets of Figures 2 and 3 represent time periods of relatively strong ligand-receptor interaction during dissociation. Positive slopes in force profiles characterized these periods. Clearly, disruption of the initial interactions between the ligands and orthosteric binding site residues, which mainly include the electrostatic interaction with D3.32 and hydrogen bonds with N7.39 and S5.42 (superscript numbers correspond to the Ballesteros & Weinstein general numbering scheme [33]), requires a maximal force (ramp symbols in the force insets). After this primary unbinding event (∼0.5 ns), forces fall as the ligands displace further from the orthosteric binding site towards the solvent through the exit channels. Then, subsequent regions of increasing forces indicate secondary interaction sites along the channels.

thumbnail
Download:
Figure 3. PMF and force profiles of ligand extraction along the extracellular channel C2.

Panel a, dissociation of cyanopindolol from β1AR. Panel b, dissociation of carazolol from β2AR. The star symbols correspond to snapshots depicted in Figures 6. See legend of Figure 2 for details.

https://doi.org/10.1371/journal.pone.0023815.g003

The extraction of cyanopindolol through channel C1 in β1AR reveals two major retention events, (Figure 2a). In an initial step at ∼1.2 ns, cyanopindolol is stabilized by an ionic interaction between the protonated amine of the ligand and D217 in ECL2 and a hydrogen bonding interaction between the β-OH group and D3567.32 (Figure 4a). Later, in the final steps of its movement toward the extracellular solvent (∼2.0 ns), increasing forces are required to break a salt bridge between E205 in ECL2 and R351 in ECL3 (also shown in Figure 4a), in order to allow the ligand escape. Conversely, the extraction of carazolol from β2AR through C1 is characterized by a single retention site at ∼1.2 ns (Figure 2b). At this point, the protonated amine of the ligand interacts with D192 in ECL2 and the β-OH group with N301 in ECL3 (Figure 4b). Table 1 lists residues in the vicinity of the ligands during the dissociation process that form this extraction channel. On the other hand, ligands extraction through channel C2 shows two retention sites at ∼0.9 and ∼1.5 ns in both adrenoceptors (Figure 3). Initially, the protonated amine of cyanopindolol or carazolol interacts with either D217 or D192 in ECL2 of β1- and β2- receptors, respectively (Figures 5a and 5b). The second barrier corresponds to Van der Waals attractive forces between the aromatic moieties of the ligands and bulky residues at positions 2.64, 2.65, 3.28, 7.36, 7.39 and 7.40 in TMs 2, 3 and 7 (summarized in Table 1). In the final steps of the simulations the ligands drifted apart from the receptors with no further retention and the forces decays to zero.

thumbnail
Download:
Figure 4. Secondary binding pockets identified in the C1 channel.

Panel a shows the cyanopindolol/β1AR complex and panel b shows the carazolol/β2AR complex. The orientation of these views are the same as in Figures 1a and 1b. Circles show the approximate locations of channels C1 and C2. Ligands are shown in green sticks, and side chains within 3 Å of the ligands are shown in white sticks. Solvent-accessible surfaces of aromatic F359/Y3087.35 and F218/F193 residues are displayed in orange. Panel c depicts the sequence alignment of this region between human βARs. Residues along the extraction trajectories that interact with ligands are highlighted in black, and non-conserved residues are showed in a smaller size.

https://doi.org/10.1371/journal.pone.0023815.g004

thumbnail
Download:
Figure 5. Secondary binding pockets identified in the C2 channel.

Panel a shows the cyanopindolol/β1AR complex and panel b shows the carazolol/β2AR complex. See legend of Figure 4 for details.

https://doi.org/10.1371/journal.pone.0023815.g005

thumbnail
Download:
Table 1. Summary of the residues that interact with the ligand during the SMD pulling experiments through C1 and C2 in β1 and β2ARs.

https://doi.org/10.1371/journal.pone.0023815.t001

Physico-chemical nature and sequence conservation of the entry channels

In both β1- and β2AR, the two identified extracellular channels of ligand entry/exit differ strongly in their physico-chemical properties, as channel C1 is strongly hydrophilic (10 polar/charged residues out of 13) whereas C2 is mainly hydrophobic (7 apolar/aromatic residues out of 10) (summarized in Table 1). Despite this overall similarity between β1AR and β2AR in the fundamental nature of the ligand entry/exit routes, sequence conservation in these regions strongly differs between receptors. Comparison of conserved residues reveals that the sequence identity between β1- and β2AR in channel C1 is only 38% (5 out of 13 residues), while in C2 is 70% (7 out of 10 residues, considering Ile and Val as nearly equivalent).

Characterization of intermediate binding sites

The potential of mean force along extraction coordinates was calculated using the second cumulant expansion of Jarzynski's expression by sampling the work from repeated trajectories [34]. PMF values between starting and ending points were used to estimate free energy changes of dissociation reactions. The free energy cost of moving the ligand from the binding site crevice to bulk water is 7.0 or 6.0 kcal/mol for β1AR, and 5.6 or 6.9 kcal/mol for β2AR, via C1 or C2 channels, respectively (Figures 2 and 3). Clearly, these positive values indicate that receptor-bound states are more favorable in both receptors. Obviously, initial (ligand bound to receptor) and final (ligand in bulk water) states of the SMD simulations, via C1 or C2 channels, are the same, allowing us to estimate the procedure error. The difference in energy of 1.0 and 1.3 kcal/mol, observed for β1AR and β2AR, respectively, between channels C1 and C2, are considered small errors given the complexity of the ligand-receptor-lipid bilayer system. Although no energy minimum was found in the free energy profile, we observed a decrease in the PMF slopes in a narrow region, at distance of ∼9 to 15 Å from the orthosteric binding sites in all experiments (black stars in Figure 2 and 3). These secondary binding pockets correlate with the retention regions identified previously in the C1 and C2 channels and comprise residues located in ECL2 and ECL3, and in the outermost solvent exposed area of TMs (Figures 4 and 5). The free energy cost to move cyanopindolol from the orthosteric binding pocket of β1AR to the secondary binding pocket situated in C1 (2.9 Kcal/mol) is comparable to the value found for the C2 channel (3.2 kcal/mol) and both are located at a distance of ∼9.0 Å from the orthosteric binding site. In contrast, the secondary binding pocket in C2 (5.7 kcal/mol) of β2AR is less favorable than in the C1 channel (3.1 kcal/mol) and is located at ∼15 Å from the orthosteric binding site. In this particular case, additional energy is required to displace the bulky carbazole group of carazolol through the bulky H2.64, I2.65, W3.28 and I7.36 residues in TMs 2, 3 and 7 (Figure 5b).

Discussion

In this work, we have explored the possible exit routes of ligands in the structures of human β1AR and β2AR using SMD simulations. We have found that both receptors have two well-defined access channels from the extracellular side (C1 and C2 in Figure 1). While we explicitly simulate the process of ligand dissociation, the relatively rigid arrangement of the extracellular domains of the receptors strongly suggests that the same channels are also used in the process of ligand entry. During dissociation, ligands are retained in the boundary with the extracellular solvent (∼9–15 Å from the orthosteric binding site, Figures 4 and 5), as evidenced by the decrease in the PMF slopes and larger force values during the SMD experiments (Figure 2 and 3, black stars). We suggest that these retention sites serve as secondary binding pockets during ligand entry. Interestingly, the access channels differ strongly in their physicochemical properties and, particularly, in their degree of sequence conservation (38% identity in C1 vs. 70% identity in C2). However, our simulations produce similar PMF profiles for C1 and C2 in both receptors and, thus, both routes may serve indistinguishably for the entry and exit of inverse agonists. Importantly, all the TM residues identified in our study have been experimentally found to be involved in ligand interactions for βARs or/and other GPCRs: 2.64 [35], [36], 2.65 [37], [38], 3.28 [39], [40], 5.36 [41], 6.55 [42], 6.58 [43], [44], 7.35 [38], [45], 7.36 [46], 7.39 [47] and 7.40 [48]. Also, as the two channels are connected through the orthosteric binding site, we cannot rule out the possibility that ligands could use one route for entry and the other for exit, in the same manner as proposed for the uptake and release of retinal in rhodopsin [23].

Charged residues in ECLs 2 and 3 separate the C1 and C2 channels from each other (Table 1). These residues are D217 and D356 in β1AR and D192 and K305, forming a salt bridge, in β2AR. Importantly, D217 in β1AR and the homologous D192 in β2AR are involved in hydrogen bonding interactions with the protonated group of cyanopindolol and carazolol, respectively, during dissociation via both the C1 and C2 channels (see Figures 4 and 5). We hypothesize that these common negatively charged side chains play an important role to attract the ligand to the channels, and to provide the energy to partially desolvate the ligand. Clearly, extracellular ligands must be transferred from the extracellular aqueous environment to the binding site crevice in the TM domain, away from bulk water. Thus, a crucial contribution to the ligand-receptor binding affinity is the desolvation of the ligand. Interestingly, the corresponding residues in β3AR are non-bulky hydrophobic amino acids, A197 and G325. These remarkable differences are most likely translated into a different pattern of ligand entry in these receptors.

In addition, the C1 and C2 channels are also delineated by F218 in β1AR and F193 in β2AR, located in ECL2, and F3597.35 in β1AR and Y3087.35 in β2AR, located in TM7 (depicted by solvent surfaces in Figures 4 and 5). Previous MD simulations on β2AR have suggested that F193 is able to achieve different conformations [12], [22]. These features were reproduced in our simulations, as we observed a rotation of both the F218 and F193 side chains (black traces in Figures 6a (β1AR, ligand exit through C1), 6b (β1AR, ligand exit through C2), 7a (β2AR, ligand exit through C1) and 7b (β2AR, ligand exit through C2) that parallels the transition of the ligands from the TM bundle into the solvent. However, in contrast with previous works, we observed that the conformational changes of F218 and F193 in ECL2 correlate with an increase in the number of water molecules around ligands during dissociation (grey contour in Figures 6a, 6b, 7a and 7b). Based on this observation, we suggest a novel role for these residues: we hypothesize that in the process of ligand entry F218 and F193 serve as a floodgate by removing the water solvent shell around the compounds during binding.

thumbnail
Download:
Figure 6. Number of water molecules at a distance of 3 Å from cyanopindolol (grey solid contour) and χ1 torsion angle of F218 (black lines) from selected β1AR SMD trajectories through channels C1 (panel a) and C2 (panel b).

https://doi.org/10.1371/journal.pone.0023815.g006

thumbnail
Download:
Figure 7. Number of water molecules at a distance of 3 Å from carazolol (grey solid contour) and χ1 torsion angle of F193 (black lines) from selected β2AR SMD trajectories through channels C1 (panel a) and C2 (panel b).

https://doi.org/10.1371/journal.pone.0023815.g007

The extraordinary variability in length and amino acid composition of the extracellular loops across the GPCR superfamily generates a wide recognition space for ligands with very diverse chemical scaffolds. For instance ECL 2 of rhodopsin, formed by two β-strands, buries the binding site from the extracellular environment, whereas ECL 2 of CXCR4, also formed by two β-strands, fully exposes the binding site to the extracellular environment. In contrast, a helical segment forms ECL 2 of the β1- and β2- adrenergic receptors. This α-helix is probably not conserved in the other members of the biogenic amine receptor family, as it was not found in the structure of the dopamine D3 receptor. It was recently shown that ECLs 2 and 3 of the β2-adrenergic receptor exist in three distinct conformations depending on the type of ligand bound to the TM core (neutral antagonists, agonists, or inverse agonists) [12], [22]. Thus, this extracellular domain of the receptor plays a key role in receptor activation. We hypothesize that small molecules binding to these secondary-binding pockets, in the extracellular domain, might act as allosteric modulators.

Methods

Molecular models and identification of ligand access channels

The high-resolution crystal structures of the β1AR [19] and β2AR [20] were obtained from the Protein Data Bank (PDB entries 2VT4 and 2RH1 respectively). MODELER [49] was used to transform the starting coordinates of the turkey β1AR (UniProtKB: P07700) to the human sequence (UniProtKB: P08588). It is important to note that major differences between turkey and human sequences are present in the N- and C-terminal regions (e.g. human β1AR have an N-terminal domain 17 residues longer). The notation of the β1AR amino acids in the manuscript corresponds to the human sequence. CAVER [29] was used to determine channels connecting the ligand binding site to the extracellular surface in snapshot structures (every 0.5 ns) along the equilibration period (see below). The initial state for cavities search was at the center of mass of the ligands and a grid spacing of 0.5 Å was used. This approach leads to the identification of two channels in both receptors (C1 and C2 in Figure 1). In addition, we include two inter-helical channels (C3 and C4) calculated by the same procedure for the GPCR opsin [23] and superimposed onto the βARs coordinates. These “rhodopsin-like” channels, however, were not detected by CAVER in the βARs structures.

Molecular dynamics (MD) simulations

The β1AR and β2AR human receptors in ligand bound conformation and nine internal water molecules in the P6.50/D2.50/N7.49/Y7.53 environment [50] were embedded in a pre-equilibrated lipid bilayer consisting of 282 molecules of 1-palmitoyl-2-oleoyl-sn-glycerol-3-phosphatidylcholine (POPC). These crystallographic water molecules did not displace significantly from their starting positions during the simulations (data not shown). Electroneutrality of the system was achieved by adding chloride ions to fulfill a net charge of zero; then, additional sodium and chloride ions were added to a final concentration of 0.1 mol/L. Simulations were carried out using the NAMD version 2.7 MD package [51] using the TIP3 water model and the CHARMM27 all-hydrogen force field [52]. Atomic charges for carazolol and cyanopindolol were calculated with HF/6-31G* and RESP [53], and compared against the corresponding atom types in the CGenFF [54]. In all cases, we only observed small differences in values, while the signs of the charges were always maintained. Long-range electrostatic interactions were calculated using the particle mesh Ewald (PME) method [55]. Initial coordinates were optimized by energy minimization. After geometry optimization, the temperature of the systems was raised in 30.000 steps by temperature reassignment method followed by 10 ns of equilibration at 300 K and constant pressure.

Steered molecular dynamics (SMD) simulations

The SMD method implemented in NAMD [24] was used to simulate ligands dissociation. The directions of the applied forces (reaction coordinate) were vectors with origin in the center of mass of the ligands and having minimal standard deviation from the path graph nodes defined by CAVER. SMD simulations were performed at constant velocity of 10 Å/ns and the spring constant was set to 250 pN/Å. These parameters were similar to those used previously in biological systems and sufficient to ensure that the work distribution is Gaussian [56]. Each trajectory was carried out until the ligands were displaced towards the receptor surface, and was repeated 6 times. The pulling force F at time t was calculated according to:(1)where k is the spring constant, v is the constant velocity of pulling, r0 and r(t) are the ligand center of mass position at initial and current time t respectively, is the direction of the pulling vector. The potential of mean force (PMF) along the reaction coordinate was calculated by the second-order cumulant expansion of the irreversible work measurements [34] according to:(2)(3)where 〈W〉 is the mean work averaged from the six trajectories, kB is Boltzmann's constant and T is the bulk temperature.

The Jarzynski's equality applied in this study is relative easy to implement compared to other free energy methods such as umbrella sampling. However, it is not exempt of the inaccuracies inherent to the insufficient sampling of the configuration space [57]. Thus, we have only used the PMF profiles as a guideline for the identification of residues involved in interaction with the ligands during the extraction process. Specifically, we do not to aim to compare the theoretical energy values with experimental binding affinities.

Supporting Information

Figure S1.

Rmsd values of the backbone atoms of β1AR (a) and β2AR (b) along the trajectories of the MD equilibrium simulations of the receptor-membrane systems.

https://doi.org/10.1371/journal.pone.0023815.s001

(TIF)

Figure S2.

Representative force profiles of ligand extraction along the C1–C4 channels. Panel a corresponds to the cyanopindolol/β1AR complex and panel b corresponds to the carazolol/β2AR complex. C1 and C2 correspond to extracellular routes whereas C3 and C4 correspond to routes that lead to the membrane core.

https://doi.org/10.1371/journal.pone.0023815.s002

(TIF)

Author Contributions

Conceived and designed the experiments: XD LP. Performed the experiments: AG. Analyzed the data: AG XD. Contributed reagents/materials/analysis tools: TPA. Wrote the paper: AG XD LP.

References

  1. 1. Fredriksson R, Lagerstrom MC, Lundin LG, Schioth HB (2003) The g-protein-coupled receptors in the human genome form five main families. Phylogenetic analysis, paralogon groups, and fingerprints. Mol Pharmacol 63: 1256–1272.
  2. 2. Fagerberg L, Jonasson K, von Heijne G, Uhlen M, Berglund L (2010) Prediction of the human membrane proteome. Proteomics 10: 1141–1149.
  3. 3. Smit MJ, Vischer HF, Bakker RA, Jongejan A, Timmerman H, et al. (2007) Pharmacogenomic and structural analysis of constitutive g protein-coupled receptor activity. Annu Rev Pharmacol Toxicol 47: 53–87.
  4. 4. Imming P, Sinning C, Meyer A (2006) Drugs, their targets and the nature and number of drug targets. Nat Rev Drug Discov 5: 821–834.
  5. 5. Palczewski K, Kumasaka T, Hori T, Behnke CA, Motoshima H, et al. (2000) Crystal structure of rhodopsin: A G protein-coupled receptor. Science 289: 739–745.
  6. 6. Lagerstrom MC, Schioth HB (2008) Structural diversity of G protein-coupled receptors and significance for drug discovery. Nat Rev Drug Discov 7: 339–357.
  7. 7. Gkountelias K, Tselios T, Venihaki M, Deraos G, Lazaridis I, et al. (2009) Alanine scanning mutagenesis of the second extracellular loop of type 1 corticotropin-releasing factor receptor revealed residues critical for peptide binding. Mol Pharmacol 75: 793–800.
  8. 8. Samson M, LaRosa G, Libert F, Paindavoine P, Detheux M, et al. (1997) The second extracellular loop of CCR5 is the major determinant of ligand specificity. J Biol Chem 272: 24934–24941.
  9. 9. Avlani VA, Gregory KJ, Morton CJ, Parker MW, Sexton PM, et al. (2007) Critical role for the second extracellular loop in the binding of both orthosteric and allosteric G protein-coupled receptor ligands. J Biol Chem 282: 25677–25686.
  10. 10. Klco JM, Wiegand CB, Narzinski K, Baranski TJ (2005) Essential role for the second extracellular loop in C5a receptor activation. Nat Struct Mol Biol 12: 320–326.
  11. 11. Scarselli M, Li B, Kim SK, Wess J (2007) Multiple residues in the second extracellular loop are critical for M3 muscarinic acetylcholine receptor activation. J Biol Chem 282: 7385–7396.
  12. 12. Bokoch MP, Zou Y, Rasmussen SGF, Liu CW, Nygaard R, et al. (2010) Ligand-specific regulation of the extracellular surface of a G-protein-coupled receptor. Nature 463: 108–112.
  13. 13. Tate CG, Schertler GFX (2009) Engineering G protein-coupled receptors to facilitate their structure determination. Curr Opin Struct Biol 19: 386–395.
  14. 14. Rosenbaum DM, Rasmussen SG, Kobilka BK (2009) The structure and function of G-protein-coupled receptors. Nature 459: 356–363.
  15. 15. Deupi X, Standfuss J (2011) Structural insights into agonist-induced activation of G-protein-coupled receptors. Curr Opin Struct Biol.
  16. 16. Wu B, Chien EY, Mol CD, Fenalti G, Liu W, et al. (2010) Structures of the CXCR4 chemokine GPCR with small-molecule and cyclic peptide antagonists. Science 330: 1066–1071.
  17. 17. Chien EY, Liu W, Zhao Q, Katritch V, Han GW, et al. (2010) Structure of the human dopamine D3 receptor in complex with a D2/D3 selective antagonist. Science 330: 1091–1095.
  18. 18. Shimamura T, Shiroishi M, Weyand S, Tsujimoto H, Winter G, et al. (2011) Structure of the human histamine H(1) receptor complex with doxepin. Nature.
  19. 19. Warne T, Serrano-Vega MJ, Baker JG, Moukhametzianov R, Edwards PC, et al. (2008) Structure of a beta1-adrenergic G-protein-coupled receptor. Nature 454: 486–491.
  20. 20. Rasmussen SG, Choi HJ, Rosenbaum DM, Kobilka TS, Thian FS, et al. (2007) Crystal structure of the human beta2 adrenergic G-protein-coupled receptor. Nature 450: 383–387.
  21. 21. Baker JG (2005) The selectivity of beta-adrenoceptor antagonists at the human beta1, beta2 and beta3 adrenoceptors. Br J Pharmacol 144: 317–322.
  22. 22. Wang T, Duan Y (2009) Ligand entry and exit pathways in the beta2-adrenergic receptor. Journal of Molecular Biology 392: 1102–1115.
  23. 23. Hildebrand PW, Scheerer P, Park JH, Choe H-W, Piechnick R, et al. (2009) A ligand channel through the G protein coupled receptor opsin. PLoS ONE 4: e4382.
  24. 24. Isralewitz B, Baudry J, Gullingsrud J, Kosztin D, Schulten K (2001) Steered molecular dynamics investigations of protein function. J Mol Graph Model 19: 13–25.
  25. 25. Isralewitz B, Izrailev S, Schulten K (1997) Binding pathway of retinal to bacterio-opsin: a prediction by molecular dynamics simulations. Biophys J 73: 2972–2979.
  26. 26. Jensen MØ, Park S, Tajkhorshid E, Schulten K (2002) Energetics of glycerol conduction through aquaglyceroporin GlpF. Proc Natl Acad Sci USA 99: 6731–6736.
  27. 27. Fishelovitch D, Shaik S, Wolfson HJ, Nussinov R (2009) Theoretical characterization of substrate access/exit channels in the human cytochrome P450 3A4 enzyme: involvement of phenylalanine residues in the gating mechanism. J Phys Chem B 113: 13018–13025.
  28. 28. Yang LJ, Zou J, Xie HZ, Li LL, Wei YQ, et al. (2009) Steered molecular dynamics simulations reveal the likelier dissociation pathway of imatinib from its targeting kinases c-Kit and Abl. PLoS ONE 4: e8470.
  29. 29. Petrek M, Otyepka M, Banas P, Kosinova P, Koca J, et al. (2006) CAVER: a new tool to explore routes from protein clefts, pockets and cavities. BMC Bioinformatics 7: 316.
  30. 30. Cherezov V, Rosenbaum DM, Hanson MA, Rasmussen SG, Thian FS, et al. (2007) High-resolution crystal structure of an engineered human beta2-adrenergic G protein-coupled receptor. Science 318: 1258–1265.
  31. 31. Park JH, Scheerer P, Hofmann KP, Choe HW, Ernst OP (2008) Crystal structure of the ligand-free G-protein-coupled receptor opsin. Nature 454: 183–187.
  32. 32. Liu X, Xu Y, Wang X, Barrantes FJ, Jiang H (2008) Unbinding of nicotine from the acetylcholine binding protein: steered molecular dynamics simulations. J Phys Chem B 112: 4087–4093.
  33. 33. Ballesteros JA, Weinstein H (1995) Integrated methods for the construction of three dimensional models and computational probing of structure-function relations in G-protein coupled receptors. Methods in Neurosciences 25: 366–428.
  34. 34. Park S, Khalili-Araghi F, Tajkhorshid E, Schulten K (2003) Free energy calculation from steered molecular dynamics simulations using Jarzynski's equality. J Chem Phys 119: 3559.
  35. 35. Xie KQ, Cao Y, Zhu XZ (2006) Role of the second transmembrane domain of rat adenosine A1 receptor in ligand-receptor interaction. Biochem Pharmacol 71: 865–871.
  36. 36. Simpson MM, Ballesteros JA, Chiappa V, Chen J, Suehiro M, et al. (1999) Dopamine D4/D2 receptor selectivity is determined by A divergent aromatic microdomain contained within the second, third, and seventh membrane-spanning segments. Mol Pharmacol 56: 1116–1126.
  37. 37. Hoffmann SH, ter Laak T, Kuhne R, Reilander H, Beckers T (2000) Residues within transmembrane helices 2 and 5 of the human gonadotropin-releasing hormone receptor contribute to agonist and antagonist binding. Mol Endocrinol 14: 1099–1115.
  38. 38. Sugimoto Y, Fujisawa R, Tanimura R, Lattion AL, Cotecchia S, et al. (2002) Beta(1)-selective agonist (-)-1-(3,4-dimethoxyphenetylamino)-3-(3,4-dihydroxy)-2-propanol [(-)-RO363] differentially interacts with key amino acids responsible for beta(1)-selective binding in resting and active states. J Pharmacol Exp Ther 301: 51–58.
  39. 39. Hogan K, Peluso S, Gould S, Parsons I, Ryan D, et al. (2006) Mapping the binding site of melanocortin 4 receptor agonists: a hydrophobic pocket formed by I3.28(125), I3.32(129), and I7.42(291) is critical for receptor activation. J Med Chem 49: 911–922.
  40. 40. Gerber BO, Meng EC, Dotsch V, Baranski TJ, Bourne HR (2001) An activation switch in the ligand binding pocket of the C5a receptor. J Biol Chem 276: 3394–3400.
  41. 41. Gether U, Nilsson L, Lowe JA 3rd, Schwartz TW (1994) Specific residues at the top of transmembrane segment V and VI of the neurokinin-1 receptor involved in binding of the nonpeptide antagonist CP 96,345 [corrected]. J Biol Chem 269: 23959–23964.
  42. 42. Wieland K, Zuurmond HM, Krasel C, Ijzerman AP, Lohse MJ (1996) Involvement of Asn-293 in stereospecific agonist recognition and in activation of the beta 2-adrenergic receptor. Proc Natl Acad Sci U S A 93: 9276–9281.
  43. 43. Spalding TA, Burstein ES, Henderson SC, Ducote KR, Brann MR (1998) Identification of a ligand-dependent switch within a muscarinic receptor. J Biol Chem 273: 21563–21568.
  44. 44. Hovelmann S, Hoffmann SH, Kuhne R, ter Laak T, Reilander H, et al. (2002) Impact of aromatic residues within transmembrane helix 6 of the human gonadotropin-releasing hormone receptor upon agonist and antagonist binding. Biochemistry 41: 1129–1136.
  45. 45. Bonner G, Meng F, Akil H (2000) Selectivity of mu-opioid receptor determined by interfacial residues near third extracellular loop. Eur J Pharmacol 403: 37–44.
  46. 46. Jarnagin K, Bhakta S, Zuppan P, Yee C, Ho T, et al. (1996) Mutations in the B2 bradykinin receptor reveal a different pattern of contacts for peptidic agonists and peptidic antagonists. J Biol Chem 271: 28277–28286.
  47. 47. Suryanarayana S, Kobilka BK (1993) Amino acid substitutions at position 312 in the seventh hydrophobic segment of the beta 2-adrenergic receptor modify ligand-binding specificity. Mol Pharmacol 44: 111–114.
  48. 48. Roth BL, Shoham M, Choudhary MS, Khan N (1997) Identification of conserved aromatic residues essential for agonist binding and second messenger production at 5-hydroxytryptamine2A receptors. Mol Pharmacol 52: 259–266.
  49. 49. Sali A, Blundell TL (1993) Comparative protein modelling by satisfaction of spatial restraints. J Mol Biol 234: 779–815.
  50. 50. Pardo L, Deupi X, Dolker N, Lopez-Rodriguez ML, Campillo M (2007) The role of internal water molecules in the structure and function of the rhodopsin family of G protein-coupled receptors. Chembiochem 8: 19–24.
  51. 51. Phillips JC, Braun R, Wang W, Gumbart J, Tajkhorshid E, et al. (2005) Scalable molecular dynamics with NAMD. J Comput Chem 26: 1781–1802.
  52. 52. MacKerell JAD, Bashford D, Bellott M, Dunbrack RL Jr, Evanseck JD, et al. (1998) All-atom empirical potential for molecular modeling and dynamics studies of proteins. Journal of Physical Chemistry B 102: 3586–3616.
  53. 53. Bayly CI, Cieplak P, Cornell WD, Kollman PA (1993) A well-behaved electrostatic potential based method using charge restraints for deriving atomic charges. J Chem Phys 97: 10269–10280.
  54. 54. Vanommeslaeghe K, Hatcher E, Acharya C, Kundu S, Zhong S, et al. (2010) CHARMM general force field: A force field for drug-like molecules compatible with the CHARMM all-atom additive biological force fields. JOURNAL OF COMPUTATIONAL CHEMISTRY 31: 671–690.
  55. 55. Patra M, Karttunen M, Hyvonen MT, Falck E, Lindqvist P, et al. (2003) Molecular dynamics simulations of lipid bilayers: major artifacts due to truncating electrostatic interactions. Biophys J 84: 3636–3645.
  56. 56. Minh DDL, McCammon JA (2008) Springs and speeds in free energy reconstruction from irreversible single-molecule pulling experiments. J Phys Chem B 112: 5892–5897.
  57. 57. Gore J, Ritort F, Bustamante C (2003) Bias and error in estimates of equilibrium free-energy differences from nonequilibrium measurements. Proc Natl Acad Sci U S A 100: 12564–12569.