Gating Motions and Stationary Gating Properties of Ionotropic Glutamate Receptors: Computation Meets Electrophysiology.

Ionotropic glutamate receptors (iGluRs) are tetrameric ligand-gated ion channels essential to all aspects of brain function, including higher order processes such as learning and memory. For decades, electrophysiology was the primary means for characterizing the function of iGluRs and gaining mechanistic insight. Since the turn of the century, structures of isolated water-soluble domains and transmembrane-domain-containing constructs have provided the basis for formulating mechanistic hypotheses. Because these structures only represent sparse, often incomplete snapshots during iGluR activation, significant gaps in knowledge remain regarding structures, energetics, and dynamics of key substates along the functional processes. Some of these gaps have recently been filled by molecular dynamics simulations and theoretical modeling. In this Account, I describe our work in the latter arena toward characterizing iGluR gating motions and developing a formalism for calculating thermodynamic and kinetic properties of stationary gating. The structures of iGluR subunits have a highly modular architecture, in which the ligand-binding domain and the transmembrane domain are well separated and connected by flexible linkers. The ligand-binding domain in turn is composed of two subdomains. During activation, agonist binding induces the closure of the intersubdomain cleft. The cleft closure leads to the outward pulling of a linker tethered to the extracellular terminus of the major pore-lining helix of the transmembrane domain, thereby opening the channel. This activation model based on molecular dynamics simulations was validated by residue-specific information from electrophysiological data on cysteine mutants. A further critical test was made through introducing glycine insertions in the linker. Molecular dynamics simulations showed that, with lengthening by glycine insertions, the linker became less effective in pulling the pore-lining helix, leading to weaker stabilization of the channel-open state. In full agreement, single-channel recordings showed that the channel open probability decreased progressively as the linker was lengthened by glycine insertions. Crystal structures of ligand-binding domains showing different degrees of cleft closure between full and partial agonists suggested a simple mechanism for one subtype of iGluRs, but mysteries surrounded a second subtype, where the ligand-binding domains open to similar degrees when bound with either full or partial agonists. Our free energy simulations now suggest that broadening of the free energy basin for cleft closure is a plausible solution. A theoretical basis for these mechanistic hypotheses on partial agonisms was provided by a model for the free energy surface of a full receptor, where the stabilization by cleft closure is transmitted via the linker to the channel-open state. This model can be implemented by molecular dynamics simulations to predict thermodynamic and kinetics properties of stationary gating that are amenable to direct test by single-channel recordings. Close integration between computation and electrophysiology holds great promises in revealing the conformations of key substates in functional processes and the mechanisms of disease-associated mutations.

Introduction

Ionotropic glutamate receptors (iGluRs) are tetrameric ligand-gated ion channels. AMPA and NMDA receptors (AMPARs and NMDARs) are two main subtypes of postsynaptic iGluRs, and mediate the vast majority of fast excitatory neurotransmission in the vertebrate nervous system. These receptors convert transient glutamate signals, arising from presynaptic release, into postsynaptic electrical and biochemical signals. The duration of glutamate-induced channel opening determines the strength of synaptic connections between neurons and is critical to neurodevelopment and higher order processes such as learning and memory.[1] Recently, numerous missense mutations of iGluRs, in particular NMDARs, have been associated with an array of neurological disorders (e.g., autism, epilepsy, intellectual disability, and schizophrenia).[2,3] Thus, deep understanding of iGluR structure and function is critical to both basic and clinical neuroscience.

