Membrane Partitioning Characteristics of Salmeterol, Formoterol, and Salbutamol

To elucidate the differences in membrane partitioning characteristics, we performed steered MD and umbrella sampling simulations using a model membrane containing 1-palmitoyl-2-oleoyl-sn-glycero-3-phosphocholine (POPC) and cholesterol. The obtained results presented below provide remarkable insights into their distinct membrane interactions in atomistic details.

Salmeterol Has a More Favorable Free Energy of Solvation in the POPC Lipid Bilayer

The PMF curves obtained from the simulations provide the free energy profiles of solvation of each ligand from water to the membrane and along the bilayer normal (z-axis). The bilayer regions at which the solutes exhibit the most and least favorable energies are characterized by energy minima and maxima, respectively (Fig. 1B). Each of the three ligands shows a favorable (negative) free energy of transfer from water to the membrane. The free energy minima of the COM of salbutamol, salmeterol, and formoterol are located at |z_min| ≈ 18, |z_min| ≈ 15, and |z_min| ≈ 14, respectively. The free energy of partitioning for salmeterol was the most favorable (ΔG_partitioning = −3.15 ± 0.02 kcal/mol), followed by formoterol (ΔG_partitioning = −2.53 ± 0.02 kcal/mol) and salbutamol (ΔG_partitioning = −0.98 ± 0.03 kcal/mol). As the COM of each molecule moved across the membrane toward the bilayer center (z = 0), the free energy barrier for crossing the membrane (ΔG_crossing) increased in the reverse order to that of ΔG_partitioning. Salbutamol showed the largest energy barrier (ΔG_crossing = 13.21 ± 0.11 kcal/mol), followed by formoterol (ΔG_crossing = 11.09 ± 0.07 kcal/mol) and salmeterol (ΔG_crossing = 10.39 ± 0.07 kcal/mol). In all cases, the thermally accessible regions (energy barrier of RT = 0.616 kcal/mol at 310 K) extended for 2 Å on both sides around the energy minima (Table 1). It is important to recognize that, unlike salbutamol, both salmeterol and formoterol are lengthy molecules, and their long axes span ∼20 Å and ∼15 Å, respectively, and the time-average bilayer locations of their COMs, i.e., |z_min| values reported above, need to be examined carefully. To obtain precise and detailed bilayer locations and orientations of the molecules, we calculated the mass density profiles for the ligands and various components of the membrane lipids along the bilayer normal (Supplemental Fig. 1). The mass density for salmeterol clearly indicates that the molecule is present in its extended conformation, the long axis of which spans almost half of the entire bilayer. Formoterol’s density profile is similar to that of salmeterol, and it is less spread out due to its smaller size. Salbutamol has the least spread, indicating its preferred location around the lipid headgroup region and the absence of any notable interactions with the lipid core. Whereas the saligenin-ethanolamine head of each ligand is located at the lipid headgroup region (Fig. 1C), the lipophilic tails of salmeterol and formoterol are oriented perpendicular to the membrane plane (i.e., in parallel to the bilayer normal) and buried into the lipid core, clearly displaying their amphiphilic nature.

The Bilayer Orientation and H-Bond Interactions of Saligenin-Ethanolamine Polar Groups Are Distinct among the Ligands

We performed comprehensive analyses of the orientation of the saligenin-ethanolamine head of each ligand and their interactions with various functional groups of the lipids to obtain more insight into the observed differences in the membrane solvation free energies among the ligands. The tilt angle measured between the longitudinal vector of the saligenin-ethanolamine head and the bilayer normal represents the orientational angle and is presented as probability density graphs (Fig. 1D, see Materials and Methods section for details). The H-bond interactions between various polar groups (three -OH groups labeled as O1, O2, and O3, NH₃⁺ group labeled as N1, and -NH group present only in formoterol, which is labeled as N2; Fig. 1A) of the ligands and the lipid headgroups including choline, phosphate, and glyceryl carbonyls were monitored through the simulation time and presented as percent H-bond occupancy (Fig. 1E and Supplemental Fig. 2, see Materials and Methods section for details). In addition, we calculated the atom contacts (as percent occupancy) that give the fraction of total simulation time during which the functional groups from the ligand are located within 4 Å of various lipid functional groups (Supplemental Fig. 3).

Intriguingly, the orientation angles are unique for each ligand, resulting from subtle differences in interaction patterns among their common polar groups as well as specific interactions of additional functional groups that are present only in salmeterol or formoterol. Salmeterol’s saligenin-ethanolamine head orients in two distinct modes, as shown by the bimodal distribution of its orientation angle (Fig. 1D). The major orientation is one in which the head is at an acute angle to the bilayer normal (θ = 40° ± 20°). In the secondary mode, the head orients perpendicular to the bilayer normal (θ = 90° ± 20°). The bimodal distribution of the orientation angle for salmeterol, reflecting its flexible nature, can be attributed to the number of rotatable bonds, which is twice as high as that of formoterol (16 vs. 8; Supplemental Table 1). In the case of formoterol, the saligenin-ethanolamine head is oriented at an acute angle (θ = 60° ± 30°) that is approximately perpendicular to the bilayer normal. Salbutamol’s orientation is vastly different from salmeterol and formoterol and can be seen mostly at an obtuse angle to the bilayer normal (θ = 140° ± 20°) in which its phenolic hydroxyl group is oriented toward the lipid alkyl tails.

