INTRODUCTION

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION CONCLUSIONS REFERENCES

Fast synaptic transmission between nerve cells in mammals is carried out predominantly by ionotropic glutamate receptors (iGluR). These receptors are a family of ligand-gated ion channels that open in response to the binding of glutamate (Dingledine et al., 1999; Sprengel et al., 2001). Glutamate is released presynaptically and binds to a post-synaptic receptor gating a cation-selective channel, thus depolarizing the post-synaptic cell. Although glutamate is the natural ligand, the various iGluRs identified by sequence comparisons may also be classified in terms of their agonist pharmacology (Hollmann and Heinemann, 1994). Those receptors (GluR1-4) that show greatest sensitivity to the synthetic agonist -amino-3-hydroxy-5-methyl-4-isoxazole propionic acid (AMPA) are termed AMPA receptors (Borges and Dingledine, 1998). Likewise, those that show greatest sensitivity to kainate (GluR5-7 and KA1-2) are referred to as kainate receptors (Lerma et al., 1997; Chittajallu et al., 1999). Receptors activated by the synthetic N-methyl-D-aspartate (NMDAR1 and NMDAR2a-d) are called NMDA receptors (Yamakura and Shimoji, 1999) and need glycine and glutamate as the natural agonist. In all of these GluRs the agonist/antagonist binding site is located within the extracellular (EC) region of the protein (Fig. 1). Preceding the glutamate-binding domains is an amino-terminal domain (ATD). This has ~16% sequence identity to the bacterial leucine/isoleucine/valine-binding protein (O'Hara et al., 1993). Although the ATD is not directly involved in ligand binding, it has been shown to be important for subtype-specific interactions (Leuschner and Hoch, 1999). The polypeptide chain then proceeds to make up most of domain S1 of the ligand-binding site before forming two transmembrane (TM) helices (M1 and M2) and an intervening P-loop, suggestive of an inverted potassium channel TM architecture (Panchenko et al., 2001). As it leaves M2 the polypeptide chain then forms most of domain S2 in the ligand-binding cleft before forming a third TM helix, M3, followed by a short C-terminus. The x-ray structure of a water-soluble construct corresponding to the ligand-binding domains (S1S2) of the EC region of GluR2 was first solved by Armstrong et al. (1998), revealing the agonist-binding site to lie at the interface between S1 and S2.

View larger version (34K): [in this window] [in a new window]   FIGURE 1   (A) Overall topology of the Glu-R2 protein. ATD represents the amino-terminal domain (thought to be involved in tetramerization of the receptor). S1 and S2 are the two domains that make up the region that binds the ligand (L), i.e., glutamate. M1, P, M2, and M3 make up the transmembrane region of each subunit. (B) X-ray structure of the S1-S2 construct (Armstrong and Gouaux, 2000). The linker region is indicated by a black rod. The bilobal domain architecture is marked by the dashed line, with the S1 domain above and the S2 domain below this line. The ligand-binding site is shown occupied by glutamate. Figure drawn with Molscript (Kraulis, 1991) and POV-ray (http://www.povray.org).

Evidence from single-channel recordings on homomeric GluR4(i) receptors showed that currents elicited by kainate as the agonist were up to eightfold smaller than those elicited by AMPA or glutamate (Smith et al., 2000). Hence, kainate may be thought of as a partial agonist of the AMPA receptor. Recently, a series of high-resolution crystal structures of different agonists and of an antagonist (DNQX) bound to the GluR2 ligand-binding domain was published (Armstrong and Gouaux, 2000). Comparison of these structures of the AMPA- and glutamate-bound complexes revealed a domain closure of ~20° relative to the apo state. In contrast, the partial agonist kainate induced a domain closure of only ~12° .

Although it is evident that inter-domain movement occurs upon ligand binding, x-ray structures provide only static snapshots of a dynamic process. In this paper we use molecular dynamics (MD) simulations to explore changes in the conformational dynamics of the GluR2 ligand-binding construct in response to different patterns of binding-site occupancy, i.e., no ligand versus full agonist (glutamate) versus partial agonist (kainate). These simulations provide an insight into agonist dependency of the dynamics of the ligand-binding construct on a multi-nanosecond timescale.

	METHODS

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION CONCLUSIONS REFERENCES

Crystal structures (Armstrong and Gouaux, 2000) for the glutamate (1FTJ), kainate (1FWO), and the apo (1FTO) states of the S1S2 construct were downloaded from www.rcsb.org. The 1FTJ coordinate set of the glutamate-bound state was used to generate an additional apo state (hereafter referred to as Apo2) from which the glutamate was removed and the binding site solvated with five additional water molecules (Table 1). We use the same residue numbering as that employed by Armstrong and Gouaux (2000) for ease of comparison. Missing residues were added using the molecular editor in Quanta (Accelrys, San Diego, CA) and the N- and C-termini were acetylated and amidated, respectively, to mimic the continuation of the protein chain. Before MD simulations were started, we performed pK_A calculations (Adcock et al., 1998) to determine whether any of the binding-site residues were likely to adopt nonstandard ionization states. On the basis of the results of these calculations, all residues were modeled in their default ionization states.

Crystallographic waters within 4 Å of the protein were retained. Each of the three starting structures (glutamate-bound, kainate-bound, and apo) were solvated by a box of ~11,000 simple point charge (SPC) (Hermans et al., 1984) water molecules (Table 1) and the appropriate number of counterions were added to ensure overall neutrality of the system. The protein and ligands were then subjected to a restrained run of 200 ps, whereby the protein (and ligand if present) were harmonically restrained with a force of 1000 kJ mol¹. After this period all restraints were removed and the simulations run for 2 ns. For simulation of Apo2 this was extended to 5 ns. Electrostatics were calculated using particle mesh Ewald (PME) (Darden et al., 1993; Essman et al., 1995) with a 10-Å cutoff. The temperature was coupled with the Berendsen thermostat (Berendsen et al., 1984), at 300 K with _T = 0.1 ps. The time step for integration was 2 fs, and coordinates and velocities were saved every 5 ps. The LINCS (Hess et al., 1997) algorithm was used to restrain bond lengths. All simulations were performed using GROMACS (http://rugmd0.chem.rug.nl/~gmx/index.html) (Berendsen et al., 1995) running on a Silicon Graphics Origin 2000 computer.

	RESULTS

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION CONCLUSIONS REFERENCES

Domain closure and inter-domain motions

Comparison of the crystal structures of the S1S2 construct with and without agonists (Armstrong and Gouaux, 2000) suggested that the degree of domain closure between the two domains could be responsible for differences in the effectiveness of ligands as agonists. In our MD simulations we wanted to explore this aspect further in terms of the intra- and inter-domain motion and of how the motion of the domains is influenced by the ligand.

The drift of each simulated structure from its initial crystallographic conformation provides information on the quality of the simulations. The drift was measured in terms of the root mean square deviation (RMSD) of the C atoms from the initial structures as a function of time. For the three (Apo1, Glu1, and Kai) simulations after an initial rise (during the first ~300 ps) the C RMSD plateaus at ~1.25 Å (data not shown). This is indicative of a stable simulation. To obtain an initial estimate of the relative motions of the two domains during the 2 ns, we determined the C RMSD of the second domain (S2) after fitting the S1 domain C atoms to the initial structure(s) (Fig. 2 A). In this case the domain S2 C RMSDs for the Glu1 and Kai simulations plateau at ~1.5 Å after 300 ps. In contrast, the Apo1 simulation shows a C RMSD, which rises to ~3 Å and exhibits much greater fluctuations than for either of the liganded states.

View larger version (33K): [in this window] [in a new window]   FIGURE 2   (A) C RMSD of the S2 domain relative to the S1 domain as a function of time for the Apo1 (), Kai (· · ·), and Glu1 (- - - -) simulations. The increased fluctuation of the Apo1 state is evident. (B) Superimposed backbone traces (saved every 200 ps) of the C atoms for these three simulations. The difference in the movement of the S2 domain relative to the S1 domain between the simulations is clearly discernible.

This comparison can be extended by fitting the domain S1 protein coordinates (saved every 20 ps) on top of the initial structure(s) and examining the backbone traces of the protein (Fig. 2 B). From this diagram it is evident that there is a progressive increase in the motion of domain S2 (relative to domain S1) from Glu1 to Kai to Apo1. This suggests that the full agonist results in greater immobilization of the domains relative to one another than does the partial agonist. This is in addition to differences in the degree of (static) domain closure observed in the x-ray structures (Armstrong and Gouaux, 2000).

A further measure of the quality of our simulations is provided by comparison of temperature (B) factors calculated from simulation root mean square fluctuations (RMSFs) with those determined experimentally in the x-ray diffraction studies. We found good agreement of the calculated B-factors from the simulation with those reported in the crystallographic structures. A couple of residues located on the outer surface loops had values higher than those reported, but generally the simulated values were consistently lower than the experimental values by ~5 Ų. Given we observed a significant difference in inter-domain motion induced by different ligands, we wished to reassure ourselves that the mean conformation of the protein in our simulations retained the domain closure seen in the crystal structures. Therefore, we determined inter-domain distance matrices averaged over the whole duration of each simulation (Fig. 3). These quite clearly reveal the gradient in inter-domain contacts between the Apo1 simulation (fewest contacts) and the Glu1 simulation (most contacts). The mean distances of these contacts are in accord with the mean degree of domain closure observed in the crystal structures, indicating that the average behavior in the MD simulations reproduces the crystallographic averages. As an additional check we compared the radius of gyration (of the C atoms) for each of the simulations to provide a measure of the overall compactness of the protein. A clear distinction can be seen between the mean radius of gyration (R_GYR) for the three different simulations with Apo1 (R_GYR  19.2 Å) > Kai (R_GYR  18.7 Å) > Glu1 (R_GYR  18.4 Å). Furthermore, the fluctuations in the radius of gyration are much greater for Apo1, again indicative of significant inter-domain motion. Thus, this analysis supports the conclusion reached by visualization of the inter-domain motions.

View larger version (45K): [in this window] [in a new window]   FIGURE 3   (A) Mean inter-domain distance matrices for the Apo1, Kai, and Glu1 simulations. Each matrix represents the C-C distances averaged across the corresponding simulation. Interactions nearer than 4 Å are indicated by a white square. The increased number of inter-domain contacts in the Glu1 simulation is apparent. The major regions of contact in the Glu1 simulation are ringed.

Apo2 simulation

We have also performed a simulation (Apo2) starting from the conformation of the protein observed in the crystal structure with glutamate bound, but with the ligand removed and replaced by water molecules (see above). Thus the protein is in the domains-closed state, but without any stabilizing ligand. This simulation was run for 5 ns. Even on this timescale we did not see any evidence of net movement of the domains relative to one another to yield a conformation resembling that in Apo1. We suspect that much longer simulations (>100 ns) might be required to see such a conformational change, although domain motions have on occasions been seen in ~2-ns simulations (Roccatano et al., 2001).

We have applied the various analyses described above to the Apo2 simulation. The drift from the initial structure seen in the RMSD plots is similar to that observed for the Glu1 simulation, implying that the protein maintains a conformation that is similar to the glutamate-bound state. The radius of gyration for the Apo2 simulation is ~18.4 Å across the whole 5 ns, which is identical to the value observed for the Glu1 simulation.

Multiple peptide bond conformations

In their analysis of the crystal forms, Armstrong et al. (1998) observed a flip in the peptide bond between Asp651 and Ser652 for the AMPA-bound and in one of three glutamate-bound protomers. They suggested that the switch of this residue from an unflipped to a flipped conformation was related to receptor activation. Our initial simulation (Glu1) was based upon protomer A, which in the crystal was observed in an unflipped conformation. During the 2-ns Glu1 simulation we did not observe any change in the backbone torsion ( and ) angles in this region. We investigated this aspect further by running a simulation (Glu2; see Table 1) starting from protomer C, which has the flipped conformation of the peptide bond. Within 600 ps, we observed the  angle of Ser652 to rotate through ~50-60 ° (Fig. 4) and remain in that new conformation for the remainder of the simulation. The new conformation is closer to the protomer B conformation in the crystal. Thus, our simulations suggest that in the presence of glutamate, multiple conformations of the Ser652 can exist, as in the x-ray structure. Although it has been postulated that this flipping could be related to receptor activation it is very difficult to make any conclusion with respect to the mechanism for two reasons. The first is the timescale of this simulation with respect to the timescale of activation; such a mechanism of activation must involve conformation changes in the region of the linker, for to observe such a propagation of conformational change would require multiple long simulations. Second, it could be that the presence of the linker in this region alone infers different conformation changes than that which occur in the full-length receptor.

View larger version (35K): [in this window] [in a new window]   FIGURE 4   Backbone torsion angles as a function of time for simulation Glu2 (starting with monomer C from the crystal structure). (A)  of Ser652; (B)  of Gly653. In the crystal structure there are three monomers (see text) with corresponding torsion angles: Ser652A,  = +118°; Gly653A,  = +167°; Ser652B,  = +53°; Gly653B,  = 138°; Ser653C,  = +27°; and Gly653C,  = 125° .

The overall R_GYR for the Glu2 simulation was 18.4 Å, i.e., identical to that observed for the Glu1 simulation. Our analysis of the RMSD plots also indicated that the overall behavior of the protein in Glu2 was similar to that in the Glu1 simulation.

Binding modes

What is the main distinction between glutamate (a full agonist) and kainate (a partial agonist)? The answer must lie within their manners of binding. We thus examined in detail interactions at the binding sites in the Glu1 and Kai simulations. The binding site is comprised of an intricate network of hydrogen bonds (Fig. 5). It can be seen that the key interactions of domain S1 involving Arg485, Thr480, and Pro478 are seen with both ligands as are the key interactions of domain S2 involving Ser654, Thr655, and Glu705. We have monitored the existence and duration of these key H-bond interactions throughout the Glu1 and Kai simulations (Glu1 and Kai; Fig. 6). Distinct differences are observed between the two ligands in terms of the hydrogen-bond lifetimes.

View larger version (24K): [in this window] [in a new window]   FIGURE 5   Key interactions of the binding sites for the glutamate and kainate binding sites, respectively. The structures are the minimized ones before the start of the restrained run. The diagram was generated with the program LIGPLOT (Wallace et al., 1995).

View larger version (56K): [in this window] [in a new window]   FIGURE 6   Hydrogen bonds versus time for key ligand/protein interactions for the Glu1 (A) and Kai (B) simulations. (Atom numbers are defined in Fig. 5.)

For both ligands, Arg485 is observed to maintain the largest number of contacts for long lifetimes in agreement with the observation (Armstrong and Gouaux, 2000) that this residue is the primary anchor for the -carboxyl group of the ligand and with experiments indicating that mutations of this residue invariably lead to loss of function (Uchino et al., 1992; Wafford et al., 1995; Kawamoto et al., 1997). However, there are also some clear differences between the two ligands. In simulation, Glu1 H-bond analysis revealed two additional long-lived interactions, not present in the starting structure. These appeared after 1 ns and were between 1) the NH of Glu and the Thr480 backbone OG atom and 2) the glutamate O and the Arg485 side-chain NH. This was examined closer and was found to be due to flipping (by ~100°) around the C-C bond of the side chain of Thr480 (Fig. 7). The lengths and angles of these H-bonds are summarized in Table 2. Qualitatively, a change is seen whereby a protein/protein H-bond (Thr480-Arg485) is replaced by two H-bonds from these side chains to the ligand, glutamate. It should be noted that a similar flipping of Thr480 was not seen during the kainate simulation, possibly due to the presence of the ring in the ligand allowing less flexibility in the binding site.

View larger version (21K): [in this window] [in a new window]   FIGURE 7   Diagram of flipping of the side chain of Thr480 in the Glu1 simulation. At 0.5 ns the 1 (C-C-C-OG) torsion is 67 °. At ~1.5 ns it flips to a value of 165 °. Figure was drawn with VMD (Humphrey et al., 1996).

We have analyzed the energetics of these binding modes. The enthalpic energy difference between the two different binding modes for glutamate is not significant. An accurate value for the enthalpic binding energy between glutamate and kainate is not possible from these simulations, but an approximate calculation reveals that the difference is substantial (of the order of 100 kcal/mol).

	DISCUSSION

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION CONCLUSIONS REFERENCES

The main achievements of these simulations have been 1) to reveal a correlation between intra-domain mobility and receptor occupation, such that low mobility corresponds to occupancy by a full agonist, and 2) to show that binding of glutamate to the S1S2 protein may be strengthened beyond that present in the crystal by rotation of the Thr480 side chain. Thus, the simulations may reveal aspects of GluR activation additional to those observed in the x-ray structures.

