Received 5 November 2015; accepted 24 January 2016; published 28 January 2016
Most smokers continue tobacco use because of their nicotine addiction despite of the well-known adverse health consequences. Of the more than 40 million smokers in the United States, more than 70% of smokers express a desire to quit smoking  . Only 3% of smokers who made a quit attempt in the past year were successful   .
Nicotine from inhaled smoke rapidly moves into the brain to exert its actions by binding to neuronal nicotinic acetylcholine receptors (nAChRs), which are members of the Cys-loop family of ligand-gated ion channels  . To date, 12 neuronal-type nAChR subunits (α2-α10, β2-β4) which can assemble into functional heteromeric or homomeric (α-subunit only) pentamers have been identified  -  . In the mammalian brain, the α4β2 subtypeisone of the most commonly found nAChRs   . It plays important roles in mediating the nicotine reward and addiction process and has been targeted for nicotine addiction treatment  . The pocket of α4β2 nAChR, to which various ligands can bind, is situated at the interface of the primary α4 subunit and the complementary β2 subunit   . Understanding the molecular interactions between α4β2 nAChR and its ligands could assist development of virtual screening of tobacco constituents with addiction potential through human α4β2 nAChR binding. The three-dimensional (3D) structure of the receptor and its dynamic molecular properties are the key for elucidation of the molecular interactions  . However, the 3D structure of a human α4β2 nAChR has yet to be experimentally determined. On the other hand, the 3D crystal structures of acetylcholine binding proteins (AChBPs) from Lymnaea stagnalis (Ls)  , Aplysia californica (Ac)  , and annelid Capitella teleta (Ct) have been experimentally determined   . These AChBPs are structurally similar to the extracellular ligand-binding domain (LBD) of human AChRs and thus serve as an important surrogate for the study of nAChRs. Based on the template of these AChBP X-ray structures, homology structures of the human α4β2 receptor have previously been constructed using the 3D crystal structure data for the AChBPs    , although this information has not been made available to the scientific community.
In this study, we constructed a 3D homology model of the extracellular domain of the human α4β2 nAChR (hereafter referred to as human α4β2 nAChR) using the X-ray crystal structure of Ct-AChBP that was determined at a high resolution  . The homology structure of human α4β2 nAChR was then further optimized using molecular dynamics (MD) simulations through the exploration of the accessible structural space for the complex   . We used the eleven compounds (including compounds contained in tobacco products such as nicotine) that were co-crystalized with nAChRs or AChBPs to elucidate the ligand binding interactions because there are α4β2 binding data for these eleven compounds. These eleven α4β2-ligand complex structures were obtained from molecular docking and further optimized by MD simulations. Examining the binding modes of these compounds in the α4β2 nAChR homology model revealed that hydrophobic interactions between ligands and V96, L97 and F151 of the α4 subunit and L111, F119 and F121 of the β2 subunit and the hydrogen bonds with S153 and W154 of the α4 subunit played key roles for the ligand binding.
The 3D structures of our human α4β2 nAChR homology model, its eleven ligand-bound complexes (supplementary materials which were made publically available), and the elucidated key interactions of these ligands in the α4β2binding pocket, could be used to develop a model for rapid screening of tobacco constituents with addiction potential through human α4β2 nAChR binding.
2. Material and Methods
2.1. Study Design and Workflow
The study design and workflow were shown in Figure 1. Briefly, the human α4β2 nAChR (target) sequences were first searched in the PDB using BLAST to select a highly homologous protein (template). The 3D structure of human α4β2 nAChR was then built using homology modeling. After docking the ligand co-crystalized with the template protein, the initial α4β2-ligand complex was then optimized with MD simulations. After MD optimization, the ligand was removed from the complex. The remaining optimized 3D structure of human α4β2
Figure 1. Study design and overall workflow. BLAST was used to search the human α4β2 nAChR sequence against PDB to select a template sequence for homology modeling. The initial homology model was then optimized using MD simulations. Eleven structurally diverse ligands PDB were docked into the optimizedα4β2 and followed by MD simulations.
nAChR was used to define a binding grid for subsequent docking analyses. Eleven structurally diverse α4β2 ligands were selected from PDB database and their 3D structures were built and optimized. The optimized ligands were then docked into the α4β2 binding pocket. The α4β2-ligand complexes obtained from docking were further optimized through MD simulations. Finally, the protein-ligand interactions were analyzed by examining the trajectory of MD simulations.
2.2. Homology Modeling
The tertiary structure of a protein is better conserved than its amino acid sequence   . Homology modeling is widely used to construct atomic-resolution models for proteins without experimentally determined structures from their amino acid sequences and experimental 3D structures of related homologous proteins   . The homology model of human α4β2 nAChR was built using Maestro 9.6 in the Schrödinger Suite  . The extracellular ligand-binding domain sequences of the human α4and β2 nAChR subunits were downloaded from UniProt database  (α4: P43681(29-242), ACHA4_HUMAN and β2: P17787(26-233), ACHB2_HUMAN). The 3D structure of Ct-AChBPs (UniprotID: I6L8L2) was downloaded from PDB (PDB ID: 4AFH)  as the template structure. The template sequence was identified through searching the target human α4β2 sequences in PDB using the BLAST homology searching with default parameters in the Prime 3.4. The target and the template sequences were aligned using the ClustralW method employed in Prime 3.4. The hetero-multimer model was used to predict the secondary structure of human α4β2 dimer. The ligand co-crystalized with the template protein was removed before homology modeling.
2.3. Molecular Docking
There are many structure-activity relationship techniques such as pharmacophore modeling  -  , comparative molecular field analysis   , machine learning methods  such as classification tree model  and Decision Forest models  -  based on molecular descriptors  that are calculated using the algorithm developed for the expert systems of structure elucidation  -  , and molecular docking    . Molecular docking was used in our previous studies of interactions between ligands and estrogen receptors   for prediction of estrogenic activity  and endocrine disruptions   . A similar docking protocol was adopted in this study. The 3D structure of α4β2 nAChR obtained from homology modeling followed by MD optimizations was prepared using the “protein preparation” tool in the Schrödinger suite  . The protein structures were optimized to reach the converged RMSD of 0.3 angstrom (Å) under OPLS_2005 force field. The grid box was defined using “receptor Grid Generation” tool also in the Schrödinger suite. The grid was set as a 15 × 15 × 15 (x × y × z, Å) cubic box around the ligand thus comprising the entire binding pocket.
The two-dimensional (2D) structures of the eleven compounds co-crystalized with AChBPs in PDB are shown in Figure 2. Hereafter, the 3-letter symbols given in Figure 2 will be used to indicate the ligands for simplicity; their names and experimental activity data are also given in Supplementary Table S1. Their 3D structures of the ligands were obtained from the corresponding PDB files of the bound AChBP complexes and further optimized usingLigPrep2.8 tool and OPLS_2005 force field in Schrödinger suite. Epik 2.6 tool of Schrödinger suite was used to generate the ligands ionization states at pH 7.0 ± 2.0. The optimized ligands were docked into the docking grid in the 3D structure of α4β2 using Glide 6.1 tool in the Schrödinger suite with standard precision (SP). The binding poses with the top glide scores were selected and the corresponding complexes were output for subsequent MD simulations.
2.4. MD Simulations
The 3D structure of the human α4β2 homology model and the elevenα4β2-ligand complexes output from molecular docking were optimized using MD simulations. The MD simulations were carried out through AMBER11  using the Amber ff99SB force field  and general Amber force field (GAFF)  for the receptor and ligands, respectively. AM1-BCC charges  were added for the ligands using Antechamber. The topology and parameters of the complexes were constructed by LEAP program  . Each complex was solvated with TIP3P water molecules with a 10.0 Å truncated octahedron periodic box. Sodium ions were added to maintain charge neutrality of the systems. The SHAKE algorithm  was used to constrain the bonds involving hydrogen. The Particle-mesh Ewald (PME) algorithm  was used to deal with the long-range interactions with a 8.0 Å non-bonded interactions cutoff. Prior to MD production simulations, a total of 10,000 steps of minimizations were applied to the entire model system including solvent molecules, followed by 2000 steps of steepest-descent, and subsequently 8000 steps of conjugated gradient minimizations. Two hundred ps (picoseconds) of equilibrations were applied to the system to reach the stable state. The molecular systems were heated in NVT ensemble from 0 to 300 K gradually over 200 ps with weak harmonic potential solutes restrained. After equilibration, 50 and 10 ns (nanoseconds) NPT ensemble MD simulations were performed for the α4β2 homology model and its eleven ligand-bound complexes, respectively, at 1 atm pressure, temperature of 300 K and a time step of 2 fs (fem to seconds). The coordinates of the systems were saved every 10 ps for subsequent trajectory analyses.
Figure 2. Chemical structures of the eleven compounds in the AChBPs complexes collected from PDB and their names used in this paper.
3. Results and Discussion
3.1. Optimized 3D Structure of Human α4β2 nAChR Homology Model
The sequences of the α4 and β2 subunits of the human nAChRs were queried using the BLAST algorithm to select a template structure in the protein data bank (PDB) to be used for the homology modeling of the human α4β2 nAChR. The α4 subunit sequence was aligned with 811 sequences in PDB (see Supplementary Table S2). The β2 subunit sequence was aligned with 765 sequences in PDB (see Supplementary Table S3). Recent studies have demonstrated that 3D structures can be constructed using homology modeling based on template structures of proteins with 25% or higher sequence identity to the target proteins   . The 3D structure of marine annelid Ct-AChBP with PDB ID 4B5D at a resolution of 2.2 Å has been successfully used previously as the template for the homology modeling of α4β2 nAChR to study the structural determinants of ligand recognition inhuman α4β2nAChR  . The crystal structure of Ct-AChBP with a higher resolution of 1.88 Å had been determined recently and deposited in PDB (ID: 4AFH). The α4 and β2 subunit sequences were found to have 31% and 28% sequence identities, respectively, with the 4AFH structure. Therefore, this was used as the template to construct the homology model of the human α4β2 nAChR.
The final alignment between the template and the target with the secondary structure prediction is shown in Figure 3 (the first three residues of α4 subunit were not included in the homology modeling because they could not be aligned to the template sequence).Most of the secondary structures including the disulphide bridges aligned well between the template and the target. The core residues located at the dimer interface, include residues Y98, W154 and Y195from the α4 subunit and residues W57, V111, F119 and L121 from theβ2 subunit. These core residues played key roles in α4β2-ligand interactions, which is consistent with the interactions between the human α4β2 nAChR and the nicotine addiction treatment drug varenicline derived from analysis using the structure-guided mutagenesis  .
The 3D structure of the human α4β2 dimer shown in Figure 4(a) was constructed by homology modeling using the structure of Ct-AChBP (Figure 4(b)) as the template. The homology model superimposed well with the template structure as depicted in Figure 4(c). Consistent with our expectation, the α helices and β sheets in both the template and target structures were more conserved compared to the loops. Our homology model was
Figure 3. Sequence alignment result between the human α4β2 nAChR and the template Ct-AChBP. The predicted secondary structure from Ct-AChBP is also shown: the tube represents the helix; arrows indicate strands; and dotted lines depict loops. The identical residues are highlighted in colors according to the property of the corresponding query residue in the alignment to depict the location of charged (purple: positive charge; red: negative charge), aromatics (brown), polar (cyan), hydrophobic (green) and special (yellow: cysteine; grey: proline) query residues. The residues that interact with the ligands are marked with black boxes.
(a) (b) (c)
Figure 4. Homology model of the human α4β2 nAChR. Template structure Ct-AChBP color coded in red (a) and the homology model is shown in solid blue ribbons for α4 and green ribbons for β2 (b). Superimposition of the homology model with the template is depicted in (c).
compared with the homology model constructed using a template structure of the same protein with a slightly lower resolution  . The overall structure and key structural features were similar.
More often than not, initial 3D structures from homology modeling were further optimized using MD simulations  . MD simulations are a technique widely applied to study protein folding, protein dynamics and protein-ligand interactions    -  . Considering the ligand impact on the residues conformations in the binding pocket  of the α4β2 homology model, MD simulations were performed in the presence of the ligand (LOB) (co-crystallized in the Ct-AChBP template protein). The root mean square deviation (RMSD) values of the optimized human α4β2-LOB complex structure were calculated to characterize the dynamic property of the molecular system during the MD simulations and were plotted in Figure 5(a). The molecular system reached stable states during the 50 ns simulations. The RSMD value of human α4β2 backbone atoms was stabilized at 3.5 Å after 30 ns simulations and the ligand fluctuated stably from 43 ns with a RMSD of 2.0 Å. The optimized human α4β2 homology model is available upon request.
The root-mean-square fluctuations (RMSF) of residues of the human α4β2 homology model during the MD simulations were calculated to characterize the mobility of the individual residues and plotted in Figure 5(b). As expected, the N-terminal alpha helix and the 10 beta strands are more rigid than the loops and, thus, had smaller RMSF values. The peaks in Figure 5(b) indicated larger RMSF values for the residues in the loops: 130 - 140, 170 - 180, 260 - 270, 280 - 290, 340 - 350 and 400 - 410.
The Ramachandran plot (Figure 5(c)) showed that 99.3% of the residues in the optimized α4β2 structure were in the favorable regions-higher than the percentage of residues (96.2%, see Figure 5(d)) in the initial homology model without MD simulations―thus indicating that the homology model of human α4β2 was optimized to an energy stable state.
To examine the structural changes of the α4β2 homology model following the MD simulations, the structures before and after the optimization were superimposed together and depicted in Figure 5(e). The two structures superimposed well with each other, especially for the β-sandwich cores, providing evidence that the homology model is a reasonably stable structure. On the other hand, much larger changes in the loop regions (some near the ligand binding site) were observed during the simulations points towards a better characterized ligand binding pocket for the α4β2. This indicates that optimization using MD simulations is a necessary step for the structural refinement of the human α4β2 homology model.
To examine the binding interactions between the human α4β2 homology model and the cognate ligand in the template, LOB, the zoomed-in view of the receptor ligand binding pocket is shown in Figure 5(f). The tertiary amine of LOB was stabilized through the π-cation interactions with the conserved aromatic residues W154 and Y195. The benzene ring of LOB also interacted with the piperidine ring of residue W154through π-π interactions. The receptor-ligand hydrophobic interactions mainly involved residues V96 and L97 from the α4 subunit and V111, F119 and L121 from the β2 subunit. The binding interactions observed for the human α4β2 homology model were similar to those between LOB and Ct-AChBP elucidated from the X-ray structure  .
3.2. Ligands Docking
Site Map  was used to generate the potential ligand binding site of the extracellular domains of human
(a) (b)(c) (d)(e) (f)
Figure 5. Optimization of the homology model by MD simulations. The RMSD values (y-axis) for atoms of protein backbone (blue line) and ligand (red line) during the optimization (x-axis) were plotted in (a). The RMSF values (y-axis) of individual residues (x-axis) were given in (b). Ramachandran plots of the optimized and initial homology models were given in (c) and (d), respectively. The homology models before (blue ribbon) and after (red ribbon) the MD simulations were superimposed (e). The binding pocket was illustrated by a zoomed surface representation (f). The core residues were color coded in blue. The arrows indicated the hydrophobic moment (residues colored in red). The red lines depicted the π-π interaction and the purple lines for the π-cation interaction.
α4β2n AChR homology model. The eleven ligands shown in Figure 2 were docked in the receptor model. For each ligand, the lowest docking score and its corresponding docking pose in the receptor were selected and output for subsequent analysis. The docking poses of all eleven ligands in the human α4β2 nAChR homology model are shown in Figure 6 and Supplementary Figure S1. The ligand binding site composed of the α4 subunit residues V96, L97, Y98, S153, W154, Y195 and Y202 and the β2 subunit residues W57, V111, F119, L121 and F157. The orientations of the eleven ligands were similar i.e. the hydrophobic parts of all the ligands extended into the major surface of the binding pocket. Aromatic and hydrophobic interactions were the major contributions to the binding.
3.3. Ligand Binding Interactions
To elucidate the binding modes of the eleven ligands in the human α4β2 receptor, the complexes of α4β2 bound with the ligands resulting from the molecular docking were again optimized using MD simulations   . MD simulations were carried out on the eleven human α4β2-ligand complexes obtained from molecular docking and the results were given in Figure 7. The RMSD values generated by the MD simulations for the proteins (Figure 7(a)) and the ligands (Figure 7(b)) indicated the complexes reached relatively stable states. In more detail, for the structures of α4β2 homology model bound with IM4, 09O, 09R, LOB and QMR ligands, the irreceptor RMSD values stabilized around 2.0 Å after 6 ns MD simulations. On the other hand, the structures of α4β2 homology model bound with NCT, 09P, 09Q, 09S, SW4 and C5E continued to show small fluctuations of the receptor (RMSD values of ~2.5 Å) after 6 ns MD simulations. Furthermore, as expected, the rigid ligands (C5E and QMR) remained highly stable as evidenced by a low RMSD of 0.3 Å during the whole simulations; the ligands with one rotatable bond (NCT and 09R) had slightly larger fluctuations of ~1 Å RMSD; the ligands with multiple rotatable bonds (e.g. 09Q, 09S, and SW4) had larger fluctuations in the simulations. In other words, the more rotatable bonds in a ligand, the larger the RMSD. Interestingly, LOB has six rotatable bonds but its conformation only changed slightly with RMSD values of about 0.5 Å during the simulations. This indicates that the fluctuation of a ligand in simulations depends not only on the ligand structure but also on its binding environment   .
To elucidate the binding interactions of the eleven ligands, the corresponding receptor-ligand complexes from the final frame of the MD simulations were examined. The binding mode of NCT was shown in Figure 8(a). On the α4 subunit side, NCT had extensive binding with residues Y98, S153, W154 and Y202 through hydrogen bond, π-π and π-cation interactions. It also bound with the hydrophobic residues V96, L97 and F151 through hydrophobic interactions. On the β2 subunit side, the ligand bound with W57 and L121 through π-sigma interactions and hydrophobic interactions. The receptor-ligand interactions observed for the human α4β2 nAChR are consistent with the ones found in AChBP crystal complex of the same ligand  .
The binding interactions between α4β2 and IM4 were shown in Figure 8(b). The pyridine ring of IM4 interacted with α4β2 through the π-π interactions with residue Y195 of the α4 subunit and the hydrophobic interact-
Figure 6. Docking poses of the eleven ligands in the binding site of human α4β2.
Figure 7. Trajectory analyses of MD simulations. The RMSD values of the protein backbones (a) and the ligands (b) were plotted again the time. The complexes are indicated by the ligand abbreviations using different colors shown in the legends.
tions with residues L111, L119 and F172 of the β2 subunit. The binding mode of IM4 in the human α4β2 elucidated in this study was similar to the one observed in the crystal structure of the Ls-AChBP complex bound with IM4  .
Ligands 09O, 09P, 09Q, 09R and 09S share a structural fragment in which a diazepane and a pyridine are connected by a single bond. Their binding interactions with the human α4β2 homology model are depicted in Figures 8(c)-(g), respectively. Although the binding modes of these ligands in the receptor were not the same due to the structural elements attached to the common fragment, the binding interactions formed were similar. These interactions were hydrogen bond, hydrophobic, π-π and π-cation interactions. The π-π and π-cation interactions mainly occurred between the diazepane and pyridine of the ligands and the residues W154, Y195, Y202, V96, L97 and Y98 of the α4 subunit and the residues W57, L121, F172 and F119 of the β2 subunit. These binding interactions as elucidated from our MD simulations are consistent with the ones identified in the AChBP X-ray complex structures bound with the respective ligands  .
The ligand binding mode ofSW4is illustrated in Figure 8(h). SW4 forms a hydrogen bond with residue Y196 of the α4 subunit and had π-cation and hydrophobic interactions with residues V111, F119 and F121 of the β2 subunit. These binding interactions were observed in the x-ray crystal structure of the Ct-AChBP bound with SW4  .
Figure 8. Ligands binding modes. Binding pocket surfaces are colored in the interpolated charge. The core residues are in cyan and the arrows indicated the hydrophobic moments (residues were labelled in red).The red lines depicted the π-π interaction and the purple lines for the π-cation interaction. The ligands in the complexes: NCT (a); IM4 (b); 09O (c); 09P (d); 09Q (e); 09R (f); 09S (g); LOB (h); QMR (i); SW4 (j) and C5E (k).
Figure 8(i) depicts the binding mode ofC5E.The ligand bound with the human α4β2 through hydrogen bonds, π-cation interactions, π-sigma interactions and hydrophobic interactions with the residues V96, L97, F151, S153, W154 and Y202 of the α4 subunit. It also has π-π interactions and hydrophobic interactions with the residues W57 and L121 of the β2 subunit. The binding interaction was confirmed in the x-ray structure of Ac-AChBP complex bound with C5E  .
The binding mode of QMR is shown in Figure 8(j). This highly rigid ligand bound to the binding pocket through π-π interactions, π-cation interactions, π-sigma interactions and hydrophobic interactions with residues V96, L97, F151, Y195 and Y202 of the α4 subunit and the side chains of residues W57, F119 and F121 of the β2 subunit. These binding interactions were in a good agreement with the Ac-AChBP complex bound with the same ligand  .
The binding mode of LOB is given in Figure 8(k). LOB interacted with α4β2 through π-cation interactions with residues W154, Y195 and Y202 and π-π interactions with W154 of the α4 subunit. It also had hydrophobic interactions with residues V96, L97 and F151 of the α4 subunit and residues V111, F119 and L121of the β2 subunit. These molecular interactions were similar to those observed in x-ray structure of LOB-CtAChBP complex  .
The binding interactions between the human α4β2 homology model and the eleven ligands were summarized in the Table 1. In general, the hydrophobic interactions mainly involved residues V96, L97 and F151 of the α4 subunit and residues L111, F119 and F121 of the β2 subunit. The π-π, π-cation, and π-sigma interactions were mainly from residues Y98, W154, Y195 and Y202 of the α4 subunit and from residue W57 of the β2 subunit.
Table 1. Summary of predicted binding interactions.
The hydrogen bonds were formed between the ligands and residues S153 and W154 of theα4 subunit. It is worth noting that residues W154 and W57 interacted with most of the eleven ligands and played key roles in the binding as confirmed in the X-ray structures of the AChBPs bound with the respective ligands   . It is worth noting that in this study, we focus on the binding of these ligands to the human α4β2 nAChR i.e. the pre-requisite step for chemicals to elicit subsequent pharmacological actions on the body-indeed, the ligands considered in this study display different pharmacology (various degree of activation and desensitization) upon binding to the human α4β2 nAChR. However, this aspect is out of the scope of our study.
The human α4β2 nAChR plays an important role in nicotine addiction and could be used for screening tobacco constituents with addiction potential. Determining the 3D structure of the human α4β2 nAChR and elucidating its ligand binding interactions would help with screening for tobacco constituents with nicotine receptor binding activity and addictive properties. This type of screening will result in a rank order of tobacco constituents that have addiction potential and may inform regulation of tobacco products. As there are no 3D structures of human α4β2 nAChR either alone or bound with ligands available to researchers today, we constructed a homology model of the human α4β2 nAChR based on the X-ray structure of the Ct-AChBP complex bound with LOB as an alternative approach for gaining insights into ligand binding characteristics of the human α4β2 nAChR. Based on the homology model of the human α4β2 nAChR, we further characterized the ligand binding interactions of human α4β2 nAChR by docking eleven ligands to the receptor followed by MD simulations. Our results provide molecular details into the ligand binding interactions with human α4β2 nAChR homology model and a modeling method for identifying nicotine receptor binding compounds that are potentially addictive.
This project was funded by the Center for Tobacco Products at the U.S. Food and Drug Administration. This research was supported in part by an appointment to the Research Participation Program at the National Center for Toxicological Research (Hui Wen Ng and Hao Ye) administered by the Oak Ridge Institute for Science and Education through an interagency agreement between the U.S. Department of Energy and the U.S. Food and Drug Administration. This publication represents the views of the author(s) and does not represent CTP position or policy.
H.H. and M.O. conceived and designed the experiments; M.S. and H.W.N. performed the experiments; M.S., H.W.N., H.L., H.Y. and W.G. analyzed the data; M.S., M.S., W.T. and H.H. wrote the paper.
Conflicts of Interest
The authors declare no conflict of interest.
Table S1. Information on the 11 ligands used in this study.
Table S2. Results of search a4 subunit on the PDB. ID consists of PDB ID (the first 4 letters) and the subunit (the last letter).
Table S3. Results of search b2 subunit on the PDB. ID consists of PDB ID (the first 4 letters) and the subunit (the last letter).
Figure S1. Docking poses of the 11 compounds in the homology model of human α4β2 nACHR.