These unique orientations and conformations are due to the specific interactions between the functional groups located on each ligand and the membrane lipids. The number of H-bonds made by salmeterol (at its low-energy well) is relatively less than the other two ligands (Supplemental Fig. 2, A and B). Whereas the aliphatic ethanolamine hydroxyl group (O3) makes H-bond interactions with oxygen atoms of the lipid phosphate and glyceryl groups in all the ligands, the interactions of both the phenolic OH (O1) and benzyl OH (O2) groups are different. Salmeterol’s O1 is located near the choline group, and its H-bond occupancy with the phosphate oxygen is much less compared with that of formoterol. Formoterol’s O1 makes H-bond interactions with both phosphate and glyceryl oxygens, although the occupancy for phosphate oxygens is higher compared with the glyceryl ones. On the other hand, salbutamol’s O1 is located near the alkyl tail and makes H-bond interaction mainly with glyceryl oxygen atoms (Fig. 1E). The benzyl hydroxy group (O2) in salmeterol and salbutamol follow similar trends to their respective O1 groups but with much higher occupancies. However, the formamide carbonyl oxygen (labeled as O2 in Fig. 1A) of formoterol does not engage in H-bond interactions with the lipid headgroups but rather is oriented toward the aqueous bulk and forms stable H-bonds with the water molecules (Fig. 1C).

The charged ethanolamine amino group (labeled as N1), in all three ligands, is located near the lipid phosphate groups and makes H-bond/electrostatic interactions with the phosphate oxygens. This amino group also makes such interactions with glyceryl oxygens (only in salmeterol and formoterol, but not in salbutamol). The formamide -NH group of formoterol (labeled as N2) is primarily located near the glyceryl/phosphate groups and engages in H-bond interactions with them. The hydrophobic tails of salmeterol and formoterol (including the oxygen atom, labeled as O4 in Fig. 1A) are mostly buried and sandwiched between the alkyl tails of the lipid core region.

The time-average bilayer locations of the ligands at the lipid headgroup/core interface and their distinct orientations clearly indicate that they remain in contact with the aqueous phase while they are partitioned into the membrane. To quantify these critical solvation characteristics, we calculated the number of water molecules within 4 Å of any heavy atom of the ligands and normalized the numbers by their respective heavy atom counts to account for the notable variation in their sizes (Table 1). The numbers are significantly different among the ligands (Supplemental Fig. 2C). Interestingly, salmeterol has significantly fewer water contacts compared with formoterol and salbutamol because only 13 (out of 30) of its heavy atoms are located at the lipid headgroup region, and the majority of them (>56%) are buried deep in the lipid core. In the case of formoterol, 14 (out of 25) of its heavy atoms are in contact with the lipid headgroup, and only 11 (44%) of its heavy atoms are in contact with the lipid core. In contrast, none of the salbutamol atoms are buried in the lipid core region. These unique solvation characteristics of the ligands appear to dictate their association paths (either by aqueous route or by diffusion through the membrane) through which the ligand access and reach the binding site of β2-AR.

Binding Pathways for Salmeterol, Formoterol, and Salbutamol

We examined the plausible receptor access and binding pathways of the three ligands, taking into consideration their membrane partitioning characteristics. All association simulations started with the ligands at their energetically favorable (time-average) locations and orientations within the membrane. For salbutamol, several simulations were run with the ligand in the aqueous phase, above and around the receptor as its starting position. Earlier studies investigating the binding of ligands to orthosteric sites of GPCRs demonstrated that such events occur spontaneously only in microsecond timescales (Dror et al., 2011, 2013). Therefore, to accelerate the ligand access and binding processes, we used WT-MetaD, an enhanced sampling technique that utilizes a biasing force and the specifics of the binding site residues to guide the ligand-receptor association (see Materials and Methods section for details).

Salmeterol Enters the Binding Pocket through TMH1 and 7 from within the Bilayer

Out of 20 WT-MetaD association simulations aimed to characterize the binding mechanism of salmeterol, we observed a complete association of the ligand in 18 simulations. The ligand association process was considered complete only when the distance between the COMs of the ligand and that of the binding site residues is within 4 Å. In 15 of 18 successful simulations, salmeterol reached the binding pocket by entry through its saligenin headfirst. Intriguingly, in the rest of the three simulations, salmeterol entered the pocket with its phenyl tail first. However, in all the cases, salmeterol approached the receptor from within the membrane through TMHs 1, 2, and 7. All the association simulations started with salmeterol partitioned within the membrane in its energetically favorable bilayer depth and orientation, at a minimum distance of ∼20 Å from the receptor binding site. Specifically, salmeterol was placed near TMHs 1 and 7 as several unbiased simulations indicated that salmeterol preferentially interacts with residues from these helices. However, a spontaneous entry of salmeterol into the binding site was not observed in any of these simulations.

The free energy surface, describing the thermodynamics of the ligand access and receptor-binding events, including multiple low-energy intermediate states (A-D), was characterized by two collective variables: the distance between the COM of the ligand (heavy atoms) and the binding site residues (x-axis) and the internal angle of the ligand accounting for the conformational changes (y-axis) (Fig. 2). As salmeterol approached the receptor, its hydrophobic aralkyl-oxy-alkyl tail made contact with several residues from TMHs 1 and 7, including W313^7.40, I309^7.36, L310^7.37, V39^1.38, W42^1.41, and I43^1.42, whereas the saligenin head appeared to explore nearby residues for potential interactions (Fig. 2A). In most of the simulations, E306^7.33 made the first contact with the upper end of salmeterol through H-bond interactions with its O1 and O2 groups. These interactions drew the ligand closer to the binding pocket and anchored it in a position favorable for entry.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 2.

