| HOME | HELP | FEEDBACK | SUBSCRIPTIONS | ARCHIVE | SEARCH | TABLE OF CONTENTS |
7 Nicotinic Receptor Ligand Binding Domain Revealed by Molecular Dynamics Simulation
* Howard Hughes Medical Institute, Department of Chemistry and Biochemistry, University of California San Diego, La Jolla, California 92093; Receptor Biology Laboratory, Department of Physiology and Biophysics, Mayo Foundation, Rochester, Minnesota 55905; and Department of Pharmacology, University of California San Diego, La Jolla, California 92093
Correspondence: Address reprint requests to Richard H. Henchman, 9500 Gilman Dr., Mail Code 0365, La Jolla, CA 92093. Tel.: 858-822-1469; Fax: 858-534-7042; E-mail: rhenchma{at}mccammon.ucsd.edu.
 |
ABSTRACT |
|---|
 TOP ABSTRACT INTRODUCTIONMETHODRESULTSDISCUSSIONCONCLUSIONACKNOWLEDGEMENTSREFERENCES |
|---|
7 nicotinic receptor is constructed based on the acetylcholine-binding protein crystal structure. This structure is refined in a 10 ns molecular dynamics simulation. The modeled structure proves fairly resilient, with no significant changes at the secondary or tertiary structural levels. The hypothesis that the acetylcholine-binding protein template is in the activated or desensitized state, and the absence of a bound agonist in the simulation suggests that the structure may also be relaxing from this state to the activatable state. Candidate motions that take place involve not only the side chains of residues lining the binding sites, but also the subunit positions that determine the overall shape of the receptor. In particular, two nonadjacent subunits move outward, whereas their partners counterclockwise to them move inward, leading to a marginally wider interface between themselves and an overall asymmetric structure. This in turn affects the binding sites, producing two that are more open and characterized by distinct side-chain conformations of W54 and L118, although motions of the side chains of all residues in every binding site still contribute to a reduction in binding site size, especially the outward motion of W148, which hinders acetylcholine binding. The Cys loop at the membrane interface also displays some flexibility. Although the short simulation timescale is unlikely to sample adequately all the conformational states, the pattern of observed motions suggests how ligand binding may correlate with larger-scale subunit motions that would connect with the transmembrane region that controls the passage of ions. Furthermore, the shape of the asymmetry with binding sites of differing affinity for acetylcholine, characteristic of other nicotinic receptors, may be a natural property of the relaxed, activatable state of 7. |
INTRODUCTION |
|---|
|
TOPABSTRACTINTRODUCTIONMETHODRESULTSDISCUSSIONCONCLUSIONACKNOWLEDGEMENTSREFERENCES |
|---|
; Changeux and Edelstein, 1998, 2001; Corringer et al., 2000; Grutter and Changeux, 2001; Sine, 2002; Karlin, 2002; Hogg et al., 2003). The lack of a nAChR structure at atomic resolution has hindered the study of the allosteric mechanism, particularly for modeling methods. The main sources of direct structural information had been the three-dimensional shape of the receptor at 4.6 Å resolution (Miyazawa et al., 1999) and the NMR structure of the M2 helix (Montal and Opella, 2002). Secondary structure determinations had also been derived from the primary sequence (Ortells, 1997; Le Novère et al., 1999). This level of structural data led to simulation studies only of the binding site (Morreale et al., 2002) or the M2 helices in the transmembrane region (Ortells et al., 1997; Tikhonov and Zhorov, 1998; Adcock et al., 2000). However, the recently determined crystal structure of the distantly homologous acetylcholine-binding protein (AChBP) (Brejc et al., 2001) has provided a useful stepping stone forward for such modeling by serving as a suitable template for modeling the nAChR LBD (Le Novère et al., 2002; Schapira et al., 2002; Molles et al., 2002; Sine et al., 2002).
Given the availability of a structure at atomic resolution, molecular dynamics (MD) simulations have proved to be an effective means of determining the structure's dynamic nature, albeit at still very short timescales. In this work, we construct a homology model of the human 7 LBD in a similar manner to that used recently to make a muscle nAChR model (Sine et al., 2002) and examine its stability and dynamic properties in a 10 ns molecular dynamics computer simulation. The 7 nAChR is a suitable choice because it is a homopentamer, the simplest form of nAChR, and its higher sequence identity (Smit et al., 2001) and relatively old ancestry make it a closer mammalian homolog to AChBP, leading to a more reliable model. Physiologically, 7 is widely expressed in the central and peripheral nervous systems and plays a role in many calcium-dependent processes by regulating calcium flow (Hogg et al., 2003).
The MD simulation provides a means to further relax the homology model. Many of the motions seen in the MD simulation would be expected to be optimizations resulting from imperfections in the homology model. The simulation also provides opportunities to observe other types of motion that may correspond with real allosteric motions from the active or desensitized states to the activatable state. Evidence suggests that the AChBP template is equivalent in conformation to the active or desensitized states of the LBD (Grutter and Changeux, 2001). This evidence includes the high affinity of AChBP for nicotinic ligands (Smit et al., 2001) and the closer resemblance to AChBP of the Torpedo nAChR electron microscopy structure with ACh bound than without, the main difference lying in small rotations of the -subunits (Unwin et al., 2002). With no ligand present in the simulation, the activated or desensitized 7 LBD would be expected to convert into the activatable state. The submicrosecond timescale of MD simulations is still too short to reproduce the full millisecond scale conformational motions. However, even for this short time, by examining the different motions occurring in the ACh binding sites, the Cys loop, and the subunits as a whole, various relationships between these types of motion may be hypothesized.
|
METHOD |
|---|
|
TOPABSTRACTINTRODUCTIONMETHODRESULTSDISCUSSIONCONCLUSIONACKNOWLEDGEMENTSREFERENCES |
|---|
7 ligand binding domain7 using version 6.0 of the program MODELLER (
ali and Blundell, 1993), together with spatial restraints provided by the AChBP structure (Brejc et al., 2001). The sequence alignment between AChBP and 7 monomers used to generate the homology model is shown in Fig. 1. To maintain complementarity between subunits at their interfaces, all five subunits were modeled simultaneously. The "patch" command in MODELLER was used to constrain the coordinates of C127 and C141, which form a disulfide bond in each subunit. The "refine1 mode" option was selected to generate the highest level of refinement using conjugate gradients coupled with simulated annealing and molecular dynamics. Modeling included all polar hydrogens to allow for main-chain hydrogen bonding but omitted nonpolar hydrogens. In this way, 100 different structures were produced and evaluated using the programs PROCHECK (Laskowski et al., 1993) and PROFILES-3D in InsightII. The majority of the 100 structures scored high in these evaluations, and we chose the one that ranked highest by PROCHECK. Unlike our previous modeling of the human muscle nAChR (Sine et al., 2002), the linkers between ß-strands 8 and 9 required no further modeling due to a gap of only three residues in the sequence of 7 rather than an insertion of 8–11 residues in muscle nAChR. Two rounds of energy minimization were applied using the program CHARMM (Brooks et al., 1983), version 27b4. The first round constrained the coordinates of all heavy atoms while allowing mobility of all hydrogens. The second round constrained the protein main chain while allowing mobility of all the side chains. In the next stage, the protonation states of all titratable residues and orientation of asparagine and glutamine amide side chains were determined using the WHAT IF program including optimization of the hydrogen bond network (Vriend, 1990). Standard protonation states were found for all residues. The position of the ionizable proton on histidine was N
1 for H114 and N
2 for H62, H104, and H140. The LBD alone numbers 16,835 atoms. Water oxygen atoms were placed close to the LBD using the GRID program (Goodford, 1985) with an energy cutoff of -8 kcal/mol-1 and otherwise default parameters. This added 552 waters. Hydrogens were added to the LBD and the water oxygen atoms so as to optimize the hydrogen bond network using the WHAT IF program (Vriend, 1990). The resulting system was solvated in a 100 Å cubic periodic box using the "tleap" module of AMBER 7 (Pearlman et al., 1995). Finally, 74 Na+ and 54 Cl- ions were randomly swapped with single water molecules, one each, to produce a neutral system and mimic 0.15 M and 0.11 M concentrations for Na+ and Cl- ions. None of the added ions was buried in the LBD. The total number of resulting water molecules was 27,248.
|
). The system was modeled with the AMBER parm99 force field (Wang et al., 2000) and TIP3P water (Jorgensen et al., 1983). Nonbonded interactions were calculated up to a cutoff of 9 Å with particle mesh Ewald (Essmann et al., 1995) for long-range interactions. The whole system was first minimized for 200 steps using the steepest descent algorithm to remove bad contacts. A short 5 ps molecular dynamics simulation at 50 K and constant volume equilibrated the solvent with the LBD fixed. The whole system was then gradually warmed up to 300 K in six simulations of 5 ps at 50 K intervals at constant pressure (1 atm) then extended to 10 ns at 300 K. All dynamics simulations were run with SHAKE (Ryckaert et al., 1977) and a 2-fs time step. A Berendsen thermostat and barostat with isotropic scaling (Berendsen et al., 1984) maintained temperature and pressure. The coupling time was 1 ps for the thermostat, 0.1 ps for the barostat during equilibration, and 1 ps for the barostat afterward. The equilibration phase was run on a dual processor Alpha EV67 and the remainder on Blue Horizon at the San Diego Supercomputer Center.
ACh docking
ACh was docked to each binding site in snapshots selected from the MD simulation to assess how intact the binding site remained during the simulation. ACh was optimized in the extended conformation using Gaussian 98 (HF/6-31G*) and parameterized with Gasteiger-Marsili charges (Gasteiger and Marsili, 1980), which use the united atom representation for nonpolar hydrogens. With all four of its dihedrals flexible, ACh was docked to each binding site in snapshots taken every 10 ps using the AutoDock 3.0.5 program (Morris et al., 1998). A standard docking protocol with default parameters and the Lamarckian genetic algorithm was used, with the following exceptions. To achieve a more thorough search, the initial population size was set to 200, and 50 independent runs were performed per binding site on a 30 Å cubic grid centered on the atom C189:S
with grid spacing 0.375 Å. The positions of the nitrogen and carbonyl carbon atoms of ACh from the 50,000 docked positions at each binding site were placed on a 30 Å cubic grid of spacing 1 Å. This relatively large grid makes docking possible outside the binding site.
|
RESULTS |
|---|
|
TOPABSTRACTINTRODUCTIONMETHODRESULTSDISCUSSIONCONCLUSIONACKNOWLEDGEMENTSREFERENCES |
|---|
Homology model
Homology models for 7 have already been published by Schapira et al. (2002) and Le Novère et al. (2002) for the very similar chick instead of the human. The alignment of Schapira and co-workers is quite different from ours regarding the position of gaps and insertions. Given that our alignment is backed up by the lysine scanning mutagenesis technique (Sine et al., 2002), we have greater confidence in the alignment used here. Regarding the other published homology model, the alignment and, consequently, backbone of our human 7 model is almost identical with the chick 7 model of Le Novère and co-workers. A trivial difference is the residue numbering. The numbering differs by one because Le Novère et al. omit E1, present in our model, such that F2 in our numbering is F1 in theirs. The more significant difference is in the alignment from residues E1 to P15, which is shifted one residue lower in their structure such that their residue F2 matches our residue E1 up to their K11 matching our V12. The alignment by Le Novère and co-workers gains sequence identity at L6 and K7 at the cost of an extra insertion in 7. Which alignment is better is uncertain, given the low sequence identity in this region. The backbones superimpose identically except with small deviations at loops 12-15, 24-26, 72-74, 94, 138-139, 162-165, and 169-170. There is quite a bit of variation in the side-chain placing throughout the model, though the binding-site side chains are very similar. The nearest side chains to the binding site that differ substantially in position are those of N52 and Y167.
Energetic stability
The LBD's potential energy and potential energy components, smoothed over 100-ps intervals, are presented in Fig. 2 for the duration of the simulation. The potential energy during the 30-ps warming-up stage has been omitted. The behavior of the potential energy is characteristic of having a homology model as the starting structure. After an initial large relaxation, the potential energy continues to decrease more sporadically. Excluding the initial relaxation in the minimization and the warming-up stage, there is a total decrease of 600 kcal/mol-1 in the LBD. This change translates to 0.06 kBT per atom, still a small fraction of the thermal energy, 3/2 kBT. Approximately half of the relaxation comes from the dihedral term, a third from the nonbonded term, and a sixth from the bond/angle term. Unlike the other two terms, the nonbonded term fluctuates quite substantially, presumably because this term, unlike the others, contains interactions with the solvent. If the solvent is included, the energy of the whole system including solvent (not shown) appears to have leveled off after 8 ns. Overall, 10 ns of simulation leads to a reasonable improvement in energy. It is unlikely that the final structure represents the activatable state after this relatively short time. A much longer simulation of probably impractical length at this stage would be required to fully relax the structure to that point.
|
s in each subunit and the average relative to the homology model, smoothed over 500-ps intervals to improve readability. The RMSD is calculated after superimposing each frame of the whole LBD on the starting homology model. The colored lines are the RMSD for each subunit and the black line is the average over all five subunits. The RMSD increases rapidly at first due to thermal vibrations and relief of unoptimized interactions in the homology model. At 
2.5 Å, the RMSD levels off a moderate amount. All subunits behave similarly with the exception of SB, which has a larger RMSD, the reason for which is apparent in the next section. Although it is not known at this stage how close the homology model is to the real 7 structure, this relatively stable behavior of RMSD further supports the homology model as a plausible structure.
|
, also from each frame every 1 ps, is illustrated in Fig. 3 B. Mapping the RMSF onto the subunit structure in Fig. 4 and coloring according to value shows where the mobility lies. Blue is least mobile, green intermediate, and red the most mobile. Overall, the structure is very stable, particularly the ß-strand regions colored blue with RMSF <1 Å. This stability is encouraging for the homology model. Excluding the chain termini, mobile parts are the helix 3-13 and the loops 23-25, 43-46, 62-74, 130-139 (Cys loop), 152-153 (loop B), 160-166 (loop G), and 186-193 (loop C). This pattern of RMSF for all five subunits is reasonably consistent among themselves and with the experimental B-factors of AChBP (Brejc et al., 2001) plotted in Fig. 3 C. Here, the MSF (mean-square fluctuation) rather than the RMSF is plotted for consistency with the B-factors. B-factors, shown in gray and scaled by 100, are plotted using the alignment in Fig. 1. The simulation of the 7 model displays flexibility in the expected places and the secondary structure remains constant and intact for the duration of the simulation. Of the seven loops that contribute to the ligand binding site, loops B, C, and (to a lesser extent) G are more flexible, whereas A, D, and F are fairly rigid. The most flexible loop is the Cys loop at the base of the LBD that most probably interacts with the transmembrane region between the M2 and M3 helices. In SB, this loop moves by as much as 7 Å toward the right as drawn in Fig. 5. The new loop position is stabilized by a salt bridge between R132 on the Cys loop and E44 on the same subunit. This motion contributes substantially toward the increase in RMSD for SB (Fig. 3 A). However, this same salt bridge still forms for three other subunits, SA, SD, and SE, although with less contact. The implications of this will be discussed later.
|
|
 | 1 |