The principal limitations of the simulations are their relatively short duration, a result of computational limitations for atomistic simulations of relatively large systems (~40,000 atoms). This short duration is problematic with respect to two distinct areas. The first of these is that the activation event (from ligand binding to channel opening) is of the order of milliseconds, and it is thus unwise to extrapolate beyond these current results. The second consideration is that of statistical sampling. We have performed reasonably long simulations, but only one flipping event around the Thr480 C-C bond is observed for example. Ideally one should run these simulations many more times to ensure greater confidence that the observed events are a significant property of the system. One would also like to be able to estimate how well conformation space has been sampled. We are quite certain that many more simulations would be needed to start to address this issue more fully.

Furthermore, even after extending the Apo2 simulation to 5 ns it was not possible to observe a conformational transition corresponding to domain opening in the absence of ligand. Comparisons of the different crystal structures provide clear evidence of domain opening/closure (Armstrong and Gouaux, 2000). However, applying similar analyses (Hingefind (Wriggers and Schulten, 1997) and DynDom (Hayward and Berendsen, 1998)) to any single simulation failed to reveal any movement of the domains that could be dissected into bilobal hinge motions. This also suggests that domain closure is on a much longer timescale than is accessible in current simulations. Indeed, a similar observation was made for nanosecond simulations of an SH2-SH3 domain system (Young et al., 2001), where domain movement was known to exist. Furthermore, it seems that rare events on this timescale may not be directly observable, but the collective motions that contribute to their behavior may well be (Tai et al., 2001).