The free energy surface for salmeterol’s access and binding to the β2-adrenergic receptor through the TMHs 1 and 7 from its energetically favorable location in the membrane. The two-dimensional free energy surface is characterized by the distance between the COM of the ligand and binding site residues (x-axis) in nm and the internal angle of the salmeterol molecule (y-axis) in degrees. The minimum energy path, as determined by well tempered metadynamics simulations, is given in bold black line, and several representative intermediate states along the path were labeled (A–D). (A) As salmeterol (licorice representation, magenta color) approaches the receptor from within the membrane, its aryl-alkoxy-alkyl tail comes into contact with the hydrophobic residues, including W313^7.40, I309^7.36, L310^7.37, V39^1.38, W42^1.41, and I43^1.42 from TMHs 1 and 7. (B and C) Salmeterol undergoes a significant conformational change from its extended to bent form and engages in H-bond interactions with D192^ECL2 through its O1 and O2 groups, which destabilize the salt bridge between D192^ECL2 and K305^7.32. The breakage of the salt bridge allows salmeterol to move further, disrupting the hydrophobic lock (between F193^ECL2 and Y308^7.35) while the tail is still anchored by the hydrophobic residues. (D) Outward movement of the side chain indole ring of W313^7.40 releases the tail and facilitates salmeterol’s entry into the pocket, which then assumes its final bound pose, mostly similar to that of the crystal structure pose (PDB ID 6MXT) (Masureel et al., 2018). In its final bound pose, salmeterol engages in polar H-bond and salt bridge interactions with S203^5.42, S207^5.46, D113^3.32, and N312^7.39. Salmeterol (crystal pose in dark blue) and the binding site residues (green) are illustrated in licorice representation. The receptor is illustrated in secondary structure representation (light blue).

As salmeterol prepared to enter the binding pocket through the helices, it underwent a significant conformational change from its extended form to a bent one, with the change in the internal angle of the ligand from ∼150° to ∼50° (from A to B) as can be seen in the free energy surface diagram (Fig. 2). The rearrangement caused the saligenin head to proceed toward the site forming H-bonds with D192^ECL2 through O1 and O2 groups (Fig. 2B). As the saligenin head of salmeterol moved further toward the pocket, its N1 group formed an H-bond with E306^7.33, whereas O1 and O2 groups engaged with the main chain -NH group of F193^ECL2, D192^ECL2, and K305^7.33, causing the breakage of the D192^ECL2-K305^7.33 salt bridge if the bond was already formed. These interactions drew salmeterol’s head further into the binding site, whereas the tail remained anchored by the residues from TMHs 1 and 7, including V39^1.38, I43^1.42, I309^7.36, L310^7.37, and W313^7.40 (Fig. 2C). As the saligenin head entered the binding pocket, the hydrophobic lock between F193^ECL2 and Y308^7.35 was also broken, which reformed briefly after breaking earlier in the simulation. Eventually, salmeterol slid into the binding pocket and advanced toward its final bound pose with its saligenin head O1 and O2 groups engaging in H-bond interactions to the side chain -OH groups of S207^5.46 and S203^5.42, respectively (Fig. 2D). The phenyl end was able to reach the pocket only when the side chain indole ring of W313^7.40 rotated outward and pushed the phenyl ring away from the anchoring point while closing the transmembrane entry. In its final bound pose, the N1 and O3 groups of the ligand formed polar interactions with D113^3.32 and N312^7.39, respectively. Interestingly, the phenyl tail adopted a slightly different position relative to the crystal structure pose. The phenyl tail drifts to the other side of the ECL2 region, close to R175^ECL2 and Y174^ECL2. However, in the crystal structure pose, there is an H-bond between the ether oxygen of salmeterol and the main chain of F193^ECL2.

Once the saligenin head reached the pocket and participated in interactions similar to the crystal structure bound pose, we switched off the bias and continued the simulations unbiased for another 200 ns to observe whether salmeterol’s tail would assume a similar conformation as that of the crystal pose. Although not identical, we observed a binding pose in which the conformation of the tail was similar to that of the crystal pose, with an overall ligand RMSD of 2.56 Å, as shown in Fig. 2D. Another notable binding pose of salmeterol observed during this extended unbiased simulation was one in which the ligand tail adopted an extended conformation with its phenyl ring positioned around K97^2.63, K305^7.32, and H93^2.64.

Salmeterol’s Lipophilic Tail Acts as a Passkey for the Receptor Entry

Although in the majority of the simulations, salmeterol made entry into the binding pocket with its saligenin head first, in three out of the 18 successful simulations, salmeterol entered in a strikingly different manner (Fig. 3). Initially, in its bent conformation as described above, salmeterol engaged with D192^ECL2 and K305^7.33 through the saligenin -OH groups O1 and O2, whereas its aryl-alkyl tail was anchored by the residues from TMHs 1 and 7 (W313^7.40, V317^7.44, I43^1.42, and V39^1.38). As the salt bridge interactions were broken (between D192^ECL2 and K305^7.32), the saligenin head proceeded to pass between the residues, which was prevented by a newly formed bond between the side chain amino group of K305^7.33 and the main chain oxygen of F193^ECL2. Salmeterol’s attempt to enter the site with its saligenin head was also blocked by the phenyl ring side chains of F193^ECL2 and F194^ECL2, which appeared to act as a gate (Fig. 3B). In a fascinating manner, salmeterol rearranged, withdrawing its saligenin head from the gate and presenting its aryl-alkyl tail as a passkey instead (Fig. 3C and Supplemental Video 1).

• Download figure
• Open in new tab
• Download powerpoint

Fig. 3.

Salmeterol enters the β2-AR binding site from within the membrane through the transmembrane helices 1 and 7 using its lipophilic tail as a passkey. (A) Hydrophobic residues from TMHs 1 and 7 anchor the lipophilic tail of salmeterol and facilitate salmeterol’s entry into the pocket from within the membrane. (B) However, in several simulations, salmeterol’s entry by its saligenin head is blocked by the side chain phenyl ring of F193^ECL2 that acts as a gate. (C) Fascinatingly, salmeterol flips 180° and presents its aryl-alkyl tail as a passkey to the gate, which immediately opens up the channel by flipping F193^ECL2 upward to gain entry into the pocket.

