Abstract
Some natural alkaloids from medicinal plants, such as yohimbine and its derivatives, have been reported with adrenoceptor (AR) α2 subtypes inhibiting activity. In trying to address the possible mechanism of the action, a set of homology models of AR α2 was built based on MOE. After that, docking and molecular dynamics methods were used to investigate the binding modes of yohimbine and its 2 derivatives in the active pocket of adrenoceptor α2 subtype A, B, and C. The key interactions between the 3 ligands and the 3 receptors were mapped. Binding mode analysis presents a strong identity in the key residues in each subtype. Only a few differences play the key role in modulating selectivity of yohimbine and its derivatives. These results can guide the design of new selective AR α2 inhibitors.
1. Introduction
Adrenoceptor (AR) is a signal receptor of epinephrine and norepinephrine and is mainly distributed in myocardial, vascular, and nervous system. Till now, 9 AR subtypes have been reported. They are classified into subtypes α1, α2 and β according to the pharmacological functions and sequences. α1 receptors have subtype A, B, and D. α2 receptors have subtypes A, B, and C. β receptors have β1, β2, β3 [1, 2].
Since discovered in 1948, AR family has always been one of the most important targets for drug discovery. Both the agonist and the antagonist on different AR subtypes can be leading compounds with pharmacological effects, particularly the significant cardiovascular effect. Looking for the candidate with high selectivity among the AR subtypes is the emphasis in drug development. The α2 ARs are widely distributed throughout the peripheral and CNS. Agonists acting at α2 ARs have analgesic properties following supraspinal [3], spinal [4], peripheral [5], and systemic administration [6].
Some natural products acting on the AR receptor have been found as the active ingredients of traditional Chinese medicine. All the compounds are alkaloids and are the AR blockers. According to the skeleton, the compounds can be divided into three categories: rauwolfia alkaloids, such as yohimbine; protoberberine, such as berberine; and aporphine, such as xylopine (Table 1).
Yohimbine (CAS: 146-48-5) was first isolated from bark of a West African tree Pausinystalia yohimbe (Rubiaceae) in 1945. The plant has a long history of human use as a stimulant and aphrodisiac, especially as a herbal supplement to improve erectile function. To date, Yohimbine was also found in Alchornea floribunda Müll. Arg. (Euphorbiaceae) [14], Rauvolfia vomitoria, Rauwolfia nitida, Rauwolfia serpentine, Rauvolfia yunnanensis, and others [15].
The pharmacological research on Yohimbine has been started since 1949 [16]. The reported activities include antidiuretic, antiadrenaline, mydriatic, serotonin antagonist, and others. The most remarkable activity is selective inhibition on ARs, high affinity for the AR α2 subtypes, and moderate affinity for the α1 subtypes. Binding affinities among the AR α2 subtypes are α2C (0.88 nM) > α2A (1.4 nM) > α2B (7.1 nM) [13, 17, 18]. Because of the high selectivity on AR subtypes, Yohimbine attracted more attention [11, 12]. Lots of pharmacological data have been reported, while the mechanism of the action is still unclear.
In this paper, three-dimensional homology models were built for each AR α2 subtype. By using molecule stimulation method, the binding modes of Yohimbine, Yohimbine monoglycine ester (YME), and Yohimbine alkyl amine (YAA) were analyzed (Figure 1). The mechanism of selectivity and the key residues were investigated.
(a)
(b)
(c)
2. Data and Methods
2.1. Homo Modeling
The AR receptor belongs to heterotrimeric guanine nucleotide-binding protein coupled receptors (GPCRs). The proteins in GPCR superfamily usually have low sequence identities but with high structure similarities. The binding sites usually locate among the seven transmembrane (TM) helices; together with three extracellular loops. GPCR receptors play important roles in drug discovery [19]. Due to the hardness in the crystallization of GPCRs, homology modeling is a main approach for the structure-based drug design of GPCRs [20].
The secondary structure assignments are annotated with color-coded bars above the sequences. Helixes are in red, strands are in yellow, 1–4 turns are in blue, and 1–5 turns are in green. The conserved residues were marked with dark color.
To date, our understanding of GPCR structure is mainly based on the high-resolution crystal protein structure. Computational modeling has also predicted that GPCRs share a seven membrane-spanning a-helix topology as a common structural property. Likewise, mutagenesis studies have tentatively identified individual transmembrane domains with specific roles in signal transduction function that are conserved throughout the GPCR family. Since 2007, several Homo AR β-subtypes were reported [21–24]. Among the crystal data, 2RH1, a complex of β2 subtype with agonist carazolol, is an ideal template for the modeling of AR family [25, 26].
The alignment errors increase rapidly when the sequence identity of two proteins is less than 30% [27]. For GPCR homology, the sequence identity is usually less than 25% [20]. Therefore, automated homology modeling of GPCRs is likely to result in more errors [28]. The errors in alignments cannot be justified by energy optimizations, molecular dynamics, or distance geometry refinements. The alignments become critical for the quality of the homology model. Multiple sequence alignments are adopted in order to enhance the alignment quality. We have used 8 AR subtypes in the sequence alignment. As shown in Figure 2, all the subtypes keep a high identity in the TM helices. For the binding site is exactly constituted by the TMs, the homo model may be of high reliability to estimate the interaction between the ligand and AR.
Homo models of AR α2A, B and C were constructed with the homology modeling module of the molecular operating environment (MOE) [29]. The AMBER99 force field was applied to the homology modeling. The carboxyl-terminal and amino-terminal modeling was disabled. The scoring method was GB/VI [30]. Ten models were built by MOE for each subtype, kept in a MOE database, and ranked by the GB/VI scores. The best model was kept for the following energy optimizing.
Partial charges for all atoms were calculated. Then, the energy was minimized by means of AMBER99 force field, and the ending condition of the calculation was the root mean square (RMS) of the two conformer energies being less than 0.1 kcal mol−1 Å−1. Furthermore, molecular mechanics (MD) simulations were applied, in MOE. The process consisted of two steps.(1)The MD simulations were performed with all Cα atoms being fixed. The temperature was increased from 0 K to 300 K within 60 ps (picoseconds); an NVT (constant volume/constant-temperature) molecular dynamic simulation of 100 ps at 300 K was carried out with a time-step length of 0.002 ps.(2)The atoms in the TM helixes were set flexible for conformation changes. The MD experiments simulated conformation changes in 1 ns timescale. The final model was obtained by energy minimization till an RMS of 0.05 kcal mol−1 Å−1 was reached.
2.2. Docking
Docking experiments were carried out to explore the binding modes of three antagonists, Yohimbine, YME, and YAA. All the ligands were built in MOE Builder. Their three-dimensional structures were energy minimized by MOE using MMFF94x force field. The initial poses were assigned by Flexible Alignment module of MOE, with the agonist carazolol in 2RH1 as the template. The 3 ligands were docked into 3 α2 subtypes, respectively. Totally 9 docking experiments were carried out.
Docking used the algorithm MOE Dock. The default settings were used. For each position, 700 iterations were generated to optimize the interactions at a location. A database of 25 complexes for each receptor with each ligand was generated. MMFF94 force field was used. Interactive forces between the ligand and the receptor include hydrogen bonding, aromatic interactions, and hydrophobic interactions. Of these, hydrogen bonding is the strongest interaction followed by aromatic interactions and last hydrophobic interactions. Therefore, the best complex was judged to be the complex with the greatest number of hydrogen bonds to side chain atoms, aromatic interactions, and placement in the helical bundle. No explicit water molecules were added during docking simulation. Solvation and entropic effects were not taken into account either.
The top 3 reasonable binding modes were picked, based upon compatibility with the reported mutation results and binding data, as well as energetic considerations, from the 10 resulting receptor-ligand complexes. A further refinement was done in MOE. The residues within 6 Å from the ligand were set to be flexible. Energy minimization was carried out by conjugated gradient minimization with the MMFF94x force field, until an RMSD of 0.1 kcal mol−1 Å−1 was reached.
2.3. Dynamic Simulations on Docking Results
The top 3 complexes for each ligand and receptor were then subjected to two-stage molecular dynamics simulations. The simulation was performed for 1 fs time step with 200 ps of gradually heating phase from 0 K to 300 K, followed by 2000 ps equilibrium (data collection) phase. The output databases contained 11000 entries collected during the equilibrium phase. The simulation used NVT parameters (holding constant moles, volume, and temperature) and calculated potential energy ( in kcal mol−1), temperature ( in kelvins), pressure ( in Kpa), total energy (, kinetic and potential in kcal mol−1), and enthalpy (, in kcal mol−1). The complex with the least RMSD fluctuation was considered as the best one.
3. Result and Discussion
Figure 3 shows the template protein 2RH1 and the three established homology models. There are little structural differences among the template and the models in the 7 TM helices. The three-dimensional structure of binding pocket, surrounded by TM3, TM5, TM6, and eLP2, is especially close to each other among the 3 models. Because of the lack of conserved residues in loop domain, there are much more difference in the 3 endocellular loops and 3 extracellular loops. Since the active pocket is located in the upper half of the TMs, the differences of the loops do not affect the study on the binding mode of ligands. Therefore, the model can be used for the following molecular docking studies.
The 3 ligands have been docked to the homology structures of the 3 AR subtypes A, B, and C, respectively. For each combination, the most reasonable result was chosen as the start conformation for a further molecule dynamics optimization. The RMSD trajectories generated from the MD simulation were shown in Figure 4. The final model was obtained by averaging the structures from the 1.5~2 ns trajectories of MD simulations. The model energy was minimized till an RMS of 0.05 kcal mol−1 Å−1 was reached.
Figure 5 shows the binding mode of the Yohimbine in the α2A three-dimensional structure. The 3 ligands’ 9 binding modes are depicted in Figures 6 and 7. The three-dimensional views are listed in Figure 6 and the ligand interaction maps are listed in the Figure 7. All the antagonists have the same orientation. The hydrophobic ends point to the opening of the pocket, and interact, through arene-H interaction, with Phe, Leu, or Ile residues on TM6 and eLP3. The polar ends face the bottom of the pocket; hydrogen bonds are formed with Tyr5.38, Ser5.42, and Ser5.46. The “hot” residues were identified by summarizing the interactions between the 3 ligands and the 3 α2 subtypes by 2D contact diagram of MOE (Figure 7, Table 2).
Table 2 summarized the interaction, hotspot residues, binding free energy, and reported Ki [11]. As can be seen, the binding mode of yohimbine is very close to that of YME and YAA. The interaction between the 3 ligands and the AR2 receptors mainly includes 4 aspects: (a) hydrogen bonds between oxygen atom 3 and 4 with polar residues Y5.38 and S5.46; (b) nitrogen atom 1 as a hydrogen bond donor formed a hydrogen bond with Y6.55; (c) ring A and the ring B contact with the aromatic/aliphatic residues, such as F6.51, I/L5.32, by arene-H interaction; (d) S5.42 may account for the selective inhibition of yohimbine and its derivatives.
As common structural characteristics through the three AR α2 subtypes, Y5.38, S5.42, S5.46, I/L5.32, and F6.51 are “hot” residues, the most interactive with the ligands (Table 2). The hydrogen bond interactions are the main contributors to binding affinities, which mainly located at the bottom of the AR binding pocket. The aromatic/aliphatic interactions are located at the opening of the pocket, close to the solvent. This area, the ligand tends to form arene-H bindings and hydrophobic interactions with “hot” residues, such as F6.51 and I/L5.32. Therefore, a good ligand should have at least one hydroxyl group at one end and two aromatic rings at the other end.
The root of the selectivity of the yohimbine and its derivatives may be the hydrogen bonds on S5.42. As shown in Table 2, the complexes with a hydrogen bond on S5.42 usually have a relatively higher Ki value. That also presents an explanation on yohimbine higher selectivity on AR α2C.
All the active pockets of AR α2 were formed by TM3, TM5, TM6, TM7, and eLP2. The pocket can be divided into two compartments (Figure 8). Compartment I is hydrophobic and close to the entrance of the pocket. Most residues in compartment I are aromatic/aliphatic, such as Phe, Leu, Ile and others. The most important of them are Phe6.51 and the Ile190 (α2A, eLP2), Leu166 (α2B, eLP2), and Leu204 (α2C, eLP2). Compartment II is a hydrophilic bottom of the binding pocket and has more hydrogen bonding opportunities, hosting “hot” residues as well: Tyr5.38, Ser5.42 and Ser5.46 on TM5.
Consequently, as to the structural characteristics of yohimbine and its derivatives, the 3 and 4 oxygen atoms and A, B rings (Figure 1) play the key role in the contact with ARs. The two oxygen atoms form “W” type 4 hydrogen bonds with Tyr5.38, Ser5.42, and Ser5.46, which all locate on TM5. The hydrogen bonds on Tyr5.38 and Ser5.46 are relatively stronger than those on Ser5.42. The two nitrogen atoms in the scaffold of yohimbine also interact with the residues on TM6 in some situations more or less. No. 1 nitrogen atom usually forms hydrogen bond with Tyr6.55; thus, it contributes more for the three ligands’ affinity on AR α2. The hydrogen bonds between O-3 and S5.42 may contribute more in the formation of the selectivity of the ligands.
4. Conclusions
A set of homology models was built for AR α2 subtypes, based on multisequence alignments and secondary structure analyses. The 3 ligands were docked into 3 AR α2 subtype models. Molecular dynamics was carried to find the stable binding mode. The “hot” residues and the mechanism of the action are identified by comparing the 9 ligand-receptor complexes. These discoveries can be important for new AR α2 selective antagonist design.
Acknowledgments
The authors would like to thank Dr. Weiyu Lin for assistance in MOE SVL programing. The work in this paper was supported by the Ministry of Science and Technology of China, Special Project SB2007FY020.