Because of their physiological importance and pathophysiological implications, iGluRs have been under intense investigations. AMPAR subunits GluA1–4 form both homo- and heterotetramers, whereas NMDARs form obligatory heteromers consisting typically of two GluN1 and two GluN2A-D subunits. Each subunit harbors an extracellular agonist-binding site. Glutamate is the agonist for all AMPAR subunits and GluN2 subunits, whereas glycine is the agonist for the GluN1 subunit. The central mechanistic question is how agonist binding triggers the opening of the ion channel in the cell membrane. For several decades before the turn of the century, electrophysiology was the primary means for characterizing the function of iGluRs and gaining mechanistic insight, both at synapses and on reconstituted or expressed receptors.[4,5] Whole-cell recordings of AMPARs show fast activation and deactivation as well as rapid and strong desensitization. In contrast, NMDARs show relatively slow gating kinetics, with weak or no desensitization depending on the specific subunits (e.g., GluN2A versus GluN2C).

Single-channel recordings of AMPARs showed a channel-closed state and up to four subconductance levels that were roughly separated by quantized increments.[6−9] A simple kinetic model was proposed, where transitions occurred only between adjacent conductance levels and occupancies of these levels were binomial.[7] According to this model, the four subunits open and close the channel independently and contribute additively to the unitary current. However, more complex behaviors were observed recently, including deviation from the binomial distribution at subsaturating agonist concentrations,[8] varying subunit open probabilities in different segments of a given current trace, and two to three substates within each conductance level.[9] NMDARs again exhibit very different behaviors in single-channel recordings.[10−13] A single conductance level was observed, but lifetime distributions revealed multiple components within both the channel-closed and open states. The prevailing model for kinetic analysis of single-channel currents assumes five closed substates and two open substates; three of the closed substates (C3, C2, and C1) and the two open substates (O1 and O2) are connected sequentially, while two longer lived substates (C5 and C4) branch off from C3 and C2, respectively.[12,14−19]

It has long been known that iGluR subunits have a modular architecture containing four distinct structural domains (Figure ): extracellular amino-terminal (ATD) and ligand-binding (LBD) domains; a transmembrane domain (TMD) forming the ion channel upon tetramer assembly, with the M3 helix lining the channel pore and harboring the activation gate near the helix C-terminus; and an intracellular, disordered C-terminal domain (CTD).[20] The first structure of an isolated LBD, that of GluA2 bound with a partial agonist, was determined in 1998 by X-ray crystallography, identifying the ligand binding site within the cleft between two subdomains (referred to as D1 and D2).[21] Since then crystal structures of GluA2 LBDs bound with an assortment of ligands have been deposited in the Protein Data Bank (PDB), reaching an astounding total of 153 entries at present. Interestingly, AMPAR full agonists induce tight closure of the LBD clefts, partial agonists induce intermediate closure, and competitive antagonists lead to wide clefts, suggesting a correlation between degree of LBD closure and agonist efficacy.[7,22] Another interesting observation is that ligand size apparently dictates the degree of AMPAR LBD closure (Figure ): small ligands, usually agonists, fit snugly in a tightly closed cleft, whereas large ligands, usually antagonists, pry open the cleft.[23] In 2009, the first crystal structure of a CTD-truncated GluA2 AMPAR bound with a competitive antagonist was published (PDB 3KG2; Figure ).[24] This structure confirmed the highly modular architecture of iGluR subunits, with well-separated ATD, LBD, and TMD connected by flexible linkers. By superimposing the D1 subdomains of an isolated LBD bound with glutamate, it was possible to imagine how agonist-induced D2 closure upon D1 could lead to an outward movement of the M3 helix.

Figure 1

Modular architecture of iGluRs. (A) Domains within a GluA2 subunit (CTD missing). Subdomains of the ATD (L1/2) and LBD (D1/2) and transmembrane helices of the TMD (M1/3/4) are indicated; a bound ligand is shown in space-filling mode and the M3-D2 linker is in green. Inset: enlarged view of the region around the linker. (B) Subunit organization within the homotetramer (PDB 3KG2).[24] The A, B, C, and D subunits are in marine, magenta, light blue, and light pink, respectively. In NMDARs, the A/C and B/D subunits are GluN1 and GluN2, respectively.

Figure 2