Surprisingly, the gate opened as the phenyl ring of F193^ECL2 flipped upward, allowing salmeterol to enter the pocket. It is worth noting that F193^ECL2 has been reported to be important for the higher selectivity of salmeterol for the β₂-AR over β₁-AR relative to other commonly used β₂-AR agonists (Carter and Hill, 2005; Baker, 2010). The tail appears to push apart TMHs 7 and 2, breaking the newly formed interactions between K305^7.33 and F193^ECL2 as well as pushing the hydrophobic lock (Y308^7.35-F193^ECL2) further apart to facilitate the entry. As the molecule stretched across the binding pocket, the tail end was seen exploring regions between TMH4 and TMH5, which allowed the saligenin head to slide into the pocket to assume its final bound pose identical to other simulations.

Formoterol Preferentially Interacts with Aqueous Bulk before Entering the Binding Pocket

In the case of formoterol, a clearly predominant association path was observed in the majority of the simulations. As the ligand approaches the receptor from within the membrane, it undergoes a slight conformational change such that the methoxyphenyl alkyl tail makes notable hydrophobic contacts with residues I94^2.65, L311^7.38, and W313^7.40 of TMHs 2 and 7, respectively (Fig. 4A). After these transient interactions, the ligand moves up to a region just above the transmembrane helices, formed by residues from TMHs 2 and 6 as well as ECL2, and then enters the binding pocket in a manner that is similar to other polar ligands. Intriguingly, the entry of formoterol from this region followed two distinct pathways, remarkably governing the pace at which the ligand reaches the binding pocket.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 4.

The free energy surface for formoterol’s access and binding to the β2-adrenergic receptor from its energetically favorable location in the membrane. The two-dimensional free energy surface is characterized by the distance between the COM of the ligand and binding site residues (x-axis) and the internal angle of formoterol (y-axis). The minimum energy path, as determined by well tempered metadynamics simulations, is given in bold black line, and several representative low-energy intermediate states along the path were labeled (A–D). (A) As formoterol (licorice representation, yellow color) approaches the receptor from within the membrane, its alkoxy-aryl-alkyl tail comes into contact with TMHs 2 and 7, engaging in hydrophobic contacts with I94^2.65, L310^7.37, and W313^7.40. (B and C) Formoterol moves up over the transmembrane helices and remains in a region where it makes extensive H-bond interactions with several polar residues and bulk water molecules. Specifically, the H-bond with D192^ECL2 through its O3 appears to destabilize the salt bridge between D192^ECL2 and K305^7.32. The breakage of the salt bridge allows formoterol to slide down from this region into the binding pocket. As formoterol enters the pocket, the hydrophobic lock (between F193^ECL2 and Y308^7.35) breaks, resulting in the formation of a π-π interaction between F193’s phenyl ring and the aromatic ring at the tail end of the ligand. (D) Formoterol is seen in its final bound pose, slightly different from its X-ray structure pose (RMSD 2.1Å) bound to turkey β1-AR [PDB ID 6IBL (Lee et al., 2020)]. In its final bound pose, formoterol engages in polar H-bond and salt bridge interactions with S203^5.42, S207^5.46, D113^3.32, and N312^7.39. Formoterol (crystal pose in cyan) and the binding site residues (green) are illustrated in licorice representation. The receptor is illustrated in secondary structure representation (light golden yellow).

In 11 out of 15 simulations, prior to entering the binding pocket, formoterol stays in the region for approximately 25 ns, during which it engages in extensive H-bond interactions with many residues. Most importantly, the H-bond with D192^ECL2 through its O3 appears to destabilize the salt bridge between D192^ECL2 and K305^7.32 (Fig. 4B). At the same time, the oxygen (O4) of its methoxyphenyl tail forms another H-bond with Y316^7.43. In this extended orientation, formoterol also forms H-bonds through its carbonyl oxygen (O2) and hydroxyl group (O1) of the saligenin head with residues E180^ECL2 and N187^ECL2, respectively. The breakage of the salt bridge allows formoterol to slide down from this region into the binding pocket while forming a transient interaction between its positively charged amino group (N1) and D192^ECL2. As this interaction breaks, formoterol proceeds to the binding site underneath the hydrophobic lock between F193^ECL2 and Y308^7.35 (Fig. 4C). As the ligand passes down, this lock is broken, resulting in the phenyl ring of F193^ECL2 rotating toward the tail end of formoterol, forming a π-π interaction with its aromatic ring enabling the ligand to assume its final pose (Fig. 4D) similar to the crystal bound pose (PDB ID 6IBL) with an RMSD of 2.1Å.

On the other hand, in four out of 15 simulations, formoterol was seen engaging in extensive contacts with water in the region but much fewer interactions with the protein residues before entering the binding pocket. Interestingly, the hydrophobic lock (between F193^ECL2 and Y308^7.35) was not closed, and the phenyl ring of F193^ECL2 was seen rotated back toward Y174^ECL2, allowing formoterol to quickly enter into the binding pocket. Remarkably, in these simulations, the entire process of ligand access and reaching the binding pocket to assume its final bound pose took less than 25 ns. Most intriguingly, once formoterol settles into the binding pocket, the residues forming the hydrophobic lock come closer and remain in proximity for the rest of the simulation (see Supplemental Fig. 4).

