Molecular Dynamics Simulations of the Ligand-binding Domain of an N-Methyl-d-aspartate Receptor*

The mechanism of partial agonism at N-methyl-d-aspartate receptors is an unresolved issue, especially with respect to the role of protein dynamics. We have performed multiple molecular dynamics simulations (7 × 20 ns) to examine the behavior of the ligand-binding core of the NR1 subunit with a series of ligands. Our results show that water plays an important role in stabilizing different conformations of the core and how a closed cleft conformation of the protein might be stabilized in the absence of ligands. In the case of ligand-bound simulations with both full and partial agonists, we observed that ligands within the binding cleft may undergo distinct conformational changes, without grossly influencing the degree of cleft closure within the ligand-binding domain. In agreement with recently published crystallographic data, we also observe similar changes in backbone torsions corresponding to the hinge region between the two lobes for the partial agonist, d-cycloserine. This observation rationalizes the classification of d-cycloserine as a partial agonist and should provide a basis with which to predict partial agonism in this class of receptor by analyzing the behavior of these torsions with other potential ligands.

The mechanism of partial agonism at N-methyl-D-aspartate receptors is an unresolved issue, especially with respect to the role of protein dynamics. We have performed multiple molecular dynamics simulations (7 ؋ 20 ns) to examine the behavior of the ligand-binding core of the NR1 subunit with a series of ligands. Our results show that water plays an important role in stabilizing different conformations of the core and how a closed cleft conformation of the protein might be stabilized in the absence of ligands. In the case of ligand-bound simulations with both full and partial agonists, we observed that ligands within the binding cleft may undergo distinct conformational changes, without grossly influencing the degree of cleft closure within the ligand-binding domain. In agreement with recently published crystallographic data, we also observe similar changes in backbone torsions corresponding to the hinge region between the two lobes for the partial agonist, D-cycloserine. This observation rationalizes the classification of D-cycloserine as a partial agonist and should provide a basis with which to predict partial agonism in this class of receptor by analyzing the behavior of these torsions with other potential ligands.
N-Methyl-D-aspartate (NMDA) 4 receptors are ligand-gated ion channels that play a major role in learning and memory (1). They have also been implicated in a number of injury or disease states including for example, schizophrenia (2)(3)(4)(5). NMDA receptors are a distinct subtype of ionotropic glutamate receptor (6) in that they require both glutamate and glycine for activation and membrane depolarization (7). The receptor is a heterotetrameric cation channel typically comprised of two NR1 subunits that contain the glycine-binding site along with two NR2A-D glutamate-binding sites. The heterogeneity of the complex is further extended if the third subtype (NR3), which also binds glycine, is also taken into account.
The architecture of each subunit of NMDA receptors is similar to that of non-NMDA receptors. Each subunit is comprised of an aminoterminal binding domain that shows homology to the LIVBP (leucineisoleucine-valine-binding protein) (8), three transmembrane helices (M1, M2, and M3), a re-entrant P-loop, and a ligand-binding domain that is comprised of two discontinuous polypeptide chains, S1 and S2, that form two distinct lobes or subdomains, D1 and D2 (Fig. 1A). An isolated D1D2 ligand-binding core has been shown to bind ligands with similar affinity as wild-type receptors (9,10).
Recently, the crystal structure of the ligand-binding domain of the NR1 subunit has been solved by x-ray diffraction with both full (9) and partial (10) agonists. The overall fold and structure (Fig. 1B) is very similar to those observed for crystal structures of the ligand-binding domains for the (non-NMDA) GluR2 (11)(12)(13), GluR5 (14,15), and GluR6 receptors (Refs. 15 and 16; for recent reviews see Refs. [17][18][19]. The major difference is an insertion of 30 residues into loop 1, which also contains two disulfide bridges. Interestingly, the mechanism of partial agonism in NMDA receptors has been proposed to differ from the mechanism in non-NMDA receptors (10). Partial agonists for NMDA receptors appear capable of promoting a degree of domain closure comparable with that observed for full agonists. Furthermore, whether they are full or partial agonists can be discriminated by the conformation of the second strand (residues 750 -755) of the hinge region. Although the crystallographic structures have been a leap forward in our understanding of these receptors, it is important to remember that proteins at physiological temperatures can exhibit important motions that determine water mobility in binding pockets (20) through to large domain motions (21). Our results show not only the importance of receptor dynamics but also how the mechanism of the partial agonist D-cycloserine can be reconciled in light of recent data.

EXPERIMENTAL PROCEDURES
System Preparation-Coordinates for the crystal structures of the NR1 with glycine (Protein Data Bank code 1PB7), D-serine (Protein Data Bank code 1PB8), D-cycloserine (Protein Data Bank code 1PB9), 1-aminocyclopropane-1-carboxylic acid (ACPC; Protein Data Bank code 1Y20), and 5,7-dichlorokynurenic acid (DCKA; Protein Data Bank code 1PBQ) were downloaded from www.rcsb.org (22). A Closed-Apo state was generated by removing the ligand from the glycine-bound structure, and an Open-Apo state was generated by removing the ligand from the DCKA-bound structure. The chemical structures of the ligands are shown in Fig. 1C. Table 1 summarizes the simulations.
The amino and carboxyl termini were acetylated and amidated, respectively, to mimic the continuation of the peptide chain. Missing side chains were built in using the WHATIF "complete a structure" web service (swift.cmbi.kun.nl/WIWWWI/). Because the initial structures did not have loop 1 resolved in the crystal structure, we used Modeler 7v7 (23) to generate the missing residues. Since beginning this research a crystal structure with an intact loop 1 has been solved (Protein Data Bank code 1Y20); comparison of the modeled loop is very similar to this crystal conformation for the closed cleft simulations (RMSD Յ 1.15 Å) but not for the open cleft ones (RMSD ϭ 4.3 Å). This can be explained by the resolved crystal loop structure being from a closed ligand-binding domain. Presumably this loop was not originally observed because it has high mobility compared with the rest of the structure. Indeed, the B-factors for loop 1 from the 1Y20 structure are generally higher than for the rest of the protein. The ligands D-cycloserine, ACPC, and DCKA were parameterized using PRODRG (24).
Simulation Parameters-The resulting structures were solvated in a cubic box with ϳ25,000 simple point charge (25) waters and counter ions added to ensure overall electrostatic neutrality. All of the simulations were performed with Gromacs 3.1.4 (26), with the ffgmx force field and the NVT ensemble at 300 K using the Berendsen thermostat. The pressure was coupled with the Berendsen barostat with a coupling constant, tau of 1.0 ps. Electrostatics were calculated using the particle mesh Ewald method (27,28) with a short range cut-off of 10 Å. The time step for integration was 2 fs. The LINCS algorithm (29) was used to restrain bond lengths. Each system was subjected to a 200-ps dynamics run with the protein restrained (10 kjmol Ϫ1 Å Ϫ2 on all heavy atoms). This was followed by 20 ns of free simulation. Except for RMSD calculations, only the last 18 ns was used for analysis. Hydrogen bonds were calculated with 30°angle and 3.5-Å distance cut-offs. All of the calculations were performed on a Linux cluster.

RESULTS AND DISCUSSION
Protein Motion-An initial evaluation of structural drift is provided by analysis of the C␣ atom RMSDs from the initial structures as a function of time. We omitted loops 1 and 2 (residues 413-450 and 486 -497, respectively) from this analysis to focus on the overall dynamics of the fold. Each simulation shows an initial jump of ϳ1.5 Å followed by slower relaxation (Fig. 2A). The closed cleft simulations, Gly, DS, DCS, ACPC, and Closed-Apo (only the Gly simulation is shown in Fig. 2A for clarity), all seem to converge to a similar RMSD (1.75 Å), indicating that the closed cleft form of the protein exhibits little conformational drift on these time scales. The RMSD of our simulations reveals that the open cleft forms of the protein (DCKA and Open-Apo simulations) exhibit a much larger conformational drift from the initial structure (up to 4 Å in the case of the Open-Apo simulation). Visual inspection of the Open-Apo simulation revealed that it appears to undergo a transition from an open to a more closed form of the protein between 5 and 8 ns. This is examined in more detail later.
In addition to conformational drift, we can also examine flexibility as a function of residue number. Eight residues in this region were disordered in the crystal structure of all but ACPC, presumably because of this high flexibility. In the complete receptor, loop 1 has been postulated to interact with a subunit on an adjacent dimer (9), which would stabilize it.
Interestingly, in the Open-Apo simulation, residues 533-544 and 663-737 show a larger degree of fluctuation compared with the other simulations. These residues are located in the D2 subdomain (residues 533-537 are in the intersubdomain strands). If one considers the C␣ root mean square fluctuation of this simulation after the cleft has closed (8 -20 ns; data not shown), the degree of fluctuation is reduced, implying that a closed cleft conformation stabilizes this region of the protein.
Unlike in other ionotropic glutamate receptors (iGluRs), in the NMDA receptors the extent of subdomain separation appears not to be linked to agonist efficacy. Both partial and full agonists confer essentially the same degree of subdomain closure on the ligand-binding domain (although antagonists continue to force the subdomains apart). Fig. 2C shows the intersubdomain separation (calculated as the distance between the centers of masses of the two subdomains) as a function of time. The closed cleft simulations (Gly, DS, DCS, ACPC, and Closed-Apo) all maintain a similar degree of domain closure (only DS is shown for clarity in Fig. 2C). Thus, at least on these time scales, partial agonists do not appear to lead to partial domain opening. The DCKA simulation shows a slight initial closure (about 1 Å over ϳ0.75 ns) to a separation of ϳ25.5 Å, which is maintained for the duration of the simulation. At 5 ns the Open-Apo simulation exhibits an initial rapid closing followed by a slower closure to approximately the same degree of domain separation as the DS simulation by ϳ17 ns.
Essential dynamics analysis of the Open-Apo simulation show that the mechanism of closing is comprised of two major motions. The first is a clamshell closing motion as shown in Fig. 3A. The second is analogous to a twisting of one side of D1 and D2 toward each other and is shown in Fig. 3B. These motions can occur concurrently. This has been observed before for this class of protein (30,31); thus, we are reasonably confident that this motion is genuine. Although subdomain closure has been observed before in simulations of these and other proteins (31,32) and evidence suggests that the ligand-binding domain exists in an open/closed equilibrium in the apo state (33), we are cautious not to overinterpret this motion. The simulation of the ligand-binding domain in the absence of the rest of the subunit may have led to a higher likelihood of a closure event caused by the lack of an opposing force that would be provided by the transmembrane ion channel domain.
Intersubdomain Contacts-In the crystal structures of NR1 there are a limited number of direct contacts formed across the binding left. Gln 405 interacts with Trp 731 directly and via a water with Asp 732 . Asn 499 and Gly 485 interact with Gln 686 , and Lys 493 interacts with Glu 712 . In addition to these, during the simulations, we commonly observe interactions between 13 residues in D1 and 10 residues in D2 ( Table 2). Perhaps the most interesting of these are between Thr 518 and both Ser 688 and Asp 732 , all of which form part of the ligand-binding pocket.
In the absence of a ligand in the binding cleft, it might be expected that the subdomains would open ready to accept a ligand. However, the Closed-Apo conformation that we created by removing the glycine ligand from Protein Data Bank structure 1PB7 does not undergo any significant increase in degree of subdomain closure throughout the   Water in the Binding Pocket-During the Open-Apo simulation, we observed protein motion that results in a conformation that is more closed than the Closed-Apo structure. However, the difference is small (less than 1 Å), and the extent of the closure is almost identical to that of the D-serine-bound structure (Fig. 2C). We are therefore confident that the conformation is genuine.
To characterize the closed conformation in the absence of a ligand more fully, we examined the nature of the interactions of water with the binding pocket. Because the Closed-Apo simulation is somewhat artificial, we used a section of trajectory (17-20 ns) from the Open-Apo simulations that corresponds to a closed form of the protein based on the degree of closure observed (Fig. 2C). 175 different hydrogen bonds with 55 unique waters are formed in this time frame, thus indicating that substantial water exchange is possible in the closed cleft form of the protein in the absence of ligand.
We looked at key waters in the binding pocket of the ligand-bound simulations in more detail. In particular we have considered the sites of three crystal waters (sites A, B, and C; see Fig. 4A) that were found in the same location in the glycine-, D-serine-, and D-cycloserine-bound structures and therefore might play a key role in ligand binding. Sites A and B are positioned between the ligand and D2, and C is located proximal to the ligand amino group between D1 and D2.
To define the sites of the crystal waters, we considered the charged atoms within 4 Å of each crystal water. Table 3 shows the number of unique waters that enter each of these sites during 20 ns of the Gly, DS, DCS, and ACPC simulations.
From these results it is clear that Gly undergoes far more water exchanges than either DS, DCS, or ACPC (although it is still a rare event). This is probably due to the larger size of D-serine, D-cycloserine, and ACPC.
In the DS simulation both waters A and B quickly leave the binding pocket and are not replaced in 20 ns. For Gly, waters at sites A and B are present but undergo several exchanges during the 20 ns. Exchanges for site A in the Gly simulation are shown in Fig. 4B. This is in contrast to the DCS simulation where waters are bound to both sites and undergo fewer exchanges but that do appear to coincide with ligand movement. Fig. 4C shows the water behavior corresponding to site A. The behavior of water in this site may also be linked to alternative orientations of DCS within the binding pocket (see Fig. 6B). Overall, this suggests that waters in sites A and B are not important in binding D-serine and that, at least in the case of Gly, the waters may act as a substitute to the H-bonding side chain of D-serine. There are no waters in site A in the ACPC crystal, but water in site B remains tightly bound with no exchanges throughout 20 ns. C remains tightly bound in all four simulations (data not shown).
The Ligand-binding Pocket-The binding pocket was defined by considering all residues in the D1 and D2 lobes that interact (interatomic distance Ͻ 3.2 Å as also described in Ref. 9) with the ligand in the crystal structures 1PB7, 1PB8, 1PB9, 1Y20, and 1PBQ. This results in a binding pocket defined by Pro 516 , Thr 518 , and Arg 523 , which reside in D1, and Ser 688 and Asp 732 , which are part of D2. In addition to these residues we identified seven new residues that are also capable of interacting with the ligand ( Table 4). Three of these (Gln 405 , Phe 484 , and Thr 486 ) are located in D1, while the others (Ser 687 , Val 689 , Trp 731 , and Ser 756 ) are in D2. None of these contacts were long-lived, but it is of interest that Ser 756 , which is in the hinge region, is capable of interacting directly with  For the analysis of the trajectories, the subsites were defined as protein atoms that were within 4 Å of the crystal waters. We then measured the distance of water oxygens to these center of masses. B shows results from the A subsite of the Gly simulation. The crystal water diffuses away from the binding pocket during the equilibration phase and is replaced by W388 (at t ϭ 940 ps), which is a crystallographic water located on the surface of the protein ϳ 28 Å away from subsite A in the crystal structure. W2121 and W4204 are waters added during the solvation process that diffuse into the binding pocket during the course of the simulation. C shows results for the A subsite for DCS simulation. W387 is the crystallographic water located in subsite A, and W319 is the crystallographic water originally observed in subsite B. the ligand because this may suggest a method through which the ligand can directly influence the hinge. The ligand remains in the binding site during all the simulations, although in both the glycine and D-cycloserine simulations there is movement of the ligand within the binding pocket that affects the hydrogen bonding patterns of the ligand. These ligands were able to adopt different orientations within the binding pocket that did not change the degree of subdomain separation (Fig. 2C). The transitions for both of these ligands were sharp and involved a loss of interaction between the ligand and Pro 516 .
When in the starting conformation (Fig. 5A), the glycine ligand interacts with Asp 732 , Ser 688 , Arg 523 , occasionally Pro 516 , and the amide group of Thr 518 . In addition, the side chain of Thr 518 interacts with Asp 732 . During the first ϳ13 ns of simulation, we observe two types of ligand movement before settling back into the crystal conformation. First, a rotation of the amino group so that it faces D2, during which the carbonyl group continues to interact with Thr 518 (Fig. 5, B and D). Second, we observe a complete flip where the ligand acts as a rigid body so that both amino and carbonyl groups are facing D2 (Fig. 5C). When this occurs contacts between the ligand and Thr 518 are broken, and often a water molecule moves into the gap between Thr 518 and the ligand, preventing Thr 518 interacting with Asp 732 in D2.
In the DCS simulation the ligand acts as a rigid body. In the starting conformation D-cycloserine interacts with the amino groups of Thr 518 and Arg 523 via its carbonyl oxygen and with Pro 516 , Asp 732 , and the side chain of Thr 518 via its amino group (Fig. 6A). As in the glycine-bound structure, the Thr 518 side chain also interacts with Asp 732 . When the ligand moves at ϳ4 ns, it continues to interact with all residues in the binding pocket except with Pro 516 . The amino group of the ligand is able to interact with Ser 688 , but the positioning of the ligand means that Thr 518 and Asp 732 no longer interact directly but via the ligand instead (Fig. 6B). This flipped conformation persists to the end of the 20-ns simulation.
Hydrogen Bonds in the Hinge-The hinge region of the ligand-binding domain was examined in closer detail as differences between the   conformation adopted in the presence of full and partial agonist or antagonist were reported. The hinge region is comprised of two strands (residues 537-539 and 750 -755), and interstrand hydrogen bonds in this region differ according to the conformation adopted by the hinge (9).
One key hydrogen bond forms between the carbonyl group of Leu 538 and the amino group of Phe 753 . This is present in all NR1 crystal structures to date and is independent of the ligand bound. The second interstrand hydrogen bond is variable. In the presence of the full agonists glycine and D-serine, Leu 538 forms a second hydrogen bond, between its amino group and the carbonyl group of Phe 754 (similar to that shown in Fig. 7A). However, in the presence of partial agonists ACPC, ACBC, and cycloleucine and the antagonist DCKA, the amino group of Leu 538 forms a hydrogen bond with the carbonyl group of Phe 753 (similar to that shown in Fig. 7B). Interestingly, the crystal structure with the partial agonist D-cycloserine bound (Protein Data Bank code 1PB9) has a hinge hydrogen bond pattern like those seen in the presence of full agonists (Fig. 7A).
In all ligand-bound cases (except DCS) the hinge conformation was unchanged from the starting structure. As might be anticipated in the DCS simulation, the hinge stably adopts a more partial agonist-like conformation with the carbonyl group of Phe 753 facing across the interstrand region (Fig. 7B). This conformation occurs within ϳ5 ns and persists until the end of the 20 ns simulation. This rotation appears to be closely associated with the observed movement of D-cycloserine in the binding pocket. As can be seen in the plot of the backbone dihedral angle of Phe 753 (Fig. 7C), rotations around the backbone occur frequently in the first 5 ns of the simulation. Stabilization occurs within 1 ns of the ligand changing binding mode.
These observations provide us with an indication of how partial agonists may confer their mode of action: partial agonists are able to move within the binding cleft and adopt different binding modes, and our results suggest that some of these binding modes affect the conformation of the hinge, whereas others do not. It is thus possible to see how a spectrum of partial agonists might act depending on how flexible they are and how readily they can influence the conformation of the hinge. Full agonists, on the other hand, might exhibit considerable movement (as we observe here for glycine) within the binding pocket but at the same time exert no effect on the hinge region.
Conclusions-We have demonstrated here that the dynamics of the NR1 subunit of the NMDA receptor are an important consideration with respect to their mode of action and in particular how partial and full agonists might be distinguished based upon their influence of the local conformations of the binding pocket. In addition, large scale conformational changes were also observed that provide a plausible pathway for a transition between the open cleft and closed cleft forms of the domain.
The Open-Apo transition and the general motion of the trajectories as decomposed by the eigenvector analysis reveals the principal motions of the ligand-binding core to be a lobe twisting motion and a hinge motion in agreement with previous observations (30,31). The speed of the closing transition observed during the Open-Apo simulation is fast (nanosecond time scale) and is currently beyond the resolution of experimental techniques for single molecules. This may be important in the context of protein ensembles (21); for example, recent NMR experiments have indicated that such considerations might be especially relevant for enzyme catalysis (34).
The closed cleft simulations showed little drift in terms of RMSD ( Fig.  2A) and thus provided a good opportunity to analyze ligand-protein and binding pocket-water interactions. The movement of the ligands within the site indicates that there is potential for the protein-ligand system to adopt multiple conformations that induce the same degree of domain closure. For the NR1 subunit, unlike the AMPA receptor subunit GluR2, recent crystallographic evidence has suggested that partial agonism is not related to the extent of cleft closure but rather to the efficacy of the ligand to induce a change in the conformation of the hinge region. In the crystal structure of NR1 complexed with the partial agonist D-cycloserine, the conformation of the hinge region was found to be "full agonistlike." Our simulations here show that D-cycloserine may in fact be capable of inducing the change in backbone hydrogen bonding of the intersubdomain strands that was found for other partial agonists crystallized (10). We also show that partial agonists are able to maintain full domain closure in the NR1 subunit and do not lead to partial domain opening. It is also interesting to note that although the glycine ligand was able to move substantially (changing hydrogen bonds with the binding pocket), it did not induce the change in the interdomain strand FIGURE 7. A, hydrogen pattern at the start of the simulation between the Leu 538 and Phe 753 /Phe 754 . In this conformation the Leu 538 amino group is hydrogen-bonded to the Phe 754 carbonyl group. B, the change in hydrogen bonding after the D-cycloserine had flipped (at t ϭ 4 ns; see Fig. 5B). The Leu 538 amino group is now hydrogen-bonded to Phe 754 carbonyl group. C, the behavior of the torsions of Phe 753 shows a distinct transition for the phi torsion angle (dark gray) and the psi torsion angle (light gray). There is also a distinct transition for the phi torsion angle of Phe 754 (black). backbone and thus remains consistent with the profile of a full agonist. Thus, we suggest that molecular dynamics simulations of the compounds within the receptor will be a useful tool in predicting agonist efficacy, by examining the influence of ligands on the local conformational dynamics of the hinge region.