Ligand size as a determinant for the degree of GluA2 LBD closure and agonist efficacy. (A) Superposition of an agonist-bound GluA2 LBD structure (PDB 1M5C; red) and an antagonist-bound structure (PDB 3R7X; blue). The D1 subdomain is used for superposition. Cα atoms of K410 in D1 (dark red) and K695 in D2 (light red), used for defining the distance dK410–K695, are shown as yellow spheres. (B) Histograms of AMPAR agonists, partial agonists, and competitive antagonists binned according to dK410–K695. (C) Correlation between dK410–K695 and ligand molecular weight, with R2 = 0.58. Reproduced with permission from ref (23). Copyright 2012 Elsevier.

At present the PDB contains 44 entries for isolated LBDs of NMDARs, but surprisingly, they have similar degrees of closure whether bound with full or partial agonists (although competitive antagonists still induced wide cleft opening).[25−28] That the crystal structures do not offer a clue to NMDAR partial agonism is a manifestation of their limits in providing mechanistic insights. There are now also 44 PDB entries of crystal and single-particle electron cryomicroscopy (cryoEM) structures for CTD-truncated AMPARs and NMDARs (for a recent review, see ref (29)). These structures putatively represent different functional states, but while the water-soluble domains (ATD and LBD) display ligand-induced structural changes, none of them has captured an open conformation for the ion channel in the TMD layer. The inability to capture the open channel could be due to the strong influence of the solubilizing environment on TMD structures and the fact that solubilizing conditions in sample preparations for crystal and cryoEM structure determination often do not model well the key biophysical properties of cell membranes.[30] If it is challenging for structural techniques to capture stable functional states, then the challenge is much greater for shedding light on conformations of the kinetic substates that have been identified by single-channel recordings. This lack of atomic-level structural information has been a tremendous impediment to progress in iGluR physiology and pathophysiology.

In this Account, I describe our recent studies to highlight the contributions of molecular dynamics (MD) simulations and theoretical modeling in filling some of the gaps in knowledge regarding iGluR structure and function. In particular, based on MD simulations we have developed an activation model in which the M3-D2 linkers play crucial and subunit-specific roles.[18,31,32] Our free energy simulations have supported the hypothesis that, whereas a reduced degree of LBD closure, corresponding to a shift in the minimum position of the LBD free energy basin, may underlie AMPAR partial agonism, a reduced curvature of the LBD free energy basin may lead to NMDAR partial agonism.[33] To further interrogate this hypothesis, we have formulated a theoretical model for calculating iGluR thermodynamic and kinetic properties during stationary gating.[34] This theoretical model points to the exciting prospect that iGluR functional properties can be predicted from intra- and interdomain energetics and dynamics in MD simulations. Guided by this theoretical model, most recently we have used free energy simulations to explore conformations of NMDAR kinetic substates.[35]

Simulations of Gating Motions

As already noted, crystal structures of both isolated domains and CTD-truncated constructs have been invaluable in formulating mechanistic hypotheses. However, it should also be recognized that these structures only represent sparse, often incomplete snapshots during iGluR activation. MD simulations can allow for two conformational snapshots to be connected through a physically plausible path, thereby molecular motions from one functional state to another can be explored. One such simulation technique is called targeted MD, by which one region of a protein system is constrained to move from one conformation to another.[36] In our study of GluA2 AMPAR activation, we mimicked agonist binding by forcing the four LBDs to move from the cleft-open to the cleft-closed conformation (Figure A).[31] As a direct consequence of the D2 closure upon D1, the M3-D2 linkers moved outward, more so in the B/D subunits than in the A/C subunits. The different extents of outward movements were attributed to the near orthogonal orientations of the M3-D2 linkers, parallel and perpendicular, respectively, to the membrane plane in the B/D and A/C subunits (Figures and 3A). The unevenly splayed M3-D2 linkers in turn pulled on the C-termini of the M3 helices, resulting in channel opening and symmetry breaking of the TMD layer from 4-fold to 2-fold.

Figure 3

Activation model of GluA2 AMPAR and validation. (A) Model from MD simulation. The A/C and B/D subunits are in two shades of red and blue, respectively. The D1 and D2 subdomains are represented as ovals, the M3 helices as rods, and the M3-D2 linkers as lines; the latter are critial in coupling LBD and TMD motions. (B) Validation by substituted cysteine modification rates of Sobolevsky et al.[39] The functional data (blue bars) are ln(kM+/kM–), where kM+ and kM– are modification rates in the presence and absence of glutamate; arrows on top indicate that the bars present lower bounds. The MD simulation results (red bars) are differences in solvent accessible area (in Å2) of residues between the activated and resting states. Adapted with permission from ref (31). Copyright 2011 Nature Publishing Group.

We recognized that functional data from whole-cell recordings of cysteine-substituted mutants could provide a variety of residue-specific information for validating mechanistic models. In particular, correlation between agonist-induced changes in rates of modifying substituted cysteines by methanethiosulfonate (MTS) reagents and changes in solvent accessible areas of the corresponding residues in key gating elements is very useful for testing activation models.[31,32,37,38] Sobolevsky et al.[39] found that, upon AMPAR activation, cysteines substituted into the M3 C-terminal residues L610, I613, S614, T617, N619, A621, and F623 had increased modification rates, whereas L620 had a decreased rate. As shown in Figure B, these changes are reproduced well by the variations in accessible area during our simulated activation process, with a correlation coefficient of 0.84. In particular, A621 is pore-facing and became more accessible upon channel activation, whereas L620 projects outward and faces the pre-M1 helix in the same subunit. As M3 expanded outward, the spacing between L620 and the pre-M1 helix was reduced and hence L620 became less accessible. Moreover, the 2-fold symmetry of the open channel suggested by our simulations is consistent with Cd2+ coordination data on substituted cysteines at A621.[40] One caveat is that different AMPAR isoforms were studied by MD simulation (GluA2) and electrophysiology (GluA1), so a more detailed comparison is not warranted.

Our follow-up study found similar gross features in gating motions for the GluN1/N2A NMDAR.[32] However, by its heteromeric nature, the unequal contributions of the two subunit types came into focus. Our MD simulations showed that LBD closure led to greater outward movement of the GluN2 M3 helices than the GluN1 counterparts. A simple explanation for this difference is that the GluN1 and GluN2 subunits align to the A/C and B/D subunits, respectively, of the GluA2 AMPAR,[41] with near orthogonal orientations of the M3-D2 linkers in the two subunit types. The simulation results fit with a wealth of electrophysiological data indicating that the GluN1 and GluN2 subunits play unequal roles during channel activation.[42−46] In particular, in the open state, substituted cysteines in the GluN1M3 helices overall showed higher modification rates (thus implicating greater access) by MTS reagents than those in the GluN2M3 helices. In our activation simulations, the greater outward movement of the GluN2 M3 helices led to a rhombus shape for the positions of four homologous M3 residues, with the GluN1 residues at the obtuse-angle vertices and the GluN2 residues at the acute-angle vertices. The simulations thus revealed that the electrophysiological observation of greater GluN1 M3 exposure unexpectedly results from greater GluN2 M3 outward movement, thereby implicating a stronger contribution of the GluN2 subunits to channel activation. It should be notated that the initial structure for MD simulations in this study was built by homology modeling to the GluA2 AMPAR, before the first crystal structures of CTD-truncated NMDARs were published.[47,48] Recently Dutta et al.[49] performed elastic network modeling on the respective crystals structures of an AMPAR and NMDAR, and found similar global modes of motion (with minor subtype-specific differences), in line with our conclusion that the two subtypes of iGluRs share similar gross features in gating motions.

As a further test of the critical roles of the M3-D2 linkers, especially the ones on GluN2A, in transmitting the stabilization provided by agonist-induced LBD closure to the open channel, we introduced glycine insertions into the linkers.[18] Our MD simulations showed that, with lengthening by glycine insertions, the linker became less effective, more so with GluN2A insertions than with GluN1 insertions, in pulling the M3 helices, leading to weaker stabilization of the channel-open state. In full agreement, single-channel recordings showed that the channel open probability decreased progressively as the linkers were lengthened by glycine insertions. This collaboration between computation and electrophysiology enabled us to gain a much deeper understanding on the mechanical coupling between the LBDs and TMD tetramer by the linkers than either approach alone would, with significant implications for further advances (see below).

Free Energy Simulations

The conformational space of proteins can also be explored via free energy simulations. Previous computational studies have confirmed that crystal structures of isolated AMPAR LBDs bound with agonists spanning a range of efficacies and antagonists correlated well with the minimum positions of the free energy basins for cleft closure,[50] thereby lending support to the hypothesis that AMPAR partial agonism arises from partial agonists not inducing as much cleft closure as full agonists.[7,22] However, that still left the partial agonism at NMDARs as a mystery, since their LBDs bound with full and partial agonists have similar degrees of closure.[25−28] We reasoned that partial agonists must perturb the free energy basin of a full agonist-bound LBD in some way; if not the minimum position, then the next suspect is the broadness, or curvature, of the basin.[33]

To test the latter hypothesis, we computed the free energy surfaces of the GluN1 LBD bound with a variety of ligands (Figure A) using umbrella sampling (Figure B,C). The free energy minima were similarly positioned for LBDs bound with full and partial agonists (in line with their crystal structures), but partial agonists significantly broadened the free energy basin (i.e., reduced the curvature), precisely as we hypothesized (Figure D,E). Antagonists both shifted the minimum position and further reduced the curvature. The broadening of free energy basin was directly illustrated in the simulations. In umbrella sampling windows where the LBD cleft was restrained to semiclosed positions, where the degree of cleft closure is intermediate between those in crystal structures of LBDs bound with a small full agonist and a large antagonist, the glycine-bound cleft retracted to more closed conformations, i.e., toward the free energy minimum, but the clefts bound with partial agonists exhibited less tendency of retraction, and those bound with antagonists showed tendency of further opening. In other words, compared to the full agonist glycine, partial agonists were more likely to maintain the LBD in semiclosed conformations, and antagonists pushed the LBD toward open conformations.

Figure 4

Broadness of LBD free energy basin as a determinant of NMDAR partial agonism. (A) Plot of the relative efficacies of GluN1 ligands[5] against their volumes. Chemical structures are displayed for the full agonists, glycine (Gly) and d-serine (DSN); partial agonists, d-cycloserine (DCS), 1-aminocyclopropane-1-carboxylic acid (ACPC), and 1-aminocyclobutane-1-carboxylic acid (ACBC); and competitive antagonists, 1-aminocyclopentane-1-carboxylic acid (CLE) and 5,7-dichlorokynurenic acid (DCKA). Linear regression (red line) for six ligands (other than DCKA) has R2 = 0.72. (B) Two distances, ξ1 and ξ2, used to describe intersubdomain motions; D1 and D2 are in green and cyan, respectively. (C) Free energy surface of the Gly-bound GluN1 LBD, presented as contours. The free energy minimum is highlighted by a white dot; a white diagonal line passing through the minimum is shown. Stars in white, gray, and black represent crystal structures bound with Gly (PDB 1PB7), CLE (PDB 1Y1M), and DCKA (PDB 1PBQ), respectively. (D) Overlay of free energy contours for the Gly-bound case (black) and the indicated non-Gly counterpart (red). (E) Slices of the free energy surfaces along the diagonal coordinate, ξ1′ = (ξ1 + ξ2)/2, for the full and partial agonists (solid curves) and antagonists (dashed curves). Reproduced with permission from ref (33). Copyright 2015 Elsevier.

The simplest explanation for these differences among the full and partial agonists and antagonists in the probabilities of propelling the LBD into semiclosed and open conformations is ligand size. Recall that for AMPARs, ligand size affects efficacy apparently by selecting the degree of cleft closure at the LBD free energy minimum (Figure ). For NMDARs, ligand size is evidently also a determinant of efficacy (Figure A), not through shifting the minimum position but by changing the broadness of the free energy basin. We thus conclude that NMDAR partial agonism arises from partial agonists broadening the free energy basin, making the LBD more readily transition from cleft-closed to cleft-semiclosed conformations and thereby provide less stabilization to the channel-open state.

Very recently we carried out a follow-up study on both GluN1 and GluN2A LBDs,[35] partly to test the assumption of Kussius and Popescu[16] that disulfide cross-linking would trap the LBDs in cleft-closed conformations. Our free energy simulations showed that, instead of being locked in the fully closed conformation, the cross-linked LBDs sampled semiclosed conformations almost as readily as the agonist-bound LBDs. For the GluN1 LBDs, the free energy minimum of the cross-linked form shifted slightly toward a lower degree of closure, whereas the free energy curvature was similar to that of the agonist-bound form. On the other hand, for the GuN2A LBDs, not only did the free energy minimum of the cross-linked form shift toward a higher degree of closure, but also the free energy curvature increased markedly. These results suggest that the cross-linking on the GluN1 subunits would mostly preserve the gating behavior of the agonist-bound receptor, but the cross-linking on the GluN2A subunits could result in a measurable increase in channel opening. These expectations were consistent with functional data from single-channel recordings.[16]

One should always be careful about whether conclusions drawn from studies of isolated domains apply to the full receptor. Although modular architecture is a hallmark of iGluR structures, the LBDs form A/D and B/C dimers via their D1 subdomains (Figure B) and, for NMDARs in particular, the LBDs and ATDs have extensive interactions.[47,48] The D1/D1 interface is very important for gating kinetics.[19] Likewise, ATD deletion has a strong effect on NMDAR gating properties including channel open probability, although the receptors remain functional.[51,52] Future work should therefore examine the impact of interdomain and intersubunit couplings on the free energy landscapes of individual domains. A similar issue in elastic network modeling was addressed by a “subsystem-environment” approach.[49] Next I turn to a very different approach.

Theoretical Modeling

The idea that both the minimum position and curvature of the LBD free energy basin are determinants of agonist efficacy is physically sound and provides plausible explanations for AMPAR and NMDAR partial agonisms. Yet how the energetics of LBD closure is transmitted to the channel is not entirely clear. What is needed is a model for the free energy function of the full receptor. Guided by the mechanistic insight on channel activation gained from our targeted MD simulations,[18,31,32] we formulated such a theoretical model.[34]

The model was highly simplified, with minimal ingredients to capture the most essential properties regarding stationary gating. An iGluR was represented by a pair of agonist-bound LBDs linked to a TMD tetramer (Figure A). The closure of the LBD pair was assumed to be synchronous and described by coordinate y; the opening of the ion channel within the TMD tetramer was described by coordinate x; and the extension of the linker was then described by y – x. Exploiting its modular architecture, we assumed that the free energy function, W(x,y), of the iGluR consisted of terms for LBD closure, channel opening, and linker-mediated interdomain coupling, denoted by Wb(y), Wc(x), and W1(y – x), respectively:For illustration, we assumed a parabolic form for Wb(y) and W1(y – x):where kb and kl are the spring constants for LBD closure and linker stretching, and Δ = L0 – Lm, with L0 denoting the linker length when both the LBD cleft and the channel are closed (i.e., y = 0 and x = 0) and Lm is the equilibrium linker length. We assumed a double-well form for Wc(x):which has a deep and shallow minima, with a free energy difference of 4ε/3, at x = −1 and x = 1, i.e., when the channel is closed and open, respectively. The linker transmitted the stabilization provided by LBD closure to the channel-open state, such that the overall free energy function, with appropriate choices for the model parameters (e.g., Δ = 1), had two nearly evenly matched basins (Figure B). One basin, “o”, had its minimum at (x,y) = (1, 0), with the channel open and the LBD cleft closed; and the other basin, “c”, had its minimum at (−0.83, – 0.91), with x and y values corresponding to a closed channel and a semiclosed LBD cleft.

Figure 5

Theoretical model for iGluR gating thermodynamics and kinetics. (A) Conformational coordinates for the LBD (y) and channel (x). The difference y – x defines linker extension. (B) Free energy function of the iGluR, with two basins representing the channel-closed and channel-open states, respectively. (C) Left: value of x as a function of time, from a Brownian dynamics simulation of the model. Right: single-channel current trace of the GluN1/N2A NMDAR during stationery gating. Reproduced with permission from ref (34). Copyright 2015 American Chemical Society.

The stabilization effect on the channel-open state by the linker-mediated coupling to the LBD could be seen by comparing the free energy term Wc(x) for the isolated channel against the potential of mean force, Wpmf(x), in x after averaging out y in the free energy function for the full receptor,The second term on the right-hand side of eq served to stabilize the open channel (i.e., x = 1) relative to the closed channel (i.e., x = −1). Agonist efficacy can be measured by the channel open probability, which in our model is given by the normalized configurational integral for the channel-open state,where x‡ denotes the top of the barrier in the potential of mean force. A reduction in the degree of cleft closure at the LBD free energy minimum corresponds to a decrease in L0 (and hence Δ), whereas a reduction in the curvature of the LBD free energy basin corresponds to a decrease in k1. Both lead to weakened stabilization of the channel-open state provided by agonist-induced LBD closure and hence lower P0, thus providing a theoretical basis for the foregoing mechanistic hypotheses on AMPAR and NMDAR partial agonisms.

We noted that the glycine insertions in the NMDAR M3-D2 linkers[18] increase Lm, and thereby decrease Δ. In essence, these NMDAR insertions have a similar effect as APAR partial agonists. In our study of glycine insertions, we empirically defined a “pulling factor,” as the slope of a linear correlation between two effects of the insertions: ΔΔG, the change in free energy difference between two kinetic substates, and ΔLm, the change in linker length (Figure ). In our theoretical model, we proved that this slope is related to the tension in the linker. Specifically, for the two states c and o in our model,where Fc and F0 denote the linker tensions in the channel-closed and open states, respectively. The theoretical model thus clarified that the single-channel recordings using linker insertions were directly probing changes in linker tensions between kinetic substates along the NMDAR activation pathway.

Figure 6

“Pulling factor” k, empirically defined as the slope of linear correlation between ΔΔG and insertion length, for each transition along the NMDAR activation pathway. (A) Plots correlating ΔΔG (kcal mol–1) and insertion length for GluN2A in the C2 → C1 transition. Also shown is the slope k, in kcal mol–1 nm–1, for linear fits to the first three (black) or four (gray) points. (B) Pulling factors, using three-point analysis, for all the transitions along the activation pathway. C1, C2, C3, O1, and O2 denote kinetic substates along the activation pathway. Reproduced with permission from ref (18). Copyright 2014 Nature Publishing Group.

In addition to thermodynamic properties, our model was extended to stationary gating kinetics, by modeling the motions along x and y as diffusive. The transition rates between c and o were calculated by both a multidimensional reaction rate theory[53] and Brownian dynamics simulations (Figure C).

Possible Conformations of Kinetic Substates Along Activation Pathway

Saturating amounts of agonists can maintain iGluRs in stationary gating, during which the channel switches between closed and open conformations. It is most probable that the LBDs adopt cleft-closed conformations in channel open periods, but very little has been known about LBD conformations in channel closed periods. Kussius and Popescu[16] observed only modest differences in stationary gating properties between agonist-bound and LBD cross-linked NMDARs. Assuming that cross-linking locked the LBD clefts closed, they concluded that LBDs mostly stayed in cleft-closed conformations even during channel closed periods. However, as noted above, our free energy simulations have now cast doubt on their assumption.[35] Our theoretical model,[34] instead, predicts that the LBDs clefts are semiclosed in the channel-closed state. Keeping the LBD clefts closed while the channel is closed would result in overstretched linkers and hence excessive tensions; the receptor relieves such tensions by reducing the degree of LBD closure.

From the GluN1 and GluN2A LBD conformations sampled in our free energy simulations, we looked for those that would be characteristic of LBDs in the channel-closed state.[35]Figure displays cleft-semiclosed conformations of the GluN1 and GluN2A LBD as candidates. Whereas the cleft-closed conformations are stabilized mostly by interactions of the bound agonist and the D1 subdomain with the D2 subdomain, the cleft-semiclosed conformations are stabilized by alternative interactions of the bound agonist and D1 with D2, as well as by interactions of the agonist and the cleft-lining residues with water.

Figure 7

Cleft-semiclosed conformations, possibly characetristic of GluN1 and GluN2A LBDs in the channel-closed state. With the D1 subdomains superimposed, the D2 subdomains exhibit significant dispacements in putative semiclosed LBD conformations (GluN1 marine and GluN2A magenta) relative to the cleft-closed conformations (gray). Inset: distinct intersubdomain hydrogen bonds in cleft-closed and semiclosed LBDs. Conformations from ref (35).

As revealed by single-channel recordings, the channel-closed state of NMDARs contains not a single but three kinetic components outside desensitization.[10−13] It is now time to conjecture the conformations of these kinetic substates along the activation pathway, and combine computation and electrophysiology to delineate and test such conjectures.[54] The first inspiration for our current efforts in modeling the conformations of the substates came from the glycine insertion data (Figure B).[18] They suggest that the earlier transitions (C3 → C2 and C2 → C1) along the activation pathway mostly involve motions within GluN2A whereas the later transitions (C1 → O1 and O1 → O2) involve concerted motions of both types of subunits. More specifically, as clarified by our theoretical model (see eq ), the greatest decreases in M3-D2 linker tensions and hence linker extensions occur in GluN2A during the C3 → C2 and C2 → C1 transitions. Consequently our current conjecture is that the LBDs adopt cleft-semiclosed conformations in all four subunits for C3, but become cleft-closed in one of the two GluN2A subunits for C2, in both GluN2A subunits for C1, and in all the four subunits in the channel-open state. If our current approach proves successful in characterizing conformations of NMDAR kinetic substates, future efforts may target AMPARs, where single-channel electrophysiology is starting to generate hints on (sub)states in channel gating.[8,9]

Toward An Atomic-Level Formalism For Predicting Stationary Gating Properties

Our initial theoretical model[34] has provided a foundation for mechanistic hypotheses on AMPAR and NMDAR partial agonisms, clarified that glycine insertions can be used to probe linker tension, and predicted that the channel-closed state is characterized by semiclosed LBDs. Yet, the most important implication of this model is the exciting prospect for iGluR functional properties to be predicted from intra- and interdomain energetics and dynamics in MD simulations. Our current efforts are directed at developing a more realistic theoretical model specifically for NMDARs, with a free energy function that features local minima corresponding to the kinetic substates. This model can be implemented by all-atom MD simulations to predict thermodynamic and kinetic properties of stationary gating that are amenable to direct test by single-channel recordings.[12,14−19] Our free energy simulations for GluN1 and GluN2A LBDs can be considered as part of this implementation.[33,35]

The resulting atomic-level formalism for modeling stationary gating has the potential to become a powerful technology for NMDAR physiology and pathophysiology. Functional data for disease-associated NMDAR mutations are lagging far behind gene sequencing data.[2] A computational formalism could ultimately complement the traditional electrophysiological approach in filling this widening gap and help realize the promise of precision medicine.