During association, although both salmeterol and formoterol made initial interactions with residues from TMH1 and 2 before entering the pocket, there were some clear differences. The measured distance between TMH1 and TMH2, using the Cα carbon of Trp32^1.32 and Met96^2.66 was much larger for salmeterol than formoterol (Supplemental Fig. 5A), thanks to the bulkier tail group in salmeterol. There was a notable increase in this distance approximately 20 ns into the simulation, which was absent for formoterol (Supplemental Fig. 5A). Also, the change in the backbone RMSD of the entire protein was much larger during salmeterol association than for formoterol. This difference in the RMSD values (calculated using the backbone atoms) is ∼3 Å between the two compounds (Supplemental Fig. 5B). Specifically, the RMSD change comes from ECL2, which underwent significant conformational change (up to 7 Å) upon entry of salmeterol (Supplemental Fig. 6). Salmeterol’s large size and many rotatable bonds allow for the molecule to adopt a multitude of conformations that might not be suitable for entry into the binding site, and both the protein and salmeterol have to undergo conformational changes for the association process.

Salbutamol Enters the Binding Site Primarily through the Canonical Aqueous Path Described for Other Polar β2-AR Ligands

Unlike salmeterol and formoterol, the chemical polarity and hydrophilic nature of salbutamol do not allow the molecule to go deeper into the membrane (Fig. 1, B and C). Instead, salbutamol prefers to locate at the lipid headgroup-water interface, where it remains in contact with a significant number of water molecules (Supplemental Fig. 3). Also, as the solvation free energies for salbutamol in the aqueous bulk and at the membrane are quite comparable (Fig. 1B), spontaneous transfer of the ligand between these two phases can be expected. Therefore, we performed the WT-MetaD simulations for salbutamol in which the molecule was placed randomly around the receptor (a distance of at least 10 Å from any residue), with its starting location either at the aqueous bulk or in the membrane (time-average location and orientation, as determined by PMF; Fig. 1B). The FES for the receptor binding of salbutamol and the minimum energy path (black line) through which it gained access to the binding pocket predominantly is illustrated in Fig. 5.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 5.

The free energy surface and minimum energy path for salbutamol’s access and binding to the β2-adrenergic receptor. The primary binding path for salbutamol is through the common aqueous route, as determined by well tempered metadynamics simulations (minimum energy path, in bold black line). (A) In each simulation, salbutamol first entered the aqueous bulk before it established its first contact with the receptor. (B) Salbutamol first encountered K97^2.68 (through its saligenin hydroxyl groups). (C) Salbutamol slid along the extracellular vestibule, establishing interactions with F193^ECL2 through its t-butyl end and with D113^3.32 and N312^7.39 with its saligenin end. (D) Salbutamol is seen in its final bound pose, very similar to its crystal X-ray structure pose (RMSD of 2.25Å) bound to turkey β1-AR (PDB ID 2Y04) (Warne et al., 2011). In its final bound pose, salbutamol engages in polar interactions with S203^5.42, S207^5.46, D113^3.32, and N312^7.39. Salbutamol (blue) and the binding site residues (green) are illustrated in licorice representation. The receptor is illustrated in secondary structure representation (light pink).

While approaching the receptor from the aqueous bulk, salbutamol made its first contact with the side chain amino group of K97^2.68 and the carboxyl group of D300^ECL3 through its saligenin hydroxyl group and protonated ethanolamine amino group, respectively (Fig. 5B). These polar interactions drew salbutamol near the salt bridge between D192^ECL2 and K305^7.32 and enabled further contacts with D192^ECL2, F193^ECL2, and T195^ECL2. This was followed by destabilization and breaking of the salt bridge, allowing further entry of salbutamol into the binding pocket. The two hydroxyl groups of the saligenin head were seen engaging in H-bond interactions with D113^3.32 and N312^7.39 as the molecule glided down into the binding site (Fig. 5C). Subsequently, salbutamol assumed its final bound pose in which its two hydroxyl groups of saligenin head made H-bond interactions with side chain -OH groups of S203^5.42 and S207^5.46. Also, the -OH and -NH₃⁺ groups of the ethanolamine chain were engaged in H-bond/electrostatic interactions with the side chain amino group of N312^7.39 and the carboxyl group of D113^3.32, respectively (Fig. 5D).

Salbutamol Accesses the Binding Site through a Novel Polar Channel as a Secondary Route

Most surprisingly, in nearly half of the simulations (5 out of 12), salbutamol entered the binding pocket through a novel polar channel that opens up on the distal and opposite side of the common access route (Fig. 6, A and B). The entrance of this channel is located at the extracellular ends of TMH3/TMH4 and ECL2 and is lined by several polar residues, including N103^3.22, Y185^ECL2, R175^ECL2, E107^3.26, and Y174^ECL2. As soon as salbutamol approaches the channel opening, it quickly engages in polar interactions with these residues through its two hydroxyl groups of the saligenin head first. In addition, brief interactions were observed with ethanolamine and t-butyl groups of salbutamol with F104^3.23, S111^3.30, and Q170^4.62. Specifically, N103^3.22, R175^ECL2, and E107^3.26 make initial contacts with the ligand that draws the molecule inside the channel, during which Y174^ECL2 and Y185^ECL2 together act as an arm, re-engaging the two hydroxyl groups and pulling the ligand deeper into the pocket (Fig. 6, C and D). Subsequently, salbutamol interacts with T110^3.29, T195^ECL2, and Y199^5.38 as it slides deeper into the pocket before assuming its final pose akin to the crystal structure bound pose (PDB ID 2Y04) (Warne et al., 2011). Salbutamol also made a π-π interaction with the aromatic ring of Y199^5.38, the orientation of which seems to push salbutamol further into the pocket and prevent its exit back through the hole. In its final pose, the two hydroxyl groups of saligenin head made H-bond interactions with side chain -OH groups of S203^5.42 and S207^5.46. Also, the -OH and -NH₃⁺ groups of the ethanolamine chain were engaged in H-bond/electrostatic interactions with the side chain amino group of N312^7.39 and the carboxyl group of D113^3.32, respectively. To examine whether salmeterol or formoterol could also enter the binding pocket through this polar channel, we ran multiple SMD simulations (association) in which the ligands were kept near (within 5 Å) the channel entrance to facilitate the entry. Neither of the ligands entered the channel in any of the simulations, suggesting that this polar channel is unlikely to accommodate these larger lipophilic β2-AR agonists and thus signifies a unique access path for the small and polar salbutamol.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 6.

Salbutamol enters the β2-AR binding site through a novel polar channel. (A) From within the membrane view of the channel entrance (in yellow color and surface representation) displaying critical polar residues (in green color and licorice representation) that make immediate interactions with salbutamol (in blue color and licorice representation) as it enters. (B) The top view of the polar channel illustrating its position that is distal to the common access path lined by F193^ECL2. (C and D) The snapshots show the intermediate states of salbutamol before reaching the binding pocket in its final bound pose, as presented in Fig. 6D. The receptor is illustrated in secondary structure representation (white).

The Distances between the D192-K305 Salt Bridge and the Hydrophobic F193-Y308

Salmeterol seems to break the salt bridge interaction before entering the binding pocket in the majority of simulations. During association, salmeterol forms interactions with D192^ECL2, E306^7.33, and K305^7.32 during association that likely help break this interaction prior to entering between the transmembrane helices. As salmeterol slips into the binding pocket, it pushes aside the F193^ECL2-Y308^7.35 residues and breaks the hydrophobic lock to facilitate its entry. However, the rift of hydrophobic lock does not seem to follow a particular trend, as in many simulations, the residues are already far apart. For formoterol, a similar trend is seen with the salt bridge, in that the D192^ECL2-K305^7.33 is intact prior to association, and as formoterol gets closer to the binding site, it breaks these contacts (almost like a knife slicing through the bond). Formoterol forms contacts with D192^ECL2 during the association in various simulations and might attribute to the salt bridge breaking open. The F193^ECL2-Y308^7.35 interaction is similar for formoterol in its distance, and as formoterol settles in the binding site, it pushes apart the F193-Y308 interaction. Salbutamol also appears to break the D192^ECL2-K305^7.33 interaction in the aqueous pathway, but for the polar channel access, it slips underneath the F193^ECL2-Y308^7.35, and the salt bridge, for the most part, seems to be unaffected until salbutamol is in the binding site. Salbutamol tends to explore where it can enter easiest, whereas salmeterol and formoterol seem to need to break the salt bridge to enter because of their size and location.

Unbinding Pathways for the β2-AR Agonists: Salmeterol Prefers to Localize into the Membrane after Dissociation from the Receptor

To elucidate the unbinding mechanisms of the studied ligands, including the paths and critical residue interactions affecting the kinetics of the dissociation processes, we performed SMD simulations using the β2-AR-ligand complexes in their native membrane environment. For salmeterol, the ligand-bound crystal structure of β2-AR (PDB ID 6MXT) (Masureel et al., 2018) was used as the starting point. For both formoterol and salbutamol, the ligands were docked to the same β2-AR crystal structure, replacing salmeterol and using their respective β1-AR-bound crystal poses [PDB IDs 6IBL (Lee et al., 2020) and 2Y04 (Warne et al., 2011), respectively] as guidance for the docking simulations. After a 50-ns equilibration run, in which all ligands maintained stable contacts and remained in the binding pocket, for each ligand, three SMD runs of 50 ns were performed. The ligands were pulled out of the binding pocket with a minimal force in the direction of least resistance (see Materials and Methods section for details). The SMD simulations revealed unique dissociation paths from the binding site for each ligand (Fig. 7, A–C). Most interestingly, in two of three simulations, salmeterol dissociated from the binding pocket and eventually localized into the membrane (Fig. 7A). The bilayer depth and orientation of salmeterol were identical to those determined by the membrane partitioning simulations carried out without the receptor (Fig. 1, B and C). Unlike salmeterol, both formoterol and salbutamol dissociated toward the region around ECL2/TMH4/TMH5 and remained in the aqueous bulk above the receptor (Fig. 7, B and C).

• Download figure
• Open in new tab
• Download powerpoint

Fig. 7.

Unbinding paths and critical residue interactions along the dissociation pathways for salmeterol, formoterol, and salbutamol. (A–C) Dissociation paths for salmeterol, formoterol, and salbutamol, respectively. Salmeterol dissociates into the aqueous bulk and spends a transient period there before partitioning into the membrane near the receptor. Formoterol and salbutamol dissociate into the aqueous bulk surrounding the extracellular loop region. The start, intermediate, and final poses of the molecules are represented as density maps in red, gray, and blue colors, respectively. (D–F) Critical residues in contact with the ligands at their bound states and along the dissociation path. (G–I) Pie charts representing residues from various TMHs and ECL2 and the fraction of simulation time (percent occupancy) during which they are within 4 Å of the ligands.

When bound, all three ligands make conserved polar interactions (H-bonds) with the binding site residues S203^5.43, S207^5.46, D113^3.32, and N312^7.39 through their two hydroxyl groups of saligenin ring (one hydroxyl and one formyl amino group in the case of formoterol) and a hydroxyl and positively charged amino group of ethanolamine, respectively (Fig. 1A). As the ligands dissociated from the pocket, these interactions were broken at various time points for each ligand, indicating the specific strengths of those interactions attributed to the differences in their structures. Also, as the ligands drifted away from the site, they were in contact with a different set of residues from various TMHs and ECL2 along the way. The fraction of the total simulation time (50 ns) during which each residue was in contact with a ligand was given as percent occupancy in Fig. 7, G–I. For each ligand, as the simulation started, the amount of force used for pulling the ligand out of the pocket increased gradually as a function of time (Fig. 8A). It reached a maximum at which critical bonds were broken successively, leading to a decrease until the ligands were entirely dissociated. Ligands formed new interactions as they drifted away from the pocket, which caused momentary fluctuations in the required force, which was more pronounced for formoterol. The time points at which the individual polar interactions severed (bond length > 4 Å) were slightly different for the ligands.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 8.

Force profiles, critical electrostatic/H-bond interactions, and PMF profiles of salmeterol, formoterol, and salbutamol while unbinding from the β2-adrenergic receptor revealed by steered MD simulations. (A) Average force profile for each ligand calculated using three independent SMD simulations. The cumulative force applied increases until critical interactions get disrupted and then decreases to zero as bonds are broken successively. (B) The time evolution of the H-bond distances between the ligands and the binding site residues, indicating the strength and longevity of the interactions contributing to binding affinities of the ligands. (C) Salmeterol bound to the receptor in its crystal structure pose, showing all the critical residues involved in its dissociation path. (D) The PMF profile for each ligand was calculated by Jarzynski equality.

For salmeterol, the H-bond with S207^5.46 was the first one to break at approximately 8 ns, followed by S203^5.43, D113^3.32, and F193^ECL2 at around 12 ns. The H-bond between the positively charged ethanolamine amino group of salmeterol with N312^7.39 lasted until approximately 15 ns, at which point all the critical interactions with the binding site residues were broken. A similar trend was observed for formoterol, except that it made additional interactions with several other residues (Fig. 7, E and H) along the way. Interestingly, for salbutamol, the bonds with S203^5.43 and S207^5.46 broke at the same time (∼13 ns), followed by D113^3.32 and N312^7.39 at approximately 15 ns.

Notably, formoterol remained in contact with binding site residues for a longer time relative to the other two ligands. Specifically, the H-bond between N312^7.39 and the positively charged amine group on formoterol was intact until after 16 ns (Fig. 8, B and D). All the other H-bonds between binding site residues and formoterol begin breaking at 10 ns, starting with S203^5.42 and S207^5.46, and followed by D113^3.32 around 12 ns. Formoterol also made contacts with many other residues, including H296^6.58, Y308^7.35, A200^5.39, and Q197^5.36, during dissociation (Fig. 7, E and H). Similarly, salbutamol also interacted with residues T195^ECL2, A200^5.39, Y308^7.35, and F193^ECL2 during dissociation (Fig. 7, F and I). Furthermore, unbinding of the ligands disrupted the salt bridge between D192^ECL2 and K305^7.33 and the hydrophobic lock between F193^ECL2 and Y308^7.35 along the way, which occurred mostly after 10 ns into the simulations (Supplemental Fig. 10).

The free energy profile along the reaction coordinate for each ligand was constructed as the PMF from multiple dissociation SMD simulations using the JE (Jarzynski, 1997), as described in the Materials and Methods section. The PMF values obtained for salmeterol, formoterol, and salbutamol are 35.266 ± 0.796 kcal/mol, 39.859 ± 0.29462 kcal/mol, and 34.621±0.495 kcal/mol, respectively (Fig. 8D). For salmeterol, all the binding site interactions are broken after 15 ns, and the dip in the PMF curve is due to its interactions with the ECL2 residues (D192, F193, F194, and T195) and surrounding water molecules before completely dissociating from the receptor and eventually settling in the membrane environment near the receptor. The PMF profile for salbutamol shows a similar dip as salmeterol after all the binding site interactions were broken due to its interactions with ECL2/TMH5 residues (F193, T195, and A200) along the dissociation path. Once these transient interactions were broken, salbutamol entered the aqueous phase and remained there for the rest of the simulation. The PMF profile for formoterol was unique, which clearly revealed strong interactions with other residues after leaving the main binding site cleft. In summary, all the three studied ligands dissociated from the binding site around the same time. However, the PMF curves indicate formoterol as having a slightly higher binding free energy compared with salmeterol and salbutamol, which have approximately the same binding free energies (Fig. 8D).

Hydration of Ligands and Binding Site Residues during Dissociation Simulations

Desolvation (removal of water or dewetting) of a ligand and its binding site is an integral step before the receptor-ligand binding process (Abel et al., 2008). The opposite process, meaning solvation of the ligand and binding site, occurs when the ligand dissociates from the binding site. To assess the extent of solvation along the dissociation paths, we calculated the number of water contacts (number of water molecules that are within 3.5 Å from any heavy atom) for the studied ligands and the binding site residues (D113^3.32, S203^5.42, S207^5.46, F193^ECL2, F194^ECL2, T195^ECL2, D192^ECL2, N312^7.39). The plots of water contacts for the ligands and binding site against the simulation time showed unique characteristics for each ligand (Fig. 9A). As expected, during dissociation, each ligand’s water contacts increased as the ligand drifted away from the binding site. However, interestingly, the water contacts for salmeterol increased until 30 ns (from 6 to 17) but decreased afterward (from 16 to 10) as the ligand partitioned into the membrane. Formoterol and salbutamol had a similar trend in that the number of water contacts increased steadily as the ligands dissociated (from 4 to 15 for formoterol and from 5 to 10 for salbutamol) but attained steady levels afterward. However, the decrease was less pronounced for these ligands.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 9.

The number of water contacts as a measure of solvation (wetting or hydration) of the ligands (A) and binding site residues (B) during the dissociation of salmeterol, formoterol, and salbutamol. The water contacts were at a minimum at the start of the simulation (t = 0 ns) when the ligands were bound. The contacts increased gradually as the ligands drifted away from the binding site and the degree of increase was directly proportional to the ligand size in the order of salmeterol > formoterol > salbutamol. Notably, the water contacts decreased for salmeterol, indicating its relocation into the surrounding membrane. Similarly, the number of water contacts with the binding site residues increased gradually, converging to approximately the same number of contacts (∼28).

The water contacts of the binding site residues exhibited similar differences among the three ligands (Fig. 9B), clearly reflecting the ligand sizes. While bound, salbutamol had more water in the binding site than formoterol and salmeterol (∼22, ∼20, and ∼18, respectively). The smaller size of salbutamol compared with the other bulky ligands can explain this difference. Once the ligands drifted away from the binding site, the number of water contacts converged to similar values (26–28 water molecules) for all the ligands, indicating an absence of significant conformational change in the receptor after ligand dissociation (Mason et al., 2012).

Binding Free Energies Calculated by Funnel Metadynamics

We used FM (Limongelli et al., 2013; Raniolo and Limongelli, 2020) to elucidate all possible binding modes and calculate the absolute binding free energies of the three β2-AR agonists. For integral membrane proteins such as β2-AR, the molecular recognition process occurs at the interface between protein, membrane, and solvent molecules, including ions. It is critical to employ an appropriate method for large and flexible ligands such as salmeterol, explicit solvent molecules, membrane lipids surrounding the receptor with a buried binding site, and regions that undergo notable conformational changes (for example, ECL2). Although a priori knowledge of the ligand-binding mode is not required for FM simulations, we used the crystal bound pose as starting conformations as used in the dissociation SMD simulations described earlier. Also, the results obtained from the WT-MetaD simulations for elucidating the association mechanisms provided useful information on the proper placement of the funnel in FM simulations (see the Materials and Methods section).

For salmeterol, the funnel was placed such that the binding and unbinding states incorporate the ligand access through TMH1 and TMH7 from within the membrane. This step is critical, as salmeterol has been found to enter the binding pocket from within the membrane through the transmembrane helices in the association simulations. Figure 10 illustrates the binding free energy surface labeled A–F with the bound and multiple unbound states of the ligand (clustered) and placement of the funnel (G). The binding free energy was calculated from the energy difference between the unbound and bound region characterized by two collective variables (lp, the position of the ligand along the funnel axis, and ld, the orthogonal distance of the ligand from the funnel axis) with a correction for the restraining potential of the cylinder portion of the funnel. The bound state was characterized by lp and ld values in the range of 0.5–0.7 nm and 0.1–0.25 nm, respectively. The lp and ld values for the unbound state are 2.5–3.0 nm and 0.1 nm, respectively. The Hills potential deposited along the simulation time was monitored, and the near absence of any depositions indicated convergence of the free energy (Supplemental Fig. 9A). During the 900-ns simulation time, there were at least five crossing events between bound and unbound states (Supplemental Fig. 9B). After applying the correction for the restraining potential of the cylinder portion of the funnel, the calculated binding free energy (ΔG_b⁰) was −14.0 kcal/mol, which is reasonably in agreement with the published experimental value of −11.59 kcal/mol.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 10.

The free energy surface of salmeterol binding to β2-AR obtained from the funnel metadynamics simulation using lp, the position of the ligand along the funnel axis, and ld, the distance of the ligand from this axis, as collective variables. Two minimum energy paths, one passing through the aqueous bulk (green) and another through the transmembrane helices (red), reaching the membrane, were determined by the Minimum Energy Path Surface Analysis tool. (A–F) Multiple clusters representing the fully bound state (A) to several intermediate conformational states (B–F) were sampled by several recrossing between bound and unbound states during the entire 900-ns simulation. (G) The funnel placement was based on several association/dissociation simulations during which salmeterol was seen leaving the pocket either by aqueous routes or by a transmembrane route through TMH 1 and 7. The cone region of the funnel was defined by a vertex height Zcc of 3.0 nm from the origin and an α angle of 0.6 rad. The radius of the cylindrical portion of the funnel rcyl was 2 Å.

For formoterol, the FM simulation successfully captured the bound pose (similar to the docked pose) and unbound pose multiple times during the 450-ns simulation time (Supplemental Fig. 9B). The funnel was placed such that the simulation captures the ligand’s association through the aqueous access route. The bound region was characterized by the collective variables in the ranges of 0–0.4 nm for lp and 0–0.3 nm for ld. The values for lp and ld in the unbound region were 3 nm and 0 nm, respectively. The final calculated (see the Materials and Methods section for details) binding free energy (ΔG_b⁰) for formoterol was −10.14 kcal/mol, which is in good agreement with the published experimental value of −11.7 kcal/mol. Supplemental Fig. 7 illustrates the binding free energy surface, labeled A–F, with the bound and multiple unbound states of the ligand (clustered) and placement of the funnel (G).

For salbutamol, the placement of the funnel was to incorporate the dissociated state in the aqueous solvent. Over the course of 500 ns, there were only two recrossing events capturing bound and unbound states (Supplemental Fig. 9B). Supplemental Fig. 8 illustrates the binding free energy surface, labeled A–E, with the bound and multiple unbound states of the ligand (clustered). Considering a correction for the restraining potential of the funnel, the final calculated binding free energy (ΔG_b⁰) was −12.5 kcal mol⁻¹, which is comparable to the experimental data [Δ G_{b_exp} = −9.78 kcal mol⁻¹ (Baur et al., 2010), see detailed information in Materials and Methods]. The small discrepancy in binding free energy can be explained by fewer recrossing events captured during the simulation relative to formoterol and salmeterol. It is important to note that even though MD-based FM protocol has been shown to predict absolute binding free energies in good agreement with the experimental values, one should be cautious about inaccuracy in force-field parameters that may introduce error. Our simulations not only resulted in predicting the binding free energy in good agreement with the experimental values but also provided useful insights into the binding poses and conformations of the ligands in the unbound states.