The other major limitation, shared with the x-ray studies, is that our studies are of the water-soluble S1S2 construct rather than the intact GluR2 protein. However, as shown and discussed in Armstrong and Gouaux (2000), the binding curves for ³H-AMPA and the IC₅₀ values for ligand displacements are very similar to longer S1S2 constructs (Chen and Gouaux, 1997) and to those for the full-length receptor (Keinänen et al., 1990). Thus it seems likely that the pharmacological properties of the S1S2 construct parallel those of the GluR2 receptor per se. Intriguingly, although very speculative, we also analyzed the movement of the linker residues and surrounding regions. Although the RMSFs of the linker residues were very similar in all three simulations, the RMSF values of the surrounding residues ~20 residues either side of this linker region showed noticeably higher values for the apo simulation compared with the ligand-bound simulations. Remaining residues had indistinguishable RMSF values. In the full-length protein the linker region would connect to the transmembrane domain, and it is thus tempting to try to speculate how this may be related to channel gating, but this will require much longer simulation times to be able to investigate this with any confidence.

	CONCLUSIONS

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION CONCLUSIONS REFERENCES

Our main conclusion is that it is important to consider the conformational dynamics of a receptor and to determine how ligand binding may modulate the dynamics. In particular, for the GluR2 it seems that agonist binding not only results in domain closure (as revealed crystallographically) but also results in a decrease in domain mobility (as revealed in simulations) and that the reduction in mobility is greater for a full than for a partial agonist. At a more microscopic level our simulations reveal an important side-chain flipping motion of a residue (Thr480) in the binding site that changes the mode of binding of the agonist (glutamate). This side-chain flip generates an alternative mode of binding not observed for kainate and not observed for glutamate in the crystal structures. Such an observation may have important consequences in the consideration of drug design. However, the timescale of the inter-domain motion is significantly longer than that which is obtainable by conventional MD. We are currently extending this work by using techniques more suited to longer timescales.