Molecular level activation insights from a NR2A/NR2B agonist

N-methyl D-aspartate receptors (NMDARs), a subclass of glutamate receptors have broad actions in neural transmission for major brain functions. Overactivation of NMDARs leading to “excitotoxicity” is the underlying mechanism of neuronal death in a number of neurological diseases, especially stroke. Much research effort has been directed toward developing pharmacological agents to modulate NMDAR actions for treating neurological diseases, in particular stroke. Here, we report that Alliin, a sulfoxide in fresh garlic, exhibits affinity toward NR2A as well as NR2B receptors based on virtual screening. Biological activities of Alliin on these two receptors were confirmed in electrophysiological studies. Ligand-binding site closure, a structural change precluding ion channel opening, was observed with Alliin during 100 ns molecular dynamics simulation. Alliin interactions with NR2A and NR2B suggest that residues E/A413, H485, T690, and Y730 may play important roles in the conformation shift. Activation of NR2A and NR2B by Alliin can be differentiated from that caused by glutamate, the endogenous neurotransmitter. These characteristic molecular features in NR2A and NR2B activation provide insight into structural requirements for future development of novel drugs with selective interaction with NR2A and NR2B for treating neurological diseases, particularly stroke.


Introduction
Stroke is the third leading cause of death and the leading cause of adult disability in the world. Global stroke death was more than 5 million in 2005 and expected to reach 6.5 million in 2015 and 7.8 million in 2030 based on projections by the World Health Organization. Approximately, 80% of stroke is caused by cerebral ischemia. Therapy for acute ischemic stroke is currently limited to thrombolysis with recombinant tissue plasminogen activator (tPA). tPA has a short therapeutic window of up to 4.5 h from stroke onset and is associated with a 10-fold increase in symptomatic intracranial hemorrhage (The National Institute of Neurological Disorders and Stroke rt-PA Stroke Study Group, 1995;Lee, Zipfel, & Choi, 1999). It is understandable that only approximately 20% of patients with ischemic stroke who arrived within the therapeutic window received tPA. Safer and more effective therapies for acute ischemic stroke are needed. Ischemic brain injury following stroke involves a multitude of pathophysiological processes (Lee et al., 1999). A major mechanism of ischemic neuronal death is excitotoxicity mediated by glutamate receptors, especially of the N-methyl D-aspartate subtype (NMDARs) (Rothman & Olney, 1995;Simon, Swan, Griffiths, & Meldrum, 1984). However, a large number of clinical trials aiming to block various glutamate receptors have failed, mostly due to low receptor subtype selectivity, and adverse side effects inherent to the receptor antagonists (Costa et al., 2010). Virtual screening of natural products offers an efficient method to identify compounds for designing novel derivatives that may be more effective for curtailing the multiple pathophysiological processes leading to ischemic brain injury. Derivatives of natural products to which humans have been exposed to for thousands of years may be less likely to cause serious side effects than the synthetic chemicals.
Recent advances on the mechanisms of excitotoxicity led to the findings that different NMDAR subtypes may have differential roles in the pathogenesis of ischemic neuronal death (Danbolt, 2001;Martin & Wang, 2010). Structurally, NMDARs are heterotetrameric complexes composed of at least one NR1 (GluN1) subunit, coupled to NR2 (GluN2) and NR3 (GluN3) subunits. In the adult forebrain, where stroke normally occurs, the two major subtypes of NMDARs are those containing NR1 and NR2A subunits (so-called the NR2A receptors) and those containing NR1 and NR2B subunits (the NR2B receptors). In particular, the NR2A receptors are thought to be pro-survival with salutary actions in response to ischemic insult. In contrast, the NR2B receptors are the primary cause of neuronal death following cerebral ischemia (Tu et al., 2010;Zhou et al., 2010), with NR2B antagonism the efficacy of NMDAR blockers in experimental studies. Thus, in minimizing adverse reactions and maximizing therapeutic efficacy, NMDAR blockers used for stroke therapy should be designed to be NR2A and NR2B receptor-selective. Development of compounds with selective NR2A agonist and/or NR2B antagonist activity remains a challenging task.
Activation of NMDARs requires simultaneous binding of glycine and glutamate at the ligand-binding domain (LBD) of NR1 and the agonist-binding domain (ABD) of NR2, respectively (Figure 1(A)). Once bound, the cleft closure mechanism of NR1-LBDs and NR2-ABDs will be initiated (Jin, Banke, Mayer, Traynelis, & Gouaux, 2003;Sun et al., 2002;Weston, Schuck, Ghosal, Rosenmund, & Mayer, 2006). During this process, N-terminal NR1-LBD and NR2-ABD will undergo conformational changes and aggregate, leading to the separation of C-terminal NR1-LBD and NR2-ABD and ion channel opening (Figure 1(B)). When blood flow to the brain is interrupted, such as in the case of ischemic stroke, elevated glutamate concentrations will activate NMDARs, leading to Ca 2+ influx and ischemic neuronal death (Lipton & Rosenberg, 1994).
This study applied computer-aided drug design techniques to elucidate molecular-level mechanisms responsible for activating NR2. We utilized traditional Chinese medicine (TCM) Database@Taiwan Tsai, Chang, & Chen, 2011) (http://tcm.cmu.edu. tw/) to rapidly screen for potential agonists, which had been widely used in the identification of novel therapies Chang, Huang, & Chen, 2011;Tou & Chen, 2012;Tou, Chang, Lee, & Chen, 2013;Yang, Chang, Chen, & Chen, 2011). Electrophysiological experiments were conducted to provide biological evidence to screening results. Molecular dynamics simulation was conducted to elucidate the agonistic mechanism of the TCM compounds.

Materials and methods
Homology modeling Due to the lack of an available complete protein structure for NMDAR NR2A, a model for the NR2A protein was constructed using homology modeling. Gene sequence of the target protein NR2A (Q12879) was obtained from Swiss-Port (Magrane & Consortium, 2011), and template proteins, 2A5S and 2A5T, were downloaded from Protein Data Bank (Furukawa, Singh, Mancusso, & Gouaux, 2005). Homology modeling was performed using Protein Modeling package of Discovery Studio 2.5 (DS 2.5; Accelrys Inc., San Diego, CA) (Discovery Studio Modeling Environment v. 2.5, 2009). Ramachandran plot and Profiles-3D in DS 2.5 were used to validate model reliability and compatibility (Kolaskar & Sawant, 1996;Wuthrich, 1990). Ramachandran plot results were plotted using RAMPAGE (Lovell et al., 2002).

Molecular docking and screening
The modeled NR2A was used as the binding site for NR2A. The glutamate binding site for NR2B was downloaded from Protein Data Bank (PDB ID: 2IPV). Glutamate was used as the control for NR2A and NR2B. Virtual screening of compounds from TCM Data-base@Taiwan was conducted under the Chemistry at HARvard Macromolecular Mechanics (CHARMm) Brooks et al., 1983 force field using the LigandFit (Lim-Wilby, Jiang, Waldman, & Venkatachalam, 2002;Venkatachalam, Jiang, Oldfield, Waldman, 2003) module of DS 2.5. Score outputs of Dock Score, ligand internal energy and Jain were selected and results were ranked according to Dock Score. TCM compounds that docked into both NR2A and NR2B were subjected to biological assays to confirm bioactivity.

HEK293 cell culture and transfection
Human embryonic kidney cells (HEK293) culture and transfection procedures were conducted as reported earlier (Liu, Wu, & Wang, 2010). In brief, HEK293 cells were cultured in Dulbecco's modified eagle medium (DMEM; Invitrogen, Carlsbald, CA, USA) supplemented with 10% fetal bovine serum (FBS). Cells were cultured at 37°C in a 5% CO 2 chamber to 60% confluence and transiently transfected using Lipofectamine 2000 (Invitrogen) according to the manufacturer's protocols. Cells were cotransfected with a combination of human NR1 and NR2A or NR1 and NR2B subunit plasmids. The transfection ratio of NR1/NR2A and NR1/NR2B plasmid was 1:1. In addition, pcDNA3-GFP was also co-transfected with NMDA receptor subunits as a transfected marker in order to facilitate the visualization of transfected cells during electrophysiological experiments. In order to achieve optimum plasmid expression, cells were re-plated on glass coverslips 15-20 h following transfection and cultured at 37°C for an additional 15-24hrs prior to whole-cell patch-clamp recordings. Since culture medium contains glutamate, AP550 uM was added to the culture medium after transfection in order to protect cell viability.

Molecular dynamics simulation
To elucidate molecular mechanisms contributing to the observed biological activities, molecular dynamics (MD) simulation was conducted under the force field of CHARMm (Fletcher, 1969) using DS 2.5. Each proteinligand complex was minimized with 1000 cycles each of Steepest Descent (Fletcher & Reeves, 1964) and Conjugate Gradient (Andersen, 1983). SHAKE algorithm (Hammonds & Ryckaert, 1991;Kabsch & Sander, 1983) with time step intervals of .002 ps was used for the duration of MD. Samples were gradually heated from 50 K to 310 K within 50 ps subsequent to minimization, equilibrated for 200 ps and produced using moles, volume and temperature (NVT) canonical ensemble for 100 ns. The MD trajectories of each protein-ligand complex and individual ligands were analyzed using the Analyze Trajectory module in DS 2.5. In addition, the MD conformations were applied to GROMACS for the analysis of secondary structure changes using the Dictionary of Secondary Structure of Proteins (DSSP) (Lin & Song, 2011) program.

Ligand migration map
The LigandPath program was used to compute the entry/ exit paths of each ligand within each protein complex (Fairfield, 1983). The program utilizes Voronoi diagram (Ginalski, 2006;Pinto & Freitas, 2011) to map the free space within the protein at every conformation. By combining individual conformations, possible ligand pathways within the dynamic structural ensemble can be identified. Clearance threshold, which designates the minimum distance between atoms for clearance, was set at 2.5 Å. Probe clearance, which determines surface probe size and is related to protein surface resolution, was set at 4.5 Å. Subsequent to ensemble generation, ligands were redocked to the portals identified by LigandPath for verification.

Results and discussion
Homology modeling for the NR2A subunit Sequence identity and similarity between ABD of the NR2A target sequence Q12879 and template structures 2A5S and 2A5T are 92.1 and 92.4%, respectively, indicating a good model (41). Regions of deviation are not located within the core ABD region ( Figure S1). Majority (96%) of the residues are located within the favoured region in the Ramachandran plot ( Figure S2A). Key ABD residues are located within the acceptable regions of the Ramachandran plot, suggesting that amino acids within the model structure retain a highly possible conformation. Profiles-3D validation of the protein structure also indicates that these key residues have compatibility scores greater than zero ( Figure S2B). The high sequence similarity scores coupled with consistent validation results confirm the robustness of the NR2A model. The generated model was used as the binding site model for subsequent studies. Table S1 lists the top five compounds with the highest Dock Scores in NR2A and NR2B. Compounds with high affinity toward both NR2A and NR2B were dencicihine and Alliin. Dencicihine and Alliin are structurally similar to glutamate, the natural substrate of NR2A and NR2B ABD.

Candidate screening and biological activity verification
To test our computational prediction, we co-overexpressed NR1/NR2A or NR1/NR2B in HEK293 cell line to test for NMDA receptor-mediated current induced by fastperfusion of the TCM candidates using whole cell patchclamp recordings. Due to the lack of commercially available dencicihine, only Alliin was tested. For both NR1/ NR2A and NR1/NR2B receptors in this system, Alliin, along with glycine (20 uM), the NMDA receptor co-agonist targeting the NR1 subunit, elicited a smaller NMDAR-mediated current, compared with that induced by the conventional NMDAR-specific agonist NMDA (50 μM) (Figure 2(A)). We then asked whether Alliin has a preference to act on NR2A-or NR2B-containing NMDARs. Under identical concentrations (100 μM), Alliin induced a larger current in cells expressing NR2A than NR2B (Figure 2(A)). When normalized against the respective NMDA-induced currents in NR2A and NR2B, Alliin induced approximately 5 and 10% of NMDA-activated currents, respectively (Figure 2(B)). These observations provide biological evidence for the agonist activity suggested by computational simulations, and shows that while Alliin may act as a weak agonist for both NR2A and NR2B, it demonstrates higher NR2B receptor efficacy.

Docking analysis
In an attempt to explore the mechanisms underlying differential receptor activities, the structures of NR2A and NR2B were investigated. Sequence alignment shows that the two proteins are highly homologous with a sequence identity of 80.1% ( Figure S3). For sequences relevant to ligand binding, differences were primarily located on Loop 1 (L1).
The apparent differences in sequence affect ligand binding in NR2A and NR2B. Both glutamate (Figure 3 (A)) and Alliin (Figure 3(B)) form hydrogen bonds (H-bond) with E413 at NR2A, but fail to do so with the corresponding residue A413 in NR2B (Figure 3(C) and (D)). The difference is likely due to the hydrophobic nature of alanine compared with glutamic acid. Extension of E413 into proximity to the ligand enables H-bond formation. H-bond with K484 was noted between glutamate and NR2A but not between glutamate and NR2B. Both glutamate and Alliin form pi interactions with Y730 in NR2B but not NR2A. All ligands form long distance pi interactions with H485 regardless of receptor type. Protein-ligand interactions formed by glutamate and Alliin in NR2A or NR2B are tabulated in Table S2.
Intriguingly, a prominent structural difference distant from the ABD is noted from P527 to A538 (Figure 3(E) and (F)). In NR2A, Backbone ξ (P527-A538) and Backbone β (V748-K767) are approximately parallel and form stable interbackbone H-bonds. These stable H-bonds are not observed in NR2B, most likely due to the increased backbone distance caused by the bended structure of Backbone ξ. These backbone differences might be important in NR2A-NR1 or NR2B-NR1 dimerization.

MD simulation
In order for simulation results to be of significance, protein equilibration should be reached prior to analysis. Alliin stabilizes after 10 ns MD regardless of receptor type, but continuous fluctuation is observed for glutamate (Figure 4(A)). Protein complexes, with the exception of NR2A-glutamate complex, reached equilibration during MD. A slight increase is observed     for the NR2A-glutamate complex implying continuous structural changes. Distinctive differences in total energy based on receptor type are detected. NR2A complexes have lower total energy profiles than their NR2B counterparts. In accordance with the complex-RMSD changes, continuous lowering of total energy is noted for the NR2A-glutamate complex. Based on the trajectories, the NR2A-glutamate complex could be continuously undergoing conformational changes to achieve a lower total energy status. H-bond distance trajectories in individual proteinligand complexes during MD are illustrated in Figure 4. H-bond formation between glutamate and NR2A is not stable (Figure 4(B(a))). This may be partially due to continuous conformational changes (Figure 4(A)). NR2A residues forming stable H-bonds with glutamate include Glu413, His485, Ser689, Thr690, and Thr730. In sharp contrast, Alliin forms stable H-bonds with NR2A at Glu413, Lys484, Tyr730, and Asp731 ( Figure 4B(b)). In NR2B, glutamate is capable of maintaining stable, optimal distance H-bonds with six residues (Figure 4B(c)). Distance fluctuations are observed during the initial 35 ns for Alliin, but level out and stabilize with increased simulation time (Figure 4(B(d))). On a larger scale, four distinct loops, I408-I418 (Loop 1, L1 (Table 1), compared with the 55.5% (L1), 93.9% (L2), 8.0% (L3), and 74.53% (L4) for NR2A-Alliin (Table 2). NR2B-Glutamate can form long-term interaction in L1 (92.1%), L2 (97.3%), L3 (99.8%), and L4 (99.0%) (Table 3), whereas NR2B-Alliin does not interact with L1 and only forms interactions with L2 (1.7%), L3 (8.03%), and L4 (98.4%) ( Table 4). A visual summary of the H-bond occupancy results in addition to respective torsion trajectories are shown in Figure 4(C) and (D).
Cross comparison of H-bond results reveals that E413(A413), H485, T690, Y730 may be involved in the activation of NR2A or NR2B (Figure 5(A)). For NR2A, glutamate can interact with E413 and Y730 with higher occurrences than can Alliin, and glutamate can form a H-bond with H485, whereas Alliin cannot. NR2B-glutatmate interactions generates high occurrence H-bonds (>95%) with NR2B-H485 and NR2B-T690. In comparison, Alliin can interact with NR2B-T690 but not NR2B-H485. Differences observed for glutamate and Alliin in their interactions with NR2A and NR2B are shown in Figure 5. In NR2A, glutamate is located at favorable distances to E413, H485, and Y730, enabling the formation of H-bonds ( Figure 5(B)). The relative orientation of Alliin to NR2A-E413, NR2A-H485, and NR2A-Y730 do not favor interactions and thus limiting the formation of H-bonds. Similarly, the orientation of the H485 in NR2B favors interaction with glutamate, but not Alliin ( Figure 5(C)). The difference between glutamate and Alliin binding is significant in terms of agonist activity and may be the distinct features in molecular dynamics that underlie the apparent difference in receptor activities derived from electrophysiological studies. H485 is a key activation residue for NR2A and NR2B where interaction strength determines agonist intensity. Sulfur in Alliin interacts with NR2B-H485 at a distance of 2.863 Å. At this low occurrence, Alliin cannot act as a strong agonist in NR2B, but the effect is stronger than in NR2A in which no interaction is observed. These results correspond to the patch-clamp recordings where Alliin elicited higher proportional current change in NR2B than NR2A (Figure 2(B)). We hypothesize that these four amino acids in four different loops are important for NR2A and NR2B receptor agonist activities.
Analysis of secondary structural changes Secondary structure changes suggest different activation mechanisms for NR2A and NR2B. In NR2A-glutamate, a bend to β-sheet change is observed in Loop 1 and Loop 2 as the cleft begins to close ( Figure S4A). Alliin can form long-term interactions with NR2A-D731 at Loop 4, initiating another structural change involving Loop 3 and Loop 4 ( Figure S4B). Loop 2 in NR2B-glutamate complex changes from a coiled structure to β-sheet and maintains the β-sheet conformation for 100 ns ( Figure S4C). On the other hand, Loop 1 and Loop 2 on the NR2B-Alliin complex only exhibited short-term cooperation ( Figure S4D). Similarly, cooperation of Loop 3 and Loop 4 resulting in the formation of a β-sheet is observed in the NR2Bglutamate complex but not in the NR2B-Alliin complex. Combining biological data from electrophysiology and secondary structural changes, the orchestrated change in Loop 1 and 2 may promote activation. Structural changes induced by Loop 3 and Loop 4 cooperation shows differential functions for NR2A and NR2B. Our results suggest that Loop3/4 cooperation inhibits activation in NR2A but promotes activation in NR2B.

Tertiary structure changes induced by Alliin
Dynamic structural changes between lobe 1 and lobe 2 in NR2A and NR2B during MD is illustrated in Figure 6.  . Illustration of the cleft closure mechanism of agonist-binding domain (ABD) of NR2A and NR2B. The cleft between lobe 1 and lobe 2 of NR2A closes from 18.5 Å to 11.7 Å (glutamate-binding domian) and 12.8 Å (Alliin-binding domian). The cleft of NR2B closed from 27.6 Å to 18.9 Å for the glutamate complex, and to 21.0 Å for the Alliin complex. The cleft of NR2A and NR2B have different sizes but both can stimulate the closure mechanism. Ligand-binding site closure initiated by glutamate and Alliin in NR2A and NR2B observed during 100 ns MD. Surfaces in cyan, purple, magenta, and green represent Loops 1, 2, 3, and 4, respectively. Ligand-binding sites are highlighted by the red circles. Clearance to NR2A-ABD is still observable in NR2A though that of NR2A-Alliin is smaller than that of NR2A-Glu. In NR2B, clearance to the glutamate-binding site remains, but that of Alliin is closed.
For NR2A, cleft opening of the glutamate complex decreases from 18.5 to 11.7 Å, whereas that of Alliin complex is reduced to 12.8 Å. The 27.6 Å cleft opening in NR2B decreases to 18.9 and 21.0 Å in the glutamate and Alliin complexes, respectively. Localized changes in ABD clearance size are also observed between glutamate and Alliin ( Figure 6). ABD clearance to NR2A is approximately .8 Å larger with bound glutamate than Alliin. For NR2B, glutamate-NR2B complex maintained an observable opening to the ABD. Surprisingly, the opening is completely closed in the Alliin-NR2B complex. The closure percentage of glutamate (37%) and Alliin (31%) in NR2A are larger than those observed in NR2B (32% for glutamate and 24% for Alliin). This trend, in addition to ABD clearance changes, match the results observed for electrophysiology shown in Figure 2 (A). The agonist effect on NR2A and NR2B seemed to be positively associated with cleft closure and clearance to the internal ABD.

Ligand migration map
Ligand migration portals were calculated to ascertain the accessibility of the ABD to ligands under dynamic conditions (Figure 7). Respective ligands are redocked at protein surface pathway portals to determine docking affinities. Redocking results are summarized in Table S3. Site 1 is the most direct pathway to the ligand-binding site of either NR2A or NR2B. Since the highest Dock Scores for glutamate and Alliin are generated with Site 1, this indicates that Site 1 is the most possible portal for ligand entry and exit. Detailed descriptions of each entry/exit pathway and respective portal Dock Scores are shown in Table S3.

Conclusion
In this study, we report that a natural substance, Alliin, elicits differential effects on NR2A-and NR2B-containing NMDARs, when compared to the endogenous NMDAR agonist: glutamate. Molecular-level insights suggest the participation of E/A413, H485, T690, Y730 and their corresponding loops in the differential activation of these NMDAR subtypes by glutamate and Alliin. Similarities and differences between glutamate and Alliin in their interactions with NR2A and NR2B receptors and the molecular mechanisms thus derived provide useful information for ligand designs in the development of novel drugs derived from natural products for therapy of ischemic stroke that involves differential NR2A (promoting neuronal survival) and NR2B (facilitating ischemic neuronal death).