atoms, mi is the mass of the atom, and 
is the distance of the atom from the central axis of the receptor (the z axis). Fig. 3 D shows how Rg of the whole LBD changes (in black) during the simulation, smoothed over 500-ps intervals. After an initial sharp drop from 25.5 Å, the LBD undergoes a gradual expansion by 0.5 Å, followed by a slight contraction toward the end. This slight expansion may arise from a number of causes. It may be relief from the confinement of crystal packing; it may result from suboptimal packing in the homology model; or it may represent natural motion of the LBD, showing that it is capable of breathing.
This breathing motion may be more precisely interpreted by breaking it down into contributions from each subunit. Fig. 3 D also shows Rg for each subunit about the same LBD axis, smoothed by 500 ps. Two subunits, SB and SE, are responsible for the initial drop in total Rg seen earlier. All five subunits then expand outward slowly before contracting slightly toward the end. After 10 ns, relative to the homology model, two subunits, SB and SE, have shifted inward; SA and SD have moved outward, whereas SC ends up close to the average. These motions lead to a slightly distorted, asymmetric structure. The sizes of these motions are small, ranging from 0.6 Å for SA, SB, and SC to 1.0 Å for SD and SE. These motions are depicted in Fig. 6, which shows the structures of the model and after 10 ns. For the model structure, the large, straight arrows indicate the motion that the center of mass (C atoms) of each subunit undergoes to move to the 10-ns structure. Likewise, the arrows on the 10-ns structure indicate how the subunits would move to return to the model structure and, correspondingly, are the reverse of those shown for the model structure. The smaller arrows represent the same motion at a finer level averaged over triplets of consecutive Cs. The large arrows are consistent with the values of Rg.
|
These different subunit motions would be expected to have an effect on their intersubunit interfaces and, consequently, the binding sites. A number of different types of interface should arise. IAB and IDE both involve a shear motion of the plus subunit moving out relative to the minus subunit and a slight separation due to the outward movement of SD and SA. IBC remains reasonably stable between two subunits moving inward together. ICD and IEA undergo compression due to the slight radially inward motion of the plus subunit. They also undergo a shear motion this time, with the plus subunit moving inward relative to the minus side. How these expected movements are borne out in the binding sites is examined in the next section.
Binding site changes
A range of different behavior is found for each of the five binding sites. Fig. 7 depicts the binding site motions at all five interfaces. Each row represents one binding site as viewed from the side with the plus side on the left and the minus side on the right. The identity of the participating subunits is indicated by the colored letters. On each row, the two structures on the left are the binding site residues of the model and after 10 ns. Shown are residues Y92, W148, Y194, and C190 on the plus side, and W54 and L118 on the minus side. Two other important residues making up the binding site, Y187 and C189, both on the plus side, are omitted for clarity and behave similarly to C190. The arrows, here marked in green, point in the direction that the C atoms of these residues move to go from the current structure to the other structure. The third and fourth structures in each row show the result of ACh docking, to be discussed later.
|
atoms is generally consistent with the global subunit motion, supporting the premise that the subunit motion occurs in real nAChRs rather than being simply the artifact of an internal localized motion. At IAB and IDE, residues on the plus side move out whereas those on the minus side move in—the same shear motion seen earlier. All residues at IBC move inward as do SB and SC. Some compression is seen at ICD and IEA, although the shear motion with the plus side moving inward relative to the minus appears to be absent. The interfaces appear unable to shear in this direction. At the ICD interface, all the residues move outward. Even though SC moves in slightly, it also rotates clockwise, explaining the outward motion at this interface. The main motion at IEA is a shift of both subunits counterclockwise and some compression of the binding site residues. The clockwise rotation of SA offsets its outward motion at IEA, leading to little net radial motion at this interface.
The time series of these motions may be seen in a selection of distance plots in Fig. 3, F–H, smoothed by 500 ps. The color refers to the plus subunit at the interface. The following distances are defined between the C atoms of the following residue pairs. The W148–L118 distance captures the width of the binding site between the two subunits. Starting at 9 Å, this distance decreases for four of the five interfaces, especially at ICD. The drop at the start is quite rapid. IDE is the exception where shearing and slight separation take place. How the atoms move to alter these distances may be seen from the arrows on these residues in Fig. 7. A particularly conspicuous motion in the binding site is the inward movement of loop C containing Y187, C189, and C190 into the binding site. The C distance, W148–C189, plotted in Fig. 3 G, decreases by 2–4 Å at IAB, IBC, and IDE, and less so at ICD and IEA. Residues Y187, C189, and C190 move fairly concertedly, and all their distances with W148 behave similarly for each interface. The decreased contraction of loop C at ICD and IEA may be explained by the shear motion described earlier whereby the minus subunits, SD and SA, move outward (Fig. 6) to counteract loop C. The fluctuation of this distance is quite large and is due to the mobility of loop C, also evident in its greater flexibility in the RMSF (Figs. 3 B and 4) and the larger arrows of C190 in Fig. 7. The backbone of W148, the other atom defining this distance, is much more stable. One final distance that differs significantly among subunits is Y194–W54, which loosely measures the height of the binding site. This increases at IDE due to the shearing and separation, and at ICD for which Y194 moves upward and marginally at IAB, mentioned earlier, yet contracts at IBC and IEA. As well as distances that change, some other clusters of residues remain more together and are separated by fairly constant distances (not shown). One unit is Y92, W148, and Y194, with the exception of ICD in which Y194 moves upward away from the other two. Another such unit is W54 and L118, which are also separated by a fairly constant distance. A third unit is Y187, C189, and C190. Interestingly, Y92 also remains relatively fixed with respect to W54 and L118 across the interface, whereas W148 does not. Y92 also remains fairly fixed with respect to Y187, except early on in SC. Hence, Y92 appears to act as an intermediary between W148 and W54. This behavior is consistent with the proximity of each residue's side chain.
The side chains of the binding site residues display a few notable differences after 10 ns. The most striking is the crowding of the binding site. Starting from a relatively open site, the side chains move in to partially fill the site. Three residues possess side chains that differ markedly between interfaces: W54, W148, and L118. The chief residue responsible for the shrinking is W148. In all cases but IBC, the lower half of the W148 indole ring swings into the site. A small dihedral change is sufficient to bring this about and does not involve a complete transition between rotameric states, as both side-chain dihedrals, 
1 and 2, remain gauche+ and gauche-, respectively. A larger change is seen at ICD, for which 1 of W148 goes from gauche+ to trans, twisting the indole ring to the side in response to the collapsing binding site. At the other extreme, the W148 indole ring at IBC remains fixed in place because W54 swings up and blocks it from doing so. This other tryptophan, W54, also undergoes quite different behavior at each interface. Most of the difference lies in 2. At IAB, ICD, and IDE, 2 remains gauche- and twists down the indole ring slightly, opening up the entrance to the binding site. However, at IBC and IEA, 2 swings up the ring to become gauche+, further crowding the binding site. In the case of IEA, W54 even manages to get behind W148, further raising that residue. Coincidentally, both of these interfaces are ones that are compressed. IDE is the only interface where 1 of W54 also changes, going from trans to gauche-, pulling the indole ring farther from the site. L118 is the third residue whose side-chain dihedrals differ among interfaces. One notable feature about this residue was that during the first stages of minimization, its side chain in each subunit adopted different conformations. In particular, 1 is trans at IAB and IDE, and gauche+ for the other three interfaces. During the simulation, both 1 and 2 largely remain gauche- in IBC and IEA, both trans in IAB and ICD, and respectively gauche- and trans in IDE. The trans position of 1 keeps the side chain more out of the binding site, whereas the gauche conformations point it in. Most of the other side chains in the binding site remain fairly constant in conformation. Their only major fluctuations are for Y92 in IAB between trans and gauche- and C189 and C190 in IBC. The occasional tyrosine ring flips also occur.
ACh docking
The results of docking ACh to each binding site are illustrated in the third and fourth columns of Fig. 7. The residue positions are those of the homology model and after 10 ns. Green regions show where the ACh's carbonyl carbon docks, and blue regions where the nitrogen docks. More specifically, each colored region is an isocontour with one docking hit per grid cube (1 Å3) for docking to the starting model, and 250 hits per grid cube for docking to the snapshots. With 50 docking runs per snapshot and 1000 snapshots, 50,000 docking conformations were generated per binding site. What is abundantly clear is that, whereas it docks well to the homology model, ACh docks poorly to the binding site for most of the structures generated by the simulation. For the starting structure, the docked positions are very similar to the docking of ACh to 7 by Le Novère et al. (2002). This is the site in which the HEPES buffer molecule lies in the AChBP crystal structure (Brejc et al., 2001). The choline end of ACh lies at the middle of the W148 indole ring nestled between the other aromatic side chains of W54 and Y92. The acetyl end of ACh lies near the top of the W148 indole ring and is bounded by N106, V107, L108, and Q116 (these residues omitted from the figures for clarity). The choline group makes cation-
interactions with the aromatic groups; the acetate methyl group lies in a hydrophobic cleft, whereas the carbonyl oxygen forms a hydrogen bond either with the backbone amide hydrogen of L118 or the polar H1 of W148. In the case of some interfaces, particularly IAB and IEA, the orientation of ACh flips by 180°. This probably arises either from a force-field deficiency, particularly in the cation- interaction, or is characteristic of the homology model.
The docking positions for the homology model are not seen in the docking results for the whole simulation since they fail to occur often enough to surpass the 250-hit threshold. Rather, docking positions predominate elsewhere. The most common cluster is on the face of W54 as it swings away from the binding site, present in all cases except IBC. Apart from a smattering of clusters on the LBD surface, two other interesting binding sites are observed. One lies in IDE below the W148 indole ring. At this interface occurs the shearing motion and slight separation of the two subunits. The outward motion of W54 and the extra space due to the interface separation contribute to increasing the space that ACh occupies. The other site lies in IBC behind W148 instead of in front. This is the interface for which the W148 indole ring remains relatively stationary and SB and SC move inward in concert. This site opens up due to the motion of the ring of F103 on the minus interface early on in the simulation that blocks this site in the other interfaces. F103 itself backs onto an indentation in the inner pore. Thus a total of four binding sites are observed including the initial one. Assuming that the site preferred in the model is the actual ACh binding site in the ligand-bound state, the side chains that move during the simulation to overlap with ACh in this position would be those that have to move aside for ACh to bind and activate the receptor. The principal residues that sterically overlap with ACh in its preferred docked position in the model are W148 and L118 (the most followed by the three tyrosines); Y194, Y187, Y92; and C189, and C190 on loop C, with W54 only clashing at IBC and less so at IEA.
|
DISCUSSION |
|---|
|
TOPABSTRACTINTRODUCTIONMETHODRESULTSDISCUSSIONCONCLUSIONACKNOWLEDGEMENTSREFERENCES |
|---|
7 based on the AChBP crystal structure to test the validity of the structure. Ideally, the simulation would be started from a crystal structure of 7 but given the current absence of knowledge of such a structure, a homology model is the next best alternative if any modeling work at the molecular level is to be attempted. The simulation leads to a large reduction in the energy of the homology model, producing an improved, more relaxed, structure. The overall structure of the LBD remains fairly constant during this relaxation at the secondary and tertiary levels of structure, supporting the validity of the initial model. The RMSF only increases to a moderate value, and what backbone flexibility is seen in the crystal structure B-factors, predominantly loop motion, is also seen in the simulation MSF. Much of the remaining relaxation is due to small atom rearrangements and side-chain repacking, an aspect not examined in this work except in the binding site. The accuracy of the model is expected to be greater in the binding site given the greater sequence identity between 7 and AChBP. Based on the continually decreasing trend in energy, 10 ns appears insufficient to fully relax the model, yet it is still adequate to assess the model quality and to glean some information on the flexibility of the LBD.
Studying the LBD flexibility is the other purpose of this investigation. 7 has pentameric symmetry, and over a long enough time, each subunit might be expected to experience the same types of motion. However, because molecular motion is stochastic and the simulation time is short, each subunit will undergo different motions during the simulation. Furthermore, by observing the types of motion occurring in the LBD and exploiting the different patterns of motion seen in different subunits, insight is gained into how motions in the binding site in response to ligand binding might correlate with other motions elsewhere that connect to the ion channel in the transmembrane region. The analysis of the simulation reveals three types of motion taking place. These are motions of the side chains in the binding site, loops, and whole subunits. Motion can arise from a number of sources. It may either be a consequence of going from the crystalline to the liquid phase, relaxation of the homology model, relaxation to the closed state, the force field, absence of the transmembrane region and lipid bilayer, or random noise due to flexibility in the structure. It is difficult at this stage, given the homology model and the short timescale, to be certain of the nature of each type of motion.
The foremost motions relevant to ligand binding are those taking place in the binding site. The main result observed here in the simulation is the shrinking of the binding site dimensions. The principal contributions to this process are the swinging in of the W148 indole ring, the inward motion of loop C, and the decrease in separation between the backbones of the Y92/Y194/W148 and W54/L118 units. Less important contributions to shrinking come from the Y92, Y187, Y194, and L118 side chains. The contribution of these three tyrosines-to-agonist binding has been noted in mutagenesis studies (Sine et al., 1994). W54 more commonly swings away from the binding site, but it may move in as well in the cases of IBC and IEA. There is strong evidence for the importance of W148 and W54 to ligand binding (Corringer et al., 1995; Zhong et al., 1998), and the different conformations of these residues among binding sites observed in the simulation suggests that they have some flexibility. The second contribution to binding site shrinkage is the inward motion of loop C.
The absence of a ligand creates space, allowing the loop to move in. The final contribution, the distance between residue backbones at the interface, appears to be correlated with subunit motion, discussed below. The trend here is not so uniform. IBC, ICD, and IEA display slight compression whereas IAB and IDE are more open. The sum effect of these three contributions still leads to shrinking of the binding site for all five cases. This shrinking all but prevents ACh from docking as would be expected in the crystal structure. Consequently, ACh prefers to dock elsewhere, either over W54, or under or behind W148. This docking site behind W148 is rather speculative. It arises due to a displacement of the F105 side chain, forming a cavity there. ACh may reach this position if it is able to burrow under W148 as it does when docking to IDE. Overall, when ACh is present in the simulation, it would presumably encourage expanding motions of the binding site to enter the site. Although the side chains and loop C appear to encroach on the binding site, their flexibility suggests that they would not act as too large an obstacle for the entrance of ACh. Their displacement upon ACh binding may even contribute to the resulting conformational change. On the other hand, the motions may be a result of imperfect side-chain packing in the homology model or inadequate solvation of the binding pocket in the initial setup. The use of GRID to place waters and some solvent-only equilibration lead to partially solvated binding sites, although the actual extent of hydration of the largely hydrophobic binding sites remains unknown. If these binding site motions are realistic, then it adds further weight to the idea that the homology model of 7 and AChBP are in the ligand-bound state, for ACh binding would only be favored when the binding site is open.
It is not known by what mechanism ligand binding induces the opening of the transmembrane pore. Based on the flexibility observed in the simulation, one possible means is via subunit motions. It should be noted that modeling the subunit motions as translations and rotations is only approximate. Subunits also have more complex and subtle internal motions that indirectly contribute to translational and rotational motion, although no clearly characterizable flexibility at an intermediate level such as within the ß-strands was identified here. Nevertheless, this analysis still serves as a useful tool to dissect motion. The whole pore undergoes an asymmetric breathing during the 10 ns whereby SA and SD move outward and SB and SE inward. Accompanying these translations are small rotations. SA and SC twist clockwise and SB and SE counterclockwise. Such subunit motion could plausibly be transmitted directly to the transmembrane helices forming the ion channel. The sizes of these motions are no larger than 1 Å, but an effect this size, when transmitted to the transmembrane region, may be sufficient to alter ion flow in the narrow pore. Similar rotations for the inner ß-sheets of muscle -subunits have been observed in electron microscopy structures (Unwin et al., 2002). In this case, the sheets are also seen to rotate counterclockwise going from the activated to the activatable state. To test if this type of sheet motion was detectable in the simulations, a similar analysis was performed on partitions of the subunits into inner and outer ß-strands. However, the translational and rotational motions of the inner and outer strands were found to behave similarly to their respective subunits, implying that the subunits remain fairly rigid during the simulation.
The other type of motion observed here that would connect with the transmembrane region is the change in the Cys loop seen in SB. This subunit happens to be one that moves inward and counterclockwise. Such motions may encourage the displacement of the loop in this direction. The formation of a salt bridge between R132 on the Cys loop and E44 would also stabilize this motion. Salt bridges involving the Cys loop are particularly noteworthy because of their proximity to the transmembrane region. In the GABAA -subunit, a variable pattern in the formation of salt bridges between K279 in the M2-M3 linker region and D57 or D149 has been observed depending on GABA concentration (Kash et al., 2003). It was found that the K279-D149 salt bridge only forms in the presence of high GABA concentration, whereas the K279-D57 bridge forms independently of GABA. A slightly different salt bridge pattern may arise for 7 with R132 choosing between E44 and possibly D265 in the M2-M3 linker region. If the R132-E44 salt bridge, characteristic of the activatable state, is preferred by SB, then this is further evidence that SB is moving toward the activatable state.
Some interesting relationships may be seen between the binding site motions and the subunit motions. The relative motion of the subunits is matched by backbone motions of residues in three of the binding sites, IAB, IBC, and IDE. IAB and IDE experience an outward shear of the plus subunit. At IBC both subunits move inward together. However, the binding-site residues at ICD and IEA, expected to have an outward shear of the minus subunit, compress slightly but largely move together, outward in the case of ICD and counterclockwise for IEA. Another characteristic of these two binding sites is that loop C moves inward much less because the minus subunit, with which loop C has contact, moves outward. Overall, backbone motion of the binding site residues parallels subunit motion to a large extent. More intriguing are the differences that arise at each interface. There is an interesting correlation between the W54 and L118 side chains and the subunit motions. At IBC and IEA, for which 1 of L118 begins as gauche+, pointing more into the binding site, W54 also adopts a gauche+ orientation and also swings into the binding site. At the same time, SB and SE move inward and rotate counterclockwise. However, at IAB and IDE, for which 1 of L118 begins as trans, then W54 moves into the gauche- conformation, swinging away from the binding site. This time, SA and SD move outward, SA rotates clockwise whereas SD experiences no net rotation. ICD is less clear-cut but resembles IAB and IDE. How these subunit and binding-site motions are connected remains subtle. It appears that when these two side chains point toward the binding site, there is a greater tendency for the interface to be tighter and the plus subunit to be drawn in. This pattern appears to repeat further down the interface at the Cys loop, which moves toward the minus interface at IBC.
The asymmetry observed and its relationship with binding site structure have implications for ACh binding and the evolution of receptors themselves. The strength of binding would be expected to be different at each type of interface, some better and some worse. The more open interfaces, IAB and IDE, with W54 more favorably placed, would be expected to bind ACh more readily. Once they do bind, the interactions of ACh with the binding-site side chains may induce the plus subunits to drift back inward to a more symmetrical structure resembling AChBP. This would then improve the binding at the other interfaces, a possible mechanism of cooperativity. It may be that the binding of ACh in these two sites alone is sufficient for nAChR activation. That only two sites might be favorable for binding and, upon binding, necessary for activation, would be consistent with Hill coefficients for receptor activation by agonist no >2 observed for 7 (Corringer et al., 1995; Gopalakrishnan et al., 1995). On a possible relationship between the observed asymmetry and evolution, the asymmetry is reminiscent of inherently asymmetrical nAChRs such as the neuronal 4ß2 and muscle 2ß with respective subunit sequences ßßß and ß. These receptors contain two ACh-binding -subunits separated by one non-ACh-binding subunit. In this simulation, SA and SD might play the role of the 4ß2 and muscle -subunits based on the outward displacements they undergo in the relaxation compared with the electron microscopy structures of Unwin et al. (2002) and the supposedly more favorable interfaces IAB and IDE for ACh binding. It might also be argued that SB and SE resemble the -subunits based on their counterclockwise rotations observed in the simulation and rotations in the same direction detected in the electron microscopy structures. Whereas this asymmetry may occur spontaneously for 7, it appears to be enforced for the 4ß2 and muscle nAChRs by the different subunit compositions. Mutations in nAChR subunits away from -subunits may have been tolerated as long as they preserved this capability of asymmetry.
|
CONCLUSION |
|---|
|
TOPABSTRACTINTRODUCTIONMETHODRESULTSDISCUSSIONCONCLUSIONACKNOWLEDGEMENTSREFERENCES |
|---|
7 LBD has demonstrated its stability. In addition, it has also revealed some of the motions that take place in the LBD and suggested how ligand binding correlates with larger-scale subunit motions that would connect with the transmembrane region. In particular, conformational changes in the side chains of W148, W54, and L118 have been correlated with larger subunit motions in the LBD and motion of the Cys loop in SB. The resulting structure may more closely resemble the activatable state, and the asymmetry with binding sites of different affinity for ACh may be a natural property of 7. These interpretations are speculative and require closer examination, as the short 10-ns simulation time would have sampled only a portion of the repertoire of motions available to the LBD, and the simultaneity of motions may be coincidental. There is also an element of uncertainty as to how representative the motions are, given the use of a homology model. Having a full crystal structure and better simulation sampling methods would provide the ideal solution to addressing these deficiencies. In the meantime, we plan to undertake further simulations of the model with a range of bound ligands including agonists, antagonists, and potentiators, and with the lipid bilayer and most of the transmembrane domain present, now made possible by the recent electron microscopy transmembrane structure minus the vestibule region between the M3 and M4 helices (Miyazawa et al., 2003). These simulations should make clearer the causes and possible functional relevance of the different motions observed. |
ACKNOWLEDGEMENTS |
|---|
|
TOPABSTRACTINTRODUCTIONMETHODRESULTSDISCUSSIONCONCLUSIONACKNOWLEDGEMENTSREFERENCES |
|---|
This work has been supported by grants from the National Science Foundation; the National Institutes of Health; the San Diego Supercomputer Center; National Biomedical Computation Resource at the University of California, San Diego; the Center for Theoretical Biological Physics; Accelrys; and the W. M. Keck Foundation.
Submitted on May 15, 2003; accepted for publication August 4, 2003.
|
REFERENCES |
|---|
|
TOPABSTRACTINTRODUCTIONMETHODRESULTSDISCUSSIONCONCLUSIONACKNOWLEDGEMENTSREFERENCES |
|---|
Arias, H. R. 1997. Topology of ligand binding sites on the nicotinic acetylcholine receptor. Brain Res. Brain Res. Rev. 25:133–191.[Medline]
Berendsen, H. J. C., J. P. M. Postma, W. F. van Gunsteren, A. Di Nola, and J. R. Haak. 1984. Molecular dynamics with coupling to an external bath. J. Chem. Phys. 81:3684–3690.
Brejc, K., W. J. van Dijk, R. V. Klaassen, M. Schuurmans, J. van der Oost, A. B. Smit, and T. K. Sixma. 2001. Crystal structure of an ACh-binding protein reveals the ligand-binding domain of nicotinic receptors. Nature. 411:269–276.[Medline]
Brooks, B. R., R. E. Bruccoleri, B. D. Olafson, D. J. States, S. Swaminathan, and M. Karplus. 1983. CHARMM: A program for macromolecular energy minimization and dynamics calculations. J. Comp. Chem. 4:187–217.
Changeux, J. P., and S. J. Edelstein. 1998. Allosteric receptors after 30 years. Neuron. 21:959–980.[Medline]
Changeux, J. P., and S. J. Edelstein. 2001. Allosteric mechanisms in normal and pathological nicotinic acetylcholine receptors. Curr. Opin. Neurobiol. 11:369–377.[Medline]
Corringer, P. J., J. L. Galzi, J. L. Eiselé, S. Bertrand, J. P. Changeux, and D. Bertrand. 1995. Identification of a new component of the agonist binding site of the nicotinic 7 homooligomeric receptor. J. Biol. Chem. 270:11749–11752.
Corringer, P. J., N. Le Novère, and J. P. Changeux. 2000. Nicotinic receptors at the amino acid level. Annu. Rev. Pharmacol. Toxicol. 40:431–458.[Medline]
Essmann, U., L. Perera, M. L. Berkowitz, T. Darden, H. Lee, and L. G. Pedersen. 1995. A smooth particle mesh Ewald method. J. Chem. Phys. 103:8577–8593.
Gasteiger, J., and M. Marsili. 1980. Iterative partial equalization of orbital electronegativity: a rapid access to atomic charges. Tetrahedron. 36:3219–3228.
Gillilan, R. E., and F. Wood. 1995. Visualization, virtual reality, and animation within the data flow model of computing. Comput. Graphics. 29:55–58.
Goodford, P. J. 1985. A computational procedure for determining energetically favorable binding sites on biologically important macromolecules. J. Med. Chem. 28:849–857.[Medline]
Gopalakrishnan, M., B. Buisson, E. Touma, T. Giordano, J. E. Campbell, I. C. Hu, D. Donnely-Roberts, S. P. Arneric, D. Bertrand, and J. P. Sullivan. 1995. Stable expression and pharmacological properties of the human 7 nicotinic acetylcholine-receptor. Eur. J. Pharmacol. 290:237–246.[Medline]
Grutter, T., and J. P. Changeux. 2001. Nicotinic receptors in wonderland. Trends Biochem. Sci. 26:459–463.[Medline]
Hogg, R. C., M. Raggenbass, and D. Bertrand. 2003. Nicotinic acetylcholine receptors: from structure to brain function. Rev. Physiol. Biochem. Pharmacol. 147:1–46.[Medline]
Humphrey, W., A. Dalke, and K. Schulten. 1996. VMD: visual molecular dynamics. J. Mol. Graph. 14:33–38.[Medline]
Jorgensen, W. L., J. Chandrasekhar, J. D. Madura, R. W. Impey, and M. L. Klein. 1983. Comparison of simple potential functions for simulating liquid water. J. Chem. Phys. 79:926–935.
Karlin, A. 2002. Emerging structure of the nicotinic acetylcholine receptors. Nat. Rev. Neurosci. 3:102–114.[Medline]
Kash, T. L., A. Jenkins, J. C. Kelley, J. R. Trudell, and N. L. Harrison. 2003. Coupling of agonist binding to channel gating in the GABAA receptor. Nature. 421:272–275.[Medline]
Laskowski, R. A., M. W. MacArthur, D. S. Moss, and J. M. Thornton. 1993. PROCHECK: a program to check the stereochemical quality of protein structures. J. Appl. Crystallogr. 26:283–291.
Le Novère, N., P. J. Corringer, and J. P. Changeux. 1999. Improved secondary structure predictions for a nicotinic receptor subunit: incorporation of solvent accessibility and experimental data into a two-dimensional representation. Biophys. J. 76:2329–2345.
Le Novère, N., T. Grutter, and J. P. Changeux. 2002. Models of the extracellular domain of the nicotinic receptors and of agonist- and Ca2+-binding sites. Proc. Natl. Acad. Sci. USA. 99:3210–3215.
Merritt, E. A., and D. J. Bacon. 1997. Raster3D: photorealistic molecular graphics. Method. Enzymol. 277:505–524.[Medline]
Miyazawa, A., Y. Fujiyoshi, M. Stowell, and N. Unwin. 1999. Nicotinic acetylcholine receptor at 4.6 Å resolution: transverse tunnels in the channel wall. J. Mol. Biol. 288:765–786.[Medline]
Miyazawa, A., Y. Fujiyoshi, and N. Unwin. 2003. Structure and gating mechanism of the acetylcholine receptor pore. Nature. 423:949–955.[Medline]
Molles, B. E., I. Tsigelny, P. D. Nguyen, S. X. Gao, S. M. Sine, and P. Taylor. 2002. Residues in the subunit of the nicotinic acetylcholine receptor interact to confer selectivity of Waglerin-1 for the - subunit interface site. Biochemistry. 41:7895–7906.[Medline]
Montal, M., and S. J. Opella. 2002. The structure of the M2 channel-lining segment from the nicotinic acetylcholine receptor. Biochim. Biophys. Acta. 1565:287–293.[Medline]
Morreale, A., F. Maseras, I. Iriepa, and E. Gálvez. 2002. Ligand-receptor interaction at the neural nicotinic acetylcholine binding site: a theoretical model. J. Mol. Graph. 21:111–118.
Morris, G. M., D. S. Goodsell, R. S. Halliday, R. Huey, W. E. Hart, R. K. Belew, and A. J. Olson. 1998. Automated docking using a Lamarckian genetic algorithm and an empirical binding free energy function. J. Comput. Chem. 19:1639–1662.
Ortells, M. O. 1997. Prediction of the secondary structure of the nicotinic acetylcholine receptor nontransmembrane regions. Proteins. 29:391–398.[Medline]
Ortells, M. O., G. E. Barrantes, C. Wood, G. G. Lunt, and F. J. Barrantes. 1997. Molecular modelling of the nicotinic acetylcholine receptor transmembrane region in the open state. Protein Eng. 10:511–517.
Pearlman, D. A., D. A. Case, J. W. Caldwell, W. S. Ross, T. E. Cheatham III, S. DeBolt, D. Ferguson, G. Seibel, and P. Kollman. 1995. AMBER, a package of computer programs for applying molecular mechanics, normal mode analysis, molecular dynamics and free energy calculations to simulate the structural and energetic properties of molecules. Comput. Phys. Commun. 91:1–41.
Ryckaert, J.-P., G. Ciccotti, and H. J. C. Berendsen. 1977. Numerical integration of the Cartesian equations of motion of a system with constraints: molecular dynamics of n-alkanes. J. Comput. Phys. 23:327–341.
ali, A., and T. L. Blundell. 1993. Comparative protein modeling by satisfaction of spatial restraints. J. Mol. Biol. 234:779–815.[Medline]
Schapira, M., R. Abagyan, and M. Totrov. 2002. Structural model of nicotinic acetylcholine receptor isotypes bound to acetylcholine and nicotine. BMC Struct. Biol. 2:1.[Medline]
Sine, S. M. 2002. The nicotinic receptor ligand binding domain. J. Neurobiol. 53:431–446.[Medline]
Sine, S. M., P. Quiram, F. Papanikolaou, H. J. Kreienkamp, and P. Taylor. 1994. Conserved tyrosines in the subunit of the nicotinic acetylcholine receptor stabilize quaternary ammonium groups of agonists and curariform antagonists. J. Biol. Chem. 269:8808–8816.
Sine, S. M., H. L. Wang, and N. Bren. 2002. Lysine scanning mutagenesis delineates structural model of the nicotinic receptor ligand binding domain. J. Biol. Chem. 277:29210–29223.
Smit, A. B., N. I. Syed, D. Schaap, J. van Minnen, J. Klumperman, K. S. Kits, H. Lodder, R. C. van der Schors, R. van Elk, B. Sorgedrager, K. Brejc, T. K. Sixma, and W. P. M. Geraerts. 2001. A glia-derived acetylcholine-binding protein that modulates synaptic transmission. Nature. 411:261–268.[Medline]
Tikhonov, D. B., and B. S. Zhorov. 1998. Kinked-helices model of the nicotinic acetylcholine receptor ion channel and its complexes with blockers: simulation by the Monte Carlo minimization method. Biophys. J. 74:242–255.
Unwin, N., A. Miyazawa, J. Li, and Y. Fujiyoshi. 2002. Activation of the nicotinic acetylcholine receptor involves a switch in conformation of the subunits. J. Mol. Biol. 319:1165–1176.[Medline]
Vriend, G. 1990. WHAT IF: a molecular modeling and drug design program. J. Mol. Graph. 8:52–56.[Medline]
Wang, J. M., P. Cieplak, and P. A. Kollman. 2000. How well does a restrained electrostatic potential (RESP) model perform in calculating conformational energies of organic and biological molecules? J. Comput. Chem. 21:1049–1074.
Zhong, W. G., J. P. Gallivan, Y. Zhang, L. Li, H. A. Lester, and D. A. Dougherty. 1998. From ab initio quantum mechanics to molecular neurobiology: a cation- binding site in the nicotinic receptor. Proc. Natl. Acad. Sci. USA. 95:12088–12093.
This article has been cited by other articles:
 |
|
 |
|
  M. Moroni, R. Vijayan, A. Carbone, R. Zwart, P. C. Biggin, and I. Bermudez Non-Agonist-Binding Subunit Interfaces Confer Distinct Functional Signatures to the Alternate Stoichiometries of the {alpha}4{bet |
