Abstract

Genetic variants in the genes GRIN1, GRIN2A, GRIN2B, and GRIN2D, which encode subunits of the N-methyl-D-aspartate receptor (NMDAR), have been associated with severe and heterogeneous neurologic and neurodevelopmental disorders, including early onset epilepsy, developmental and epileptic encephalopathy, intellectual disability, and autism spectrum disorders. Missense variants in these genes can result in gain or loss of the NMDAR function, requiring opposite therapeutic treatments. Computational methods that predict pathogenicity and molecular functional effects of missense variants are therefore crucial for therapeutic applications. We assembled 223 missense variants from patients, 631 control variants from the general population, and 160 missense variants characterized by electrophysiological readouts that show whether they can enhance or reduce the function of the receptor. This includes new functional data from 33 variants reported here, for the first time. By mapping these variants onto the NMDAR protein structures, we found that pathogenic/benign variants and variants that increase/decrease the channel function were distributed unevenly on the protein structure, with spatial proximity to ligands bound to the agonist and antagonist binding sites being a key predictive feature for both variant pathogenicity and molecular functional consequences. Leveraging distances from ligands, we developed two machine-learning based predictors for NMDA variants: a pathogenicity predictor which outperforms currently available predictors and the first molecular function (increase/decrease) predictor. Our findings can have direct application to patient care by improving diagnostic yield for genetic neurodevelopmental disorders and by guiding personalized treatment informed by the knowledge of the molecular disease mechanism.

Introduction

Pathogenic variants in the GRIN family of genes encoding the N-methyl-D-aspartate receptor (NMDAR) subunits have been found in patients with various neuropsychiatric disorders, including autism spectrum disorders, epilepsy, intellectual disability, attention-deficit/hyperactivity disorder, and schizophrenia [1–8]. NMDARs are tetrameric ligand-gated ion channels permeable to Na+, K+, and Ca2+, composed of two glycine-binding GluN1 subunits and two glutamate-binding GluN2 subunits, which can be a combination of any two of GluN2A, GluN2B, GluN2C, or GluN2D [9, 10]. GluN1 subunits are encoded by the gene GRIN1 and GluN2 subunits are encoded by the genes GRIN2A, GRIN2B, GRIN2C, and GRIN2D. Among the pathogenic variants identified in the GRIN gene family, those in GRIN2A (46%) and GRIN2B (38%) account for the vast majority, followed by GRIN1 variants (14%) [11]. Variants in these genes have been associated with a spectrum of neurodevelopmental disorders [12]. For example, GRIN1 and GRIN2B patients can present with mild or severe intellectual disabilities [1, 3, 13]. While some GRIN2A patients have severe intellectual disabilities, roughly half have no intellectual impairment [8, 13, 14]. Most patients with variants in GRIN2A have seizures whereas the majority of patients with variants in GRIN2B do not have seizures [4]. In addition, low muscle tone is rare among GRIN2A patients while common among patients involving other GRIN genes [15, 16]. All GRIN patients present with speech impairment, even those without intellectual disabilities [17]. GRIN2D patients appear to have the most severe phenotype, although there is not yet enough data to understand the full clinical spectrum [15, 16]. More recently, patients with protein-truncating variants in GRIN2A have been associated with schizophrenia [18] and are susceptible to seizures and delayed maturation of parvalbumin interneurons, both of which resolve after adolescence [19].

NMDARs serve many cellular functions, including both pre-synaptically to influence neurotransmitter release and long-term plasticity [20–22] and post-synaptically to mediate the slow component of postsynaptic currents and synaptic plasticity [23, 24]. NMDARs function as signal coincidence detectors, as their activation requires changes in membrane potential to relieve the pore blocking by Mg2+ ions as well as the synaptic release of glutamate [10, 25, 26]. Thus, NMDARs are precisely regulated by the binding of glutamate and its co-agonist glycine, of extracellular Mg2+ that blocks the channel, and of other endogenous extracellular modulators such as zinc ions (Zn2+) [27]. Genetic variants in GRIN genes can cause a heterogeneous spectrum of alterations of the NMDAR function, which can be grouped into two main types: gain of the NMDAR function (or gain-of-function effect or GoF) and partial or complete loss of the NMDAR function (or loss-of-function effect or LoF) [28]. Because of the many functions and modulators of NMDARs, a large number of NMDAR-based molecules have been developed as therapeutic options designed to mitigate dysfunction of the glutamatergic system [29, 30].

Variant interpretation in the GRIN genes, both for pathogenicity and molecular functions, is still challenging. Currently, more than 65% of missense variants in GRIN genes are classified as variant of unknown significance according to the ClinVar database (accessed July 2022) [31]. To improve variant interpretation, several exome-wide bioinformatic approaches have been developed that can identify clusters of patient variants or population variant depletion across genomic, protein sequences or 3D structures [32, 33](p3), [34–37]. These methods and targeted clinical-genetic studies for GRIN2A and GRIN2B showed enrichment of pathogenic over population variants in several UniProt defined domains [13, 14, 37]. Although these approaches can identify important regions, they typically don’t explore the underlying structure-to-function relationship. However, in the absence of functional characterization of every possible genetic variant in NMDAR encoding genes, prediction models are needed to enable precision care since different disease mechanisms have contraindicated treatment needs [38–41]. For example, among patients with developmental and epileptic encephalopathy, those with a variants that cause gain of NMDAR function represent candidates for potential treatment with NMDAR blockers, such as memantine [42–44] or GluN2B-selective inhibitors [41, 45], while those with complete or partial loss of the NMDAR function may potentially respond to positive allosteric modulators of the NMDAR [46–48].

However, the functional consequences are not known for most of the variants. Electrophysiological studies that experimentally determine the molecular functional consequences introduced by missense variants are expensive and time-consuming, and it is difficult to envision how all the possible GRIN missense variants can be functionally assessed. Machine learning (ML)-based methods may be able to take advantage of the limited available experimental data to predict the molecular functional consequences of the NMDAR variants that have not been experimentally tested, as has been proven successful for example in voltage-gated potassium channels [49] and sodium channels [50]. These methods are based on sequence and structural protein features, as discriminative gene-level and protein-level features have been found to be associated with the GoF/LoF effects of variants [51–53]. To date, few predictors are available for variant functional effects (VPatho [54], LoGoFunc https://www.biorxiv.org/content/10.1101/2022.06.08.495288v1.full.pdf) and none are specifically designed for GRIN variants. Here we seek to identify structural features of missense variants that are predictive of variant pathogenicity and of increased or decreased functional effects that are specific for the GRIN genes to develop a ML-based method to predict pathogenicity and Increased/Decreased consequences in GRIN genes. Therefore, we aggregated a unique set of 201 expert-curated patient variants and 631 population variants from the gnomAD database together with 160 functionally characterized missense variants from electrophysiological readouts across GRIN1, GRIN2A and GRIN2B. Some of these variants were functionally characterized in this study for the first time. Since previous work [12, 55–59] described how individual missense variants in NMDARs alter their ligand-induced regulation, we first sought whether spatial distance from agonists that are bound to their binding sites in the NMDAR protein structure correlates with pathogenicity and increased/decreased functional effects. With the identified distances from ligands and additional biophysical and evolutionary scores, we have built a ML-based predictor specifically trained to classify variants in GRIN genes as pathogenic or benign, and a ML-based predictor specifically trained to classify the molecular consequences of missense variants into increased or decreased effect, thus providing a valuable resource for clinical genetics.

Materials and methods

Clinical data set of GRIN missense variants

Clinical dataset

The variant data set comprises manually curated patient missense variants collected through patient registries. Clinical cases were collected using a REDCap survey with 391 fields on genetic and clinical data (REDCap version 10.9.3 https://www.project-redcap.org/) and offer the integration of retrospective and longitudinal data. All clinical cases were manually reviewed and classified in accordance with American College of Medical Genetics (ACMG) guidelines. These variants were found in the genes GRIN1, GRIN2A, and GRIN2B. Missense variants from the general population in the same genes were derived from the gnomAD database [60]. In order to map the variants on the NMDAR structure, we restricted our analysis to missense variants that were localized in the domains of the NMDAR structure that are atomically resolved in the PDB: the amino-terminal (ATD), agonist-binding (ABD), and transmembrane (TMD) domains. Missense variants that are localized in C-terminal domain (CTD) could not be considered as the CTD is not present in any resolved NMDAR structure. This comprises a set of 832 missense variants, of which 201 were from patients and 631 were from controls. The full list of the considered variants along with their clinical annotation is available and interactively accessible at https://GRIN-portal.broadinstitute.org and reported in Supplementary Tables S1S3.

Functional data set of GRIN missense variants

We identified 127 missense variants in GRIN1, GRIN2A and GRIN2B for which functional testing through electrophysiological and biochemical assays were published or completed here to an extent to allow GoF and LoF determination by the criteria of Myers et al. 2023 [28] (Supplementary Table S4). We also identified additional variants in the literature for which some functional data exists, but which lack completion of all assays needed to classify by criteria described in Myers et al. (2023; Supplementary Table S5). Among these two sets of variants, we present new data allowing completion of functional and biochemical assessment of 34 known missense variants according to the criteria of Myers et al. (2023; Supplementary Tables S6, S7). From these determinations, we categorized a total of 160 variants as either having an increased or decreased ion channel function of the response for our implementation here by evaluating the criteria described in Myers et al. (2023) [28], but omitting consideration of surface expression here, since our algorithms explore receptor function as it relates to protein structure, not factors that influence trafficking. Thus, functional data from 160 variants is the starting point for our computational analyses.

Functional analysis of missense variants

cDNAs encoding human NMDAR subunits GluN1-1a (hereafter GluN1; GenBank accession codes: NP_015566), GluN2A (GenBank accession codes: NP_000824), and GluN2B (GenBank accession codes: NM_000825) were used and site-directed mutagenesis used to introduce human variants (QuikChange; Stratagene, La Jolla, CA, United States); all mutant cDNAs were verified by dideoxy DNA sequencing (Eurofins MWG Operon, Huntsville, AL, United States). The cDNA for WT and mutant NMDAR subunits were linearized using the appropriate restriction enzyme and cRNA for each was synthesized in vitro using the mMessage mMachine T7 kit (Ambion, Austin, TX, United States).

Xenopus laevis ovaries with unfertilized oocytes (Stage V-VI) were obtained from Xenopus One Inc (Dexter, MI, United States) and digested with Collagenase Type 4 (Worthington-Biochem, Lakewood, NJ, United States; 850 μg/ml, 15 ml for a half ovary) in Ca2+-free Barth’s solution that contained (in mM) 88 NaCl, 2.4 NaHCO3, 1 KCl, 0.82, maintained at 16°C and injected with cRNA encoding either WT or variant NMDAR subunits (GluN1:GluN2A or GluN2B ratio 1:2, 5–10 ng total in 50 nl of RNAase-free water per oocyte). Injected oocytes were maintained in normal Barth’s solution at 16–19°C.

Two-electrode voltage clamp (TEVC) current recordings from Xenopus oocytes expressing NMDARs were performed as previously described [61, 62]. Oocytes were transferred to a recording chamber and were perfused with extracellular recording solution composed of (in mM) 90 NaCl, 1 KCl, 0.5 BaCl2, 10 HEPES, and 0.01 EDTA (23°C, pH 7.4 with NaOH, EDTA omitted in experiment measuring Mg2+IC50). Current responses to glutamate and glycine were recorded under voltage clamp at a holding potential of −40 mV; current and voltage electrodes were filled with 3 and 0.3 M KCl, respectively. Maximally effective concentrations of agonists (100 μM glutamate and 100 μM glycine) were used unless stated otherwise. The reagent 2-aminoethyl methanethiol sulfonate hydrobromide (MTSEA; Toronto Research Chemicals, Ontario, Canada) was made fresh and used within 30 min.

HEK293 cells (ATCC CRL-1573) were plated on glass coverslips coated with 0.1 mg/ml poly-D-lysine and maintained in Dulbecco’s modified Eagle medium (DMEM) with 10% fetal bovine serum and 10 U/ml streptomycin at 37°C (5% CO2). The cells were transfected with cDNA encoding GluN1, GluN2A, and eGFP at a ratio of 1:1:5, or GluN1, GluN2B, and eGFP at a ratio of 1:1:3 using the calcium phosphate method [61]. 12–24 h post transfection the cells were transferred to the recording chamber and perfused with recording solution composed of (in mM) 150 NaCl, 3 KCl, 1 CaCl2, 10 HEPES, 0.01 EDTA, and 2.10 D-mannitol (the pH was adjusted to 7.4 with NaOH). 3–5 MΩ fire-polished patch electrodes were made from thin-walled micropipettes (TW150F-4, World Precision Instruments, Sarasota, FL, United States) and filled with internal solution composed of (in mM) 110 D-gluconate, 110 CsOH, 30 CsCl, 5 HEPES, 4 NaCl, 0.5 CaCl2, 2 MgCl2, 5 BAPTA, 2 NaATP and 0.3 NaGTP (the pH was adjusted to 7.4 with CsOH; osmolality was 300–305 mOsmol/kg). Whole cell currents in response to application of 1.0 mM glutamate and 100 μM glycine (VHOLD -60 mV, 23°C) were recorded by an Axopatch 200B patch-clamp amplifier (Molecular Devices, Union City, CA, United States). The current responses were filtered at 8 kHz (−3 dB) with an 8-pole Bessel filter (Frequency Devices) and digitized at 20 kHz on a Digidata 1440A system controlled by Clampex 10.3 (Molecular Devices). The position of double-barreled theta-glass tubing was controlled by a piezoelectric translator to obtain rapid solution exchange (Burleigh Instruments, Newton, NJ, United States). Large current responses were corrected off-line for series resistance errors [63].

To measure receptor surface expression, HEK 293 cells grown in 96-well plates (50 000 cells/well) were transfected with cDNA encoding beta-lactamase (β-lac) fused in frame to the N-terminus of GluN1 with WT or mutant GluN2, or similarly constructed β-lac-GluN2 with WT or mutant GluN1 using Fugene6 (Promega, Madison, WI) [55]. Wells treated with Fugene6 alone without cDNA were used to determine background absorbance. NMDAR antagonists (200 μM DL-APV and 200 μM 7-CKA) were added to cultures when transfected. Six wells were transfected for each variant to determine surface and total protein levels (3 wells each). After 24 h, cells were rinsed with Hank’s Balanced Salt Solution (HBSS) that was composed of (in mM) 140 NaCl, 5 KCl, 0.3 Na2HPO4, 0.4 KH2PO4, 6 glucose, 4 NaHCO3 with 10 mM HEPES added. Subsequently, 100 μl of a 100 μM nitrocefin (Millipore, Burlington, MA, United States) solution in HBSS plus HEPES was added to each of the wells and extracellular enzymatic activity was determined. The cells in the three wells were lysed by 50 μl H2O (30 min) prior to the addition of 50 μl of 200 μM nitrocefin to determine total enzymatic activity. The absorbance at 468 nm was read every min for 30 min at 30°C, and the rate of increase in absorbance was determined from the slope of a linear fit to the data.

Structural localization of variants

In order to investigate correlations between functional effects (pathogenicity and molecular function) and structural features of NMDAR missense variants, we localized each variant onto the 3-dimensional (3D) protein structure using SIFTS tools [64] to cross reference amino acid positions between protein sequences and protein structures. We used two crystal structures of the human NMDAR available in the PDB: 7EU7 [65] (3.50 Å of resolution) which comprises two GRIN1-encoded and two GRIN2A-encoded chains, and 7EU8 [65] (4.07 Å of resolution) which comprises two GRIN1-encoded chains and two GRIN2B-encoded chains. We used chain A and B of 7EU7 [65] to localize GRIN1 and GRIN2A variants respectively, and chain B of 7EU8 to localize GRIN2B variants. Of the variants in our collected clinical dataset, 832 could be mapped onto the heterotetrametric protein structure. Of these, 201 are from patients and 631 from population.

Computation of distance features for missense variants

For each variant that could be mapped onto the 3D structure, we computed features that synthesize information about their position in the 3D structure with respect to the functionally important regions of the protein complex. In particular, we computed their distances from: the membrane center, the pore axis, and the four ligands known to primarily regulate NMDAR function (glutamate, glycine, magnesium ion Mg2+, and zinc ion Zn2+). We used the protein complex 7EU7 [65] to compute distances for GRIN1 and GRIN2 variants and 7EU8 [65] for GRIN2B variants. We annotated the membrane through the PPM server (https://opm.phar.umich.edu/ppm_server, version PPM 2.0) [66] and calculated the minimum distance of the wild-type residue of each variant from the membrane center as defined by the PPM server. We annotated the pore using the Mole2.5 webserver [67] (https://mole.upol.cz/, access 24.03.2022) and calculated the minimum distance of the wildtype residue of each variant from the pore axis using R-script and the “bio3D” package [68]. To calculate the minimum distance of the wildtype residue of each variant from the ligands (glutamate, glycine, Mg2+, Zn2+) we first mapped these ligands to the corresponding reference protein structures 7EU7 and 7EU8. Two ligands, glutamate and glycine, were already crystallized in the protein complex 7EU7. To map these two ligands to the 7EU8 protein complex, we performed a structural alignment of the two protein complexes (7EU7 and 7EU8) using mTM-align [69] (root mean square error (RMSD): 3.22 Å). To place Zn2+, we carried out two structural alignments of the amino terminal domain of the NMDAR crystallized with two Zn2+ ions bound to GluN2A (PDB-ID: 5TQ2, resolution of 3.29 Å) with 7EU7 [65] and 7EU8 [65], respectively. The structural alignments were performed with the mTM-align [69] program (RMSD7EU7 = 3.01 Å, RMSD7EU8 = 2.22 Å). Since no NMDAR structure has been crystallized with Mg2+ inside yet, we calculated the estimated the location of the Mg2+ ions by calculating the center of the described magnesium binding residues (N615 GRIN2A, N616 GRIN2B and N616 GRIN1) [2, 70]. For each missense variant, we then calculated the minimum distance of its wild-type residue from each ligand. The distance calculation was restricted to residues located in the same domain. Domain annotations are detailed in Supplementary Table S3 and are: amino terminal domain (ATD), agonist binding domain (ABD), ABD-TMB linkers (S1 and S2), and transmembrane domain (TMD) comprising M1, M2, M3, and M4 helices (Supplementary Table S9). All distances were computed considering all atoms. The considered distance features are tabulated in Supplementary Table S3.

Computation of biophysical and evolutionary features for missense variants

In addition to distance features, we also computed biophysical and evolutionary scores. For each variant, we computed three biophysical features: the relative solvent accessibility of the wild-type residue computed through the DSSP program [71], the difference in hydrophobicity between the wild-type and the substituted residues according to the Kyte-Dolittle scale and the difference in the interaction energy (computed as the Bastolla-Vendruscolo statistical potential [72]) between the wild-type and mutate residue with their structural environment (5 Å sphere centered in the mutated residue). These two latter scores have been weighted through the sequence profile as in Montanucci et al., [73, 74]. The biophysical features are listed in Supplementary Table S1. In addition to biophysical scores, we computed, for each variant in our clinical and functional dataset the BLOSUM62 [75] and EVE [76] evolutionary conservation scores. These scores for each variant are tabulated in Supplementary Table S3.

Enrichment analysis

We performed an enrichment analysis to reveal the distribution of pathogenic and population variants and of variants with increased and decreased functional consequence on the NMDAR structure. To perform the enrichment analysis, we applied the Wilcoxon rank sum test to investigate differences in patient and population variants and their 3D distances from functionally relevant protein region sites (pore axis and membrane) and the four selected ligands. Similarly, we tested differences for variants associated with an increased or decreased functional effect. We applied Bonferroni correction to account for multiple testing.

Developing of ML-based binary classifiers for pathogenicity and functional prediction

Pathogenicity predictor

In order to develop a binary ML-based predictor which classifies missense variants in the GRIN genes as benign or pathogenic, we trained each of four binary predictors, each one on a different set of features, listed in Supplementary Table S2. PP-dist, PP-evo, PP-biophys trained on only distances from ligands and pore axis, evolutionary scores, biophysical features, respectively. Finally, PP-dist&evo is trained on the combination of features of the two best performing predictors, which are distances and evolutionary features. All the predictors have been trained on the subset of the clinical dataset comprising variants that could be mapped on the 3D structure. A 5-fold cross-validation procedure was applied.

Functional predictor

In order to develop a binary ML-based predictor which classifies missense variants in the GRIN genes as increasing or decreasing functional effect, we trained four binary predictors, each one trained on a different feature set, listed in Supplementary Table S2. FP-dist is trained on distances from the four ligands: glutamate, glycine, Mg2+, and Zn2+; FP-evo is trained on evolutionary scores; FP-biophys is trained on the three biophysical features of the variant. Finally, FP-dist&evo is trained on distance and evolutionary features. The predictors have been trained on the functional dataset comprised of 160 variants, of which 89 with increased and 71 with decreased functional effect (see Supplementary Tables S4S8). Details about machine learning implementation, cross-validation and indexes to evaluate performances are described in the Supplementary Notes. In order to compare the performance with existing variant functional effects predictors, we retrieved the prediction of our 160 functionally characterized variants of LoGoFunc. At the moment of submission VPatho was not available.

Results

Spatial proximity to ligands and pore is different in pathogenic versus benign variants

Here we asked whether the spatial distance from important ligands, such as glutamate, glycine, magnesium and Zn2+ and functional important sites such as the pore-axis correlate with variant pathogenicity in the GRIN1, GRIN2A and GRIN2B determined using ACMG criteria. We calculated the 3D distances of the 832 variants that could be mapped on the deposited NMDAR structures from the pore, the membrane and the four ligands that regulate the NMDAR activity and we compared the residue distances from these sites between the two groups of variants, an expert-curated set of patient variants.

Compared to the spatial distribution of population variants, patient variants where closer to the pore (Population variants Median distance = 35 Å, Patient variants Median distance = 18 Å, P = 6.8e-45, Wilcoxon rank sum test, Fig. 1D) and to the closest ligand (Population variants Median distance = 23 Å, Patient variants Median distance = 17 Å, P = 3.7e-15, Wilcoxon rank sum test, Fig. 1E). When performing the same enrichment analysis for each domain separately, we found that in the agonist binding domain, GRIN2A and GRIN2B patient variants in the agonist binding domain are located closer to glutamate, which is bound in the cleft of the bilobed agonist binding domain, compared to population variants (PGRIN2A = 4.3e-14, and PGRIN2B = 4.1e-03, Supplementary Fig. S1), while no significant difference in proximity to glycine was found for variants in GRIN1. In the amino terminal domain, no significant difference is found in proximity to ligands between patient variants. However, this result could be due to the small number of population variants in this domain. In the transmembrane binding site, we observed that patient variants are closer to the Mg2+ binding site compared to population variants in GRIN2A (P = 2.6e-02, Supplementary Fig. S1).

Variant distance from the pore axis and from ligands significantly differs between patient and population variants in the GRIN-genes. (A) GluN1-GluN2A heterotetramer protein complex (PDB-ID: 7eu7). The four ligand binding sites are highlighted. (B) NMDAR structure with the 631 population variants of this study represented as spheres. (C) NMDAR structure with the 201 patient variants of this study represented as spheres. (D) Boxplot of the distance from the pore axis for 201 patients and 631 population missense variants in GRIN1, GRIN2A, and GRIN2B. (E) Boxplot of the minimum distance from the closest ligand (glutamate, glycine, Zn2+and Mg2+) for 201 patients and 631 population missense variants in GRIN1, GRIN2A, and GRIN2B. To quantify the differences in the distances to the ligands and the pore axis we performed the Wilcoxon-rank sum test and corrected for eight tests using Bonferroni correction.
Figure 1

Variant distance from the pore axis and from ligands significantly differs between patient and population variants in the GRIN-genes. (A) GluN1-GluN2A heterotetramer protein complex (PDB-ID: 7eu7). The four ligand binding sites are highlighted. (B) NMDAR structure with the 631 population variants of this study represented as spheres. (C) NMDAR structure with the 201 patient variants of this study represented as spheres. (D) Boxplot of the distance from the pore axis for 201 patients and 631 population missense variants in GRIN1, GRIN2A, and GRIN2B. (E) Boxplot of the minimum distance from the closest ligand (glutamate, glycine, Zn2+and Mg2+) for 201 patients and 631 population missense variants in GRIN1, GRIN2A, and GRIN2B. To quantify the differences in the distances to the ligands and the pore axis we performed the Wilcoxon-rank sum test and corrected for eight tests using Bonferroni correction.

Machine learning method based on ligand and pore proximity to predict variant pathogenicity

Given the strong variant position to pathogenicity associations that we observed, we explored whether our generated distance features could be used to develop a method for the prediction of variant pathogenicity in the GRIN genes using our collected clinical dataset of 201 patient and 631 population variants. We developed a ML-based method, PP-dist, based on only distance features (distance from pore axis, glutamate, glycine, Mg2+ion and Zn2+) to classify a GRIN variant as pathogenic or benign. As a comparison, we trained two additional pathogenicity predictors on the same clinical dataset, PP-evo trained with only evolutionary features, and PP-biophys trained on biophysical properties of the amino-acid substitution (see Supplementary Tables S1 and S2 for features descriptions and input features of each predictor). The binary classifier based on only distances from ligands, PP-dist, reaches high prediction performances, with an overall accuracy of 0.892, an area under the ROC curve (AUC) of 0.9237 and a Matthews correlation coefficient (MCC) of 0.698 (see Table 1 and Fig. 2A). The predictive power of these distances is therefore very high, considering that these input features are based on only the position of the substituted residue and do not contain any information on the properties of the alternative residue. This indicates that the distance from pore and ligands, are major features in determining pathogenicity of variants in the GRIN genes. While PP-biophys shows a low MCC of 0.156 indicative of a poor predictor, PP-evo shows a high MCC (0.534) and an overall accuracy of 0.832 (Table 1 and Fig. 2A). This suggests that biophysical features provide much less information about variant pathogenicity than evolutionary scores. When a last predictor PP-dist&evo was trained with features from the two classes of distances from pore and ligands and evolutionary, the prediction performances improved by 1% in respect to the individual predictors PP-dist and PP-evo, reaching an overall accuracy of 0.903, a MCC of 0.726 and an AUC of 0.945.

Table 1

Prediction performances of the pathogenicity and functional predictors.

Pathogenicity predictor (on the 838 clinically characterized variants)
PredictoraQ2MCCAUCTNRNPVTPRPPVNPatientNPopulation
PP-dist0.8920.6980.9370.9420.9180.7360.801159498
PP-evo0.8320.5340.8730.9250.8950.5650.724200572
PP-biophys0.7480.1560.6580.9290.7800.1790.444201631
PP-dist&evo0.9030.7260.9450.9540.9200.7410.837201631
AlphaMissense0.6850.4820.9100.5960.9820.9650.432201631
EVE0.8250.5500.8840.8760.8870.6800.657200572
REVEL0.6860.2870.7600.7090.8500.6130.405199614
Functional (Increased/Decreased function) predictor (on the 160 functionally characterized variants)
PredictoraQ2MCCAUCTNRNPVTPRPPVNincreaseNdecrease
FP-dist0.7400.4820.8000.6440.7760.8280.7166459
FP-evo0. 68410.2690.6350.5510.6130.7140.6598469
FP-biophys0.7030.4040.7400.5000.7610.8710.6798570
FP-dist&evo0.7650.5230.8090.6670.7800.8450.7558469
LoGoFunc0.5220.2170.5100.9710.4820.1610.8757087
Pathogenicity predictor (on the 838 clinically characterized variants)
PredictoraQ2MCCAUCTNRNPVTPRPPVNPatientNPopulation
PP-dist0.8920.6980.9370.9420.9180.7360.801159498
PP-evo0.8320.5340.8730.9250.8950.5650.724200572
PP-biophys0.7480.1560.6580.9290.7800.1790.444201631
PP-dist&evo0.9030.7260.9450.9540.9200.7410.837201631
AlphaMissense0.6850.4820.9100.5960.9820.9650.432201631
EVE0.8250.5500.8840.8760.8870.6800.657200572
REVEL0.6860.2870.7600.7090.8500.6130.405199614
Functional (Increased/Decreased function) predictor (on the 160 functionally characterized variants)
PredictoraQ2MCCAUCTNRNPVTPRPPVNincreaseNdecrease
FP-dist0.7400.4820.8000.6440.7760.8280.7166459
FP-evo0. 68410.2690.6350.5510.6130.7140.6598469
FP-biophys0.7030.4040.7400.5000.7610.8710.6798570
FP-dist&evo0.7650.5230.8090.6670.7800.8450.7558469
LoGoFunc0.5220.2170.5100.9710.4820.1610.8757087

aThe performances of all the predictors are computed with a threshold of 0.5. The Q2 is the overall accuracy of the predictor; MCC is the Matthews correlation coefficient; AUC is the area under the ROC curve; TNR is the true negative rate; NPV is the negative predicted value; TPR is the true positive rate; PPV is the positive predicted value (See Supplementary material for the corresponding equations). NPatient and NPopulation and Nincrease and Ndecrease are the number of variants found in patients and in the general population, and in the training dataset. The best predictor is highlighted in bold.

Table 1

Prediction performances of the pathogenicity and functional predictors.

Pathogenicity predictor (on the 838 clinically characterized variants)
PredictoraQ2MCCAUCTNRNPVTPRPPVNPatientNPopulation
PP-dist0.8920.6980.9370.9420.9180.7360.801159498
PP-evo0.8320.5340.8730.9250.8950.5650.724200572
PP-biophys0.7480.1560.6580.9290.7800.1790.444201631
PP-dist&evo0.9030.7260.9450.9540.9200.7410.837201631
AlphaMissense0.6850.4820.9100.5960.9820.9650.432201631
EVE0.8250.5500.8840.8760.8870.6800.657200572
REVEL0.6860.2870.7600.7090.8500.6130.405199614
Functional (Increased/Decreased function) predictor (on the 160 functionally characterized variants)
PredictoraQ2MCCAUCTNRNPVTPRPPVNincreaseNdecrease
FP-dist0.7400.4820.8000.6440.7760.8280.7166459
FP-evo0. 68410.2690.6350.5510.6130.7140.6598469
FP-biophys0.7030.4040.7400.5000.7610.8710.6798570
FP-dist&evo0.7650.5230.8090.6670.7800.8450.7558469
LoGoFunc0.5220.2170.5100.9710.4820.1610.8757087
Pathogenicity predictor (on the 838 clinically characterized variants)
PredictoraQ2MCCAUCTNRNPVTPRPPVNPatientNPopulation
PP-dist0.8920.6980.9370.9420.9180.7360.801159498
PP-evo0.8320.5340.8730.9250.8950.5650.724200572
PP-biophys0.7480.1560.6580.9290.7800.1790.444201631
PP-dist&evo0.9030.7260.9450.9540.9200.7410.837201631
AlphaMissense0.6850.4820.9100.5960.9820.9650.432201631
EVE0.8250.5500.8840.8760.8870.6800.657200572
REVEL0.6860.2870.7600.7090.8500.6130.405199614
Functional (Increased/Decreased function) predictor (on the 160 functionally characterized variants)
PredictoraQ2MCCAUCTNRNPVTPRPPVNincreaseNdecrease
FP-dist0.7400.4820.8000.6440.7760.8280.7166459
FP-evo0. 68410.2690.6350.5510.6130.7140.6598469
FP-biophys0.7030.4040.7400.5000.7610.8710.6798570
FP-dist&evo0.7650.5230.8090.6670.7800.8450.7558469
LoGoFunc0.5220.2170.5100.9710.4820.1610.8757087

aThe performances of all the predictors are computed with a threshold of 0.5. The Q2 is the overall accuracy of the predictor; MCC is the Matthews correlation coefficient; AUC is the area under the ROC curve; TNR is the true negative rate; NPV is the negative predicted value; TPR is the true positive rate; PPV is the positive predicted value (See Supplementary material for the corresponding equations). NPatient and NPopulation and Nincrease and Ndecrease are the number of variants found in patients and in the general population, and in the training dataset. The best predictor is highlighted in bold.

Predictors performances (A) pathogenicity predictors ROC curve for the ML-based pathogenicity predictors (PP-dist, PP-evo, PP-biophys and PP-dist&evo) and three additional pathogenicity predictors and scores (AlphaMissense [77], EVE [78], and REVEL [79]) on the dataset composed of the 832 GRIN variants of the clinical dataset. (B) Functional predictors ROC curve for the four ML-based functional predictors (PP-dist, PP-evo, PP-biophys, PP-dist&evo), the EVE score and AlphaMissense pathogenicity predictor, and the LoGoFunc LoF/GoF predictor on the dataset composed of 160 variants in the GRIN genes.
Figure 2

Predictors performances (A) pathogenicity predictors ROC curve for the ML-based pathogenicity predictors (PP-dist, PP-evo, PP-biophys and PP-dist&evo) and three additional pathogenicity predictors and scores (AlphaMissense [77], EVE [78], and REVEL [79]) on the dataset composed of the 832 GRIN variants of the clinical dataset. (B) Functional predictors ROC curve for the four ML-based functional predictors (PP-dist, PP-evo, PP-biophys, PP-dist&evo), the EVE score and AlphaMissense pathogenicity predictor, and the LoGoFunc LoF/GoF predictor on the dataset composed of 160 variants in the GRIN genes.

To compare our developed predictors with existing methods, we also report the performances of other three among the main pathogenicity scores and predictors on the clinical dataset of this study (see Table 1 and Fig. 2A). This comparison shows that the simple distance measurements from ligands and pore allow a prediction accuracy that is ~ 10% higher than the best non-gen-specific pathogenicity predictor.

Validation of the pathogenicity predictor and pathogenicity prediction for variants of uncertain significance

So far, we only used expert curated patient variants. As a further validation, we applied our best pathogenicity predictor, PP-dist&evo, to classify variants that were in ClinVar (July 2022) and, at the same time, were not part of our expert curated dataset (for prediction scores see Supplementary Table S3). The total number of ClinVar variants that were not in our clinical dataset and were classified as (likely-) benign and (likely-) pathogenic and located in our 75% most confident class assignments is 44. Of these variants, 39 were correctly classified by our method, reaching a prediction accuracy of 0.89 (Fig. 3A). We then used our best model, PP-dist&evo, to classify 100 ClinVar variants that were of uncertain or conflicting significance (VUS). Predictions for VUS were distributed across the whole spectrum of the pathogenicity score (Fig. 3B). We reclassified 95 VUS that were assigned with a prediction score within the 75% most confident class assignments. Out of these 95 VUS, we predicted 19% (n = 18 VUS) as pathogenic and 81% (n = 77 VUS) as benign.

Pathogenicity prediction for ClinVar variants that were not part of our clinical dataset (A) pathogenicity prediction for ClinVar variants that were not comprised in our clinical dataset according to PP-dist&evo. Prediction scores > 0.5 indicate pathogenic variants, while predictor scores < 0.5 indicate benign variants. The small grey line above and below 0.5 corresponds to the 75% most confident class assignments (over all possible amino acid variants). Prediction scores are displayed along the multiple sequence alignment of GRIN1, GRIN2A & GRIN2B. Blue, red, and grey dots correspond to ClinVar benign, pathogenic and VUS variants, respectively. (B) Distribution of PP-dist&evo prediction for ClinVar variants classified, in ClinVar, as benign, pathogenic and variants of unknown significance (VUS).
Figure 3

Pathogenicity prediction for ClinVar variants that were not part of our clinical dataset (A) pathogenicity prediction for ClinVar variants that were not comprised in our clinical dataset according to PP-dist&evo. Prediction scores > 0.5 indicate pathogenic variants, while predictor scores < 0.5 indicate benign variants. The small grey line above and below 0.5 corresponds to the 75% most confident class assignments (over all possible amino acid variants). Prediction scores are displayed along the multiple sequence alignment of GRIN1, GRIN2A & GRIN2B. Blue, red, and grey dots correspond to ClinVar benign, pathogenic and VUS variants, respectively. (B) Distribution of PP-dist&evo prediction for ClinVar variants classified, in ClinVar, as benign, pathogenic and variants of unknown significance (VUS).

Spatial proximity to ligands and pore is different in variants which cause an increased versus decreased functional consequence

After exploring the relationship between residue localization on the 3D-structure and pathogenicity of variants in the GRIN genes, we explored the relationship between the variant localization on the 3D-structure (Supplementary Table S3) and the differential functional outcomes of the GRIN variants. We used data from 160 variants combining data from the peer reviewed literature and 33 variants generated in this study (Supplementary Tables S4S8). In Fig. 4 the boxplot of the distance of variants with different molecular effect is shown. Variants in GRIN2A and GRIN2B with decreasing effect are significantly closest (P = 7.2e-3) to the glutamate in respect to variants with an increasing effect. Conversely, variants with increasing effect are significantly closest (P = 4.3e-2) to Mg2+ in respect to variants with a decreasing effect. We evaluated function based on potency measurements for glutamate, glycine, Mg2+, and Zn2+. Because there is no reason a priori to assume that distance to the pore or agonist binding site will necessarily correlate with trafficking, which involves other parts of the receptor, we did not include an assessment of surface expression. This means that our categorization as Increasing and Decreasing function will not necessarily predict patient Gain- or Loss-of-Function as defined by Myers et al. (2023) [28]. Rather, this predictor will indicate whether a missense variant is likely to alter function of a receptor once it reaches the cell surface. In addition, for a subset of 69 variants, we did not have a measure of variant actions on deactivation time course, however this correlates with agonist EC50 and thus this effect is in part captured by our measure of potency in these variants [80–82]. We aggregated electrophysiological readouts from published and newly recorded data (Supplementary Tables S4S7) in this study for 160 variants and classified all variants according to their molecular effect as described in the methods. We found the number of variants could be categorized as follows: NIncrease = 89, NDecrease = 71. The results for this set of variants are summarized in Supplementary Tables S2S5. Distances from ligands are reported in Supplementary Table S3. We observed that proximity to the pore, and in particular to the Mg2+ binding sites are associated with variants that “increased” function effects and depleted for variants that “decreased” function in all the GRIN genes (Fig. 4E). In GRIN2A (P-value = 0.01), variants that “decreased” function are associated with close proximity to glutamate binding sites, compared to variants that “increased” function (Fig. 4C).

Variant distances from NMDAR ligands significantly differ between variants with different electrophysiological readouts in the GRIN-genes. In plot (A and B) variants whose molecular effect was classified into either decreasing (A) or increasing (B) effect are shown on the NMDA structure. In plots (C–E) the boxplot of the distances between the variant and the ligands for the variants with different molecular effect are shown. Only variants that are located in the domain where the target ligand binds are considered. Distances from C) glutamate, D) glycine, and E) Mg2+ are shown.
Figure 4

Variant distances from NMDAR ligands significantly differ between variants with different electrophysiological readouts in the GRIN-genes. In plot (A and B) variants whose molecular effect was classified into either decreasing (A) or increasing (B) effect are shown on the NMDA structure. In plots (C–E) the boxplot of the distances between the variant and the ligands for the variants with different molecular effect are shown. Only variants that are located in the domain where the target ligand binds are considered. Distances from C) glutamate, D) glycine, and E) Mg2+ are shown.

ML-based method using ligand and pore proximity to predict functional effect of variants

With a similar procedure as that for pathogenicity predictor, we independently developed four different binary classifiers to predict the functional consequences of variants in GRIN1, GRIN2A and GRIN2B described in Supplementary Table S2. All performances are shown in Table 1 and Fig. 2. While the predictors based on only evolutionary or biophysical features (FP-evo and FP-biophys, respectively) are poor predictors (with MCC of 0.269 and 0.404, respectively), the FP-dist predictor, based only on the 3D distances of variants from each of the four considered ligands, reaches an overall accuracy of 0.740 and a MCC of 0.482. This shows here, for the first time, that distances of the variant from ligands contain information about the direction (“increasing” or “decreasing” activity) of the functional effect of the variant and that can be used for functional effect prediction.

If we combine distance and evolutionary scores, our final best functional predictor for GRIN variants, FP-dist&evo reaches an overall accuracy of 0.765 and a MCC of 0.523 (for predictions see Supplementary Table S3). As a comparison, we also reported in Table 1 the performances of the pathogenicity predictor Alphamissense and of the only other available non-protein specific functional (LoF/GoF) predictor LoGoFunc (https://www.biorxiv.org/content/10.1101/2022.06.08.495288v1.full.pdf). These performances show that pathogenicity prediction scores like Alphamissense have little to no predictive power for the functional effects of variants and FP-dist&evo clearly outperforms, on the GRIN variants, the only other available functional predictor.

Discussion

Tailoring treatment to individual patients’ genetic variants has made significant progress in many fields of medicine in recent years [83]. For disorders caused by genetic variants in GRIN genes, the possibility of successful drug treatment critically depends on the knowledge of the change in function caused by the pathogenic variant, that is gain or loss (partial or complete) of the NMDAR function. Indeed, the knowledge of the molecular mechanism affected by a variant can guide safe and effective targeted treatment [44]. To date, only 20% of all the single amino acid exchanges in GRIN genes have been characterized by electrophysiological and biochemical readouts [27]. In this work we collected the most comprehensive dataset of GRIN variants located in the structurally resolved regions of the NMDAR, with 201 patient variants curated by clinical and genetics experts together and 631 population variants and 160 variants whose molecular consequences have been characterized by in vitro electrophysiology. Taking advantage of this comprehensive clinical and functional dataset of GRIN variants and using machine learning, we identify distance from ligands as a main predictive feature for both pathogenicity and functional prediction and we provide the most accurate pathogenicity prediction specifically developed for GRIN variants to date, as well as the first computational method for the prediction of variant functional consequences in the GRIN genes (for prediction scores see Supplementary Table S3).

On the one hand, protein-unspecific pathogenicity predictors take the advantage of large training datasets and of capturing general principles of protein functioning though the use of massive evolutionary information [78, 84–86]. On the other hand, protein-specific pathogenicity predictors that can incorporate knowledge on protein-specific structure and function, have been shown to enhance the accuracy of pathogenicity prediction [49, 50, 87, 88] (p1). When allowed by data availability, protein-specific pathogenicity predictors can improve the accuracy of pathogenicity prediction, by capturing protein-specific structural patterns and constraints. Here we increase the accuracy of pathogenicity prediction for missense variants in the GRIN genes by integrating structural information on the distance of variant residues to functionally important sites, ligands and the pore.

We also showed that distance from ligands and the pore has predictive power also for functional prediction. A limitation of this work, however, is that we explicitly omit effects of variants on trafficking, which almost certainly have structural determinants beyond the ion channel pore and agonist binding pocket. Surface expression is clinically important, and a topic we will explore separately. So, our functional predictor is predicting increased and decreased function and not GoF/LoF. In line with previously reported observations in which variants with a LoF effect are predominantly located in the ABD domain, due to a disturbance of the agonist binding sites of GluN1 and GluN2A/B [4, 14] we also find an enrichment of variants with increased functional effect. For a couple of variants close to the ligand binding site it has been proposed that the amino acid substitutions lead to a reduced agonist binding [58, 59]. In contrast most variants in the TMD have been shown to have a GoF (or complex) effect [82, 89, 90]. Here we quantify for the first time the correlation of the functional effects of variants in NMDAR proteins with their spatial distance from the ligands. Hence these 3D distances can be used as a proxy to estimate the functional consequence of yet untested variants. Although this is the first step towards an accurate model specifically designed to predict the functional effect of NMDAR variants, we show that protein-unspecific models trained for pathogenicity were not sufficient to develop a strong prediction model for variant effects in NMDAR genes, and only including our newly generated distance features sufficiently boosts the performances to allow functional prediction.

This work is affected by the limitations that classification of the functional consequence of a variant in the two classes of LoF and GoF is an over-simplification of the real biophysical modifications which take place on a molecular level. Still, we could separately show that the distance of variant residues to ligands that regulate the NMDAR in particular correlate with the ligand-specific fold change potency, demonstrating its validity also on the level of individual electrophysiological readouts. Consequently, once more specific electrophysiological readouts will become available, our spatial distance to functional site annotation will become even more helpful to train models that predict the molecular consequence with higher granularity.

In summary, we introduced a potentially powerful approach to predict the directionality of the functional effects of likely pathogenic missense variants in GRIN genes. In a clinical setting like, treatment decisions must often be made before functional studies of disease-causing variants can be performed. In the future, our prediction method could be adapted and benchmarked for use in conjunction with best current clinical practices, for example, to predict which individuals with pathogenic variants may be likely to benefit from a particular treatment based on their variants’ LoF or GoF effects. Our method could potentially be refined with large-scale experimental data, for example, by introducing more specific types of predictions than the mere binary LoF and GoF classification, such as directly predicting a specific change of potency (e.g. glutamate). Because most GRIN genes are depleted for functional variants in the general population, it is likely that more GRIN genes could contribute to disease for which disease associations or mechanisms have not yet been elucidated and to which our method could potentially be applied—such as LoF variants in GRIN2D. In future iterations, also clinical and phenotypic data might be incorporated to enhance predictions for the underlying molecular defect.

Acknowledgements

We thank Rui Song, Wenshu XiangWei, Yuchen Xu, Daniel Teuscher, Lingling Xie, Ruth K. Mizu, and Wenjuan Chen for valuable assistance and sharing unpublished data. This work was supported by a grant from SFARI (732132 to SFT, TAB, JL). This work was supported by NS111619 (ST), HD082373 (HY), and AG072142 (SJM). Funding for this work was provided from the German Federal Ministry for Education and Research (BMBF) to DL, TB and PM (Treat-ION, 01GM1907D) and PM (Treat-Ion2, 01GM2210B) and the Fonds National de la Recherche Luxembourg (Research Unit FOR-2715, FNR grant INTER/DFG/21/16394868 MechEPI2) to PM.

Conflict of interest statement

H.Y. and S.F.T. are co-inventors of Emory-owned intellectual property. S.F.T. is a member of the SAB for Sage Therapeutics, Eumentis Therapeutics, Neurocrine, the GRIN2B Foundation, the CureGRIN Foundation, and CombinedBrain. S.F.T. is a consultant for GRIN Therapeutics. H.Y. is the PI on a research grant from Sage Therapeutics and GRIN Therapeutics to Emory. S.F.T. is PI on a research grant from GRIN Therapeutics to Emory. S.F.T. is cofounder of NeurOp, Inc. and Agrithera. DL received funds from the Simons Foundation and GRIN Therapeutics.

Funding

None declared.

References

1.

Lemke
 
JR
,
Geider
 
K
,
Helbig
 
KL
. et al.  
Delineating the GRIN1 phenotypic spectrum: a distinct genetic NMDA receptor encephalopathy
.
Neurology
 
2016
;
86
:
2171
2178
.

2.

Endele
 
S
,
Rosenberger
 
G
,
Geider
 
K
. et al.  
Mutations in GRIN2A and GRIN2B encoding regulatory subunits of NMDA receptors cause variable neurodevelopmental phenotypes
.
Nat Genet
 
2010
;
42
:
1021
1026
.

3.

Hamdan
 
FF
,
Gauthier
 
J
,
Araki
 
Y
. et al.  
Excess of de novo deleterious mutations in genes associated with glutamatergic systems in nonsyndromic intellectual disability
.
Am J Hum Genet
 
2011
;
88
:
306
316
.

4.

Myers
 
SJ
,
Yuan
 
H
,
Kang
 
JQ
. et al.  
Distinct roles of GRIN2A and GRIN2B variants in neurological conditions
.
F1000Res
 
2019
;
8
:
F1000
.

5.

Carvill
 
GL
,
Regan
 
BM
,
Yendle
 
SC
. et al.  
GRIN2A mutations cause epilepsy-aphasia spectrum disorders
.
Nat Genet
 
2013
;
45
:
1073
1076
.

6.

Lemke
 
JR
,
Lal
 
D
,
Reinthaler
 
EM
. et al.  
Mutations in GRIN2A cause idiopathic focal epilepsy with rolandic spikes
.
Nat Genet
 
2013
;
45
:
1067
1072
.

7.

Lemke
 
JR
,
Hendrickx
 
R
,
Geider
 
K
. et al.  
GRIN2B mutations in west syndrome and intellectual disability with focal epilepsy
.
Ann Neurol
 
2014
;
75
:
147
154
.

8.

Reutlinger
 
C
,
Helbig
 
I
,
Gawelczyk
 
B
. et al.  
Deletions in 16p13 including GRIN2A in patients with intellectual disability, various dysmorphic features, and seizure disorders of the rolandic region
.
Epilepsia
 
2010
;
51
:
1870
1873
.

9.

Laube
 
B
,
Kuhse
 
J
,
Betz
 
H
.
Evidence for a tetrameric structure of recombinant NMDA receptors
.
J Neurosci
 
1998
;
18
:
2954
2961
.

10.

Paoletti
 
P
,
Bellone
 
C
,
Zhou
 
Q
.
NMDA receptor subunit diversity: impact on receptor properties, synaptic plasticity and disease
.
Nat Rev Neurosci
 
2013
;
14
:
383
400
.

11.

XiangWei
 
W
,
Jiang
 
Y
,
Yuan
 
H
.
De novo mutations and rare variants occurring in NMDA receptors
.
Curr Opin Physiol
 
2018
;
2
:
27
35
.

12.

Benke
 
TA
,
Park
 
K
,
Krey
 
I
. et al.  
Clinical and therapeutic significance of genetic variation in the GRIN gene family encoding NMDARs
.
Neuropharmacology
 
2021
;
199
:
108805
.

13.

Platzer
 
K
,
Yuan
 
H
,
Schütz
 
H
. et al.  
GRIN2B encephalopathy: novel findings on phenotype, variant clustering, functional consequences and treatment aspects
.
J Med Genet
 
2017
;
54
:
460
470
.

14.

Strehlow
 
V
,
Heyne
 
HO
,
Vlaskamp
 
DRM
. et al.  
GRIN2A-related disorders: genotype and functional consequence predict phenotype
.
Brain
 
2019
;
142
:
80
92
.

15.

Li
 
D
,
Yuan
 
H
,
Ortiz-Gonzalez
 
XR
. et al.  
GRIN2D recurrent De novo dominant mutation causes a severe epileptic encephalopathy treatable with NMDA Receptor Channel blockers
.
Am J Hum Genet
 
2016
;
99
:
802
816
.

16.

Platzer
 
K
,
Krey
 
I
,
Lemke
 
JR
. GRIN2D-related developmental and epileptic encephalopathy. In:
Adam
 
M.P.
,
Everman
 
D.B.
,
Mirzaa
 
G.M.
. et al.
(eds.),
GeneReviews®
.
Seattle
:
University of Washington
,
2022
, .

17.

Myers
 
KA
,
Scheffer
 
IE
. GRIN2A-related speech disorders and epilepsy. In:
Adam
 
M.P.
,
Everman
 
D.B.
,
Mirzaa
 
G.M.
. et al.
(eds.),
GeneReviews®
.
Seattle
:
University of Washington
,
1993
, .

18.

Liu
 
R
,
Dang
 
W
,
Du
 
Y
. et al.  
Correlation of functional GRIN2A gene promoter polymorphisms with schizophrenia and serum D-serine levels
.
Gene
 
2015
;
568
:
25
30
.

19.

Camp
 
CR
,
Vlachos
 
A
,
Klöckner
 
C
. et al.  
Loss of Grin2a causes a transient delay in the electrophysiological maturation of hippocampal parvalbumin interneurons
.
Commun Biol
 
2023
;
6
:
952
.

20.

Wang
 
H
,
Lv
 
S
,
Stroebel
 
D
. et al.  
Gating mechanism and a modulatory niche of human GluN1-GluN2A NMDA receptors
.
Neuron
 
2021
;
109
:
2443
2456.e5
.

21.

Bouvier
 
G
,
Larsen
 
RS
,
Rodríguez-Moreno
 
A
. et al.  
Towards resolving the presynaptic NMDA receptor debate
.
Curr Opin Neurobiol
 
2018
;
51
:
1
7
.

22.

Dore
 
K
,
Stein
 
IS
,
Brock
 
JA
. et al.  
Unconventional NMDA receptor Signaling
.
J Neurosci
 
2017
;
37
:
10800
10807
.

23.

Iacobucci
 
GJ
,
Popescu
 
GK
.
NMDA receptors: linking physiological output to biophysical operation
.
Nat Rev Neurosci
 
2017
;
18
:
236
249
.

24.

Hunt
 
DL
,
Castillo
 
PE
.
Synaptic plasticity of NMDA receptors: mechanisms and functional implications
.
Curr Opin Neurobiol
 
2012
;
22
:
496
508
.

25.

Attwell
 
D
,
Gibb
 
A
.
Neuroenergetics and the kinetic design of excitatory synapses
.
Nat Rev Neurosci
 
2005
;
6
:
841
849
.

26.

Traynelis
 
SF
,
Wollmuth
 
LP
,
McBain
 
CJ
. et al.  
Glutamate receptor ion channels: structure, regulation, and function
.
Pharmacol Rev
 
2010
;
62
:
405
496
.

27.

Hansen
 
KB
,
Wollmuth
 
LP
,
Bowie
 
D
. et al.  
Structure, function, and pharmacology of glutamate receptor ion channels
.
Pharmacol Rev
 
2021
;
73
:
1469
1658
.

28.

Myers
 
SJ
,
Yuan
 
H
,
Perszyk
 
RE
. et al.  
Classification of missense variants in the N-methyl-d-aspartate receptor GRIN gene family as gain- or loss-of-function
.
Hum Mol Genet
 
2023
;
32
:
2857
2871
.

29.

Hanson
 
JE
,
Ma
 
K
,
Elstrott
 
J
. et al.  
GluN2A NMDA receptor enhancement improves brain oscillations, synchrony, and cognitive functions in Dravet syndrome and Alzheimer’s disease models
.
Cell Rep
 
2020
;
30
:
381
396.e4
.

30.

Kemp
 
JA
,
McKernan
 
RM
.
NMDA receptor pathways as drug targets
.
Nat Neurosci
 
2002
;
5
:
1039
1042
.

31.

Landrum
 
MJ
,
Kattman
 
BL
.
ClinVar at five years: delivering on the promise
.
Hum Mutat
 
2018
;
39
:
1623
1630
.

32.

Silk
 
M
,
Petrovski
 
S
,
Ascher
 
DB
.
MTR-viewer: identifying regions within genes under purifying selection
.
Nucleic Acids Res
 
2019
;
47
:
W121
W126
.

33.

Silk
 
M
,
Pires
 
DEV
,
Rodrigues
 
CHM
. et al.  
MTR3D: identifying regions within protein tertiary structures under purifying selection
.
Nucleic Acids Res
 
2021
;
49
:
W438
W445
.

34.

Quinodoz
 
M
,
Peter
 
VG
,
Cisarova
 
K
. et al.  
Analysis of missense variants in the human genome reveals widespread gene-specific clustering and improves prediction of pathogenicity
.
Am J Hum Genet
 
2022
;
109
:
457
470
.

35.

Sivley
 
RM
,
Dou
 
X
,
Meiler
 
J
. et al.  
Comprehensive analysis of constraint on the spatial distribution of missense variants in human protein structures
.
Am J Hum Genet
 
2018
;
102
:
415
426
.

36.

Pérez-Palma
 
E
,
May
 
P
,
Iqbal
 
S
. et al.  
Identification of pathogenic variant enriched regions across genes and gene families
.
Genome Res
 
2020
;
30
:
62
71
.

37.

Perszyk
 
RE
,
Kristensen
 
AS
,
Lyuboslavsky
 
P
. et al.  
Three-dimensional missense tolerance ratio analysis
.
Genome Res
 
2021
;
31
:
1447
1461
.

38.

Chiron
 
C
,
Marchand
 
MC
,
Tran
 
A
. et al.  
Stiripentol in severe myoclonic epilepsy in infancy: a randomised placebo-controlled syndrome-dedicated trial. STICLO study group
.
Lancet
 
2000
;
356
:
1638
1642
.

39.

Jen
 
JC
,
Graves
 
TD
,
Hess
 
EJ
. et al.  
Primary episodic ataxias: diagnosis, pathogenesis and treatment
.
Brain
 
2007
;
130
:
2484
2493
.

40.

Schoonjans
 
AS
,
Lagae
 
L
,
Ceulemans
 
B
.
Low-dose fenfluramine in the treatment of neurologic disorders: experience in Dravet syndrome
.
Ther Adv Neurol Disord
 
2015
;
8
:
328
338
.

41.

Mullier
 
B
,
Wolff
 
C
,
Sands
 
ZA
. et al.  
GRIN2B gain of function mutations are sensitive to radiprodil, a negative allosteric modulator of GluN2B-containing NMDA receptors
.
Neuropharmacology
 
2017
;
123
:
322
331
.

42.

Pierson
 
TM
,
Yuan
 
H
,
Marsh
 
ED
. et al.  
GRIN2A mutation and early-onset epileptic encephalopathy: personalized therapy with memantine
.
Ann Clin Transl Neurol
 
2014
;
1
:
190
198
.

43.

Chidambaram
 
S
,
Manokaran
 
RK
.
Favorable response to “Memantine” in a child with GRIN2B epileptic encephalopathy
.
Neuropediatrics
 
2022
;
53
:
287
290
.

44.

Xu
 
Y
,
Song
 
R
,
Chen
 
W
. et al.  
Recurrent seizure-related GRIN1 variant: molecular mechanism and targeted therapy
.
Ann Clin Transl Neurol
 
2021
;
8
:
1480
1494
.

45.

Han
 
W
,
Yuan
 
H
,
Allen
 
JP
. et al.  
Opportunities for precision treatment of GRIN2A and GRIN2B gain-of-function variants in Triheteromeric N-methyl-D-aspartate receptors
.
J Pharmacol Exp Ther
 
2022
;
381
:
54
66
.

46.

Tang
 
W
,
Liu
 
D
,
Traynelis
 
SF
. et al.  
Positive allosteric modulators that target NMDA receptors rectify loss-of-function GRIN variants associated with neurological and neuropsychiatric disorders
.
Neuropharmacology
 
2020
;
177
:
108247
. .

47.

Zhu
 
S
,
Paoletti
 
P
.
Allosteric modulators of NMDA receptors: multiple sites and mechanisms
.
Curr Opin Pharmacol
 
2015
;
20
:
14
23
.

48.

Addis
 
L
,
Virdee
 
JK
,
Vidler
 
LR
. et al.  
Epilepsy-associated GRIN2A mutations reduce NMDA receptor trafficking and agonist potency - molecular profiling and functional rescue
.
Sci Rep
 
2017
;
7
:
66
.

49.

Boßelmann
 
CM
,
Hedrich
 
UBS
,
Müller
 
P
. et al.  
Predicting the functional effects of voltage-gated potassium channel missense variants with multi-task learning
.
EBioMedicine
 
2022
;
81
:
104115
. .

50.

Heyne
 
HO
,
Baez-Nieto
 
D
,
Iqbal
 
S
. et al.  
Predicting functional effects of missense variants in voltage-gated sodium and calcium channels
.
Sci Transl Med
 
2020
;
12
:eeaay6848.

51.

Sevim Bayrak
 
C
,
Stein
 
D
,
Jain
 
A
. et al.  
Identification of discriminative gene-level and protein-level features associated with pathogenic gain-of-function and loss-of-function variants
.
Am J Hum Genet
 
2021
;
108
:
2301
2318
.

52.

Jung
 
S
,
Lee
 
S
,
Kim
 
S
. et al.  
Identification of genomic features in the classification of loss- and gain-of-function mutation
.
BMC Med Inform Decis Mak
 
2015
;
15
:
S6
.

53.

Liu
 
M
,
Watson
 
LT
,
Zhang
 
L
.
HMMvar-func: a new method for predicting the functional outcome of genetic variants
.
BMC Bioinformatics
 
2015
;
16
:
351
.

54.

Ge
 
F
,
Li
 
C
,
Iqbal
 
S
. et al.  
VPatho: a deep learning-based two-stage approach for accurate prediction of gain-of-function and loss-of-function variants
.
Brief Bioinform
 
2023
;
24
:
bbac535
.

55.

Swanger
 
SA
,
Chen
 
W
,
Wells
 
G
. et al.  
Mechanistic insight into NMDA receptor dysregulation by rare variants in the GluN2A and GluN2B agonist binding domains
.
Am J Hum Genet
 
2016
;
99
:
1261
1280
.

56.

Hansen
 
KB
,
Clausen
 
RP
,
Bjerrum
 
EJ
. et al.  
Tweaking agonist efficacy at N-methyl-D-aspartate receptors by site-directed mutagenesis
.
Mol Pharmacol
 
2005
;
68
:
1510
1523
.

57.

Chen
 
PE
,
Geballe
 
MT
,
Stansfeld
 
PJ
. et al.  
Structural features of the glutamate binding site in recombinant NR1/NR2A N-methyl-D-aspartate receptors determined by site-directed mutagenesis and molecular modeling
.
Mol Pharmacol
 
2005
;
67
:
1470
1484
.

58.

Maier
 
W
,
Schemm
 
R
,
Grewer
 
C
. et al.  
Disruption of interdomain interactions in the glutamate binding pocket affects differentially agonist affinity and efficacy of N-methyl-D-aspartate receptor activation
.
J Biol Chem
 
2007
;
282
:
1863
1872
.

59.

Laube
 
B
,
Hirai
 
H
,
Sturgess
 
M
. et al.  
Molecular determinants of agonist discrimination by NMDA receptor subunits: analysis of the glutamate binding site on the NR2B subunit
.
Neuron
 
1997
;
18
:
493
503
.

60.

Karczewski
 
KJ
,
Francioli
 
LC
,
Tiao
 
G
. et al.  
The mutational constraint spectrum quantified from variation in 141,456 humans
.
Nature
 
2020
;
581
:
434
443
.

61.

Chen
 
W
,
Tankovic
 
A
,
Burger
 
PB
. et al.  
Functional evaluation of a De novo GRIN2A mutation identified in a patient with profound global developmental delay and refractory epilepsy
.
Mol Pharmacol
 
2017
;
91
:
317
330
.

62.

XiangWei
 
W
,
Kannan
 
V
,
Xu
 
Y
. et al.  
Heterogeneous clinical and functional features of GRIN2D-related developmental and epileptic encephalopathy
.
Brain
 
2019
;
142
:
3009
3027
.

63.

Traynelis
 
SF
.
Software-based correction of single compartment series resistance errors
.
J Neurosci Methods
 
1998
;
86
:
25
34
.

64.

Dana
 
JM
,
Gutmanas
 
A
,
Tyagi
 
N
. et al.  
SIFTS: updated structure integration with function, taxonomy and sequences resource allows 40-fold increase in coverage of structure-based annotations for proteins
.
Nucleic Acids Res
 
2019
;
47
:
D482
D489
.

65.

Zhang
 
Y
,
Ye
 
F
,
Zhang
 
T
. et al.  
Structural basis of ketamine action on human NMDA receptors
.
Nature
 
2021
;
596
:
301
305
.

66.

Lomize
 
MA
,
Pogozheva
 
ID
,
Joo
 
H
. et al.  
OPM database and PPM web server: resources for positioning of proteins in membranes
.
Nucleic Acids Res
 
2012
;
40
:
D370
D376
.

67.

Pravda
 
L
,
Sehnal
 
D
,
Toušek
 
D
. et al.  
MOLEonline: a web-based tool for analyzing channels, tunnels and pores (2018 update)
.
Nucleic Acids Res
 
2018
;
46
:
W368
W373
.

68.

Grant
 
BJ
,
Rodrigues
 
APC
,
ElSawy
 
KM
. et al.  
Bio3d: an R package for the comparative analysis of protein structures
.
Bioinformatics
 
2006
;
22
:
2695
2696
.

69.

Dong
 
R
,
Pan
 
S
,
Peng
 
Z
. et al.  
mTM-align: a server for fast protein structure database search and multiple protein structure alignment
.
Nucleic Acids Res
 
2018
;
46
:
W380
W386
.

70.

Fedele
 
L
,
Newcombe
 
J
,
Topf
 
M
. et al.  
Disease-associated missense mutations in GluN2B subunit alter NMDA receptor ligand binding and ion channel properties
.
Nat Commun
 
2018
;
9
:
957
.

71.

Kabsch
 
W
,
Sander
 
C
.
How good are predictions of protein secondary structure?
 
FEBS Lett
 
1983
;
155
:
179
182
.

72.

Bastolla
 
U
,
Farwer
 
J
,
Knapp
 
EW
. et al.  
How to guarantee optimal stability for most representative structures in the protein data Bank
.
Proteins
 
2001
;
44
:
79
96
.

73.

Montanucci
 
L
,
Capriotti
 
E
,
Birolo
 
G
. et al.  
DDGun: an untrained predictor of protein stability changes upon amino acid variants
.
Nucleic Acids Res
 
2022
;
50
:
W222
W227
.

74.

Montanucci
 
L
,
Capriotti
 
E
,
Frank
 
Y
. et al.  
DDGun: an untrained method for the prediction of protein stability changes upon single and multiple point variations
.
BMC Bioinformatics
 
2019
;
20
:
335
.

75.

Henikoff
 
S
,
Henikoff
 
JG
.
Amino acid substitution matrices from protein blocks
.
Proc Natl Acad Sci USA
 
1992
;
89
:
10915
10919
.

76.

Hopf
 
TA
,
Ingraham
 
JB
,
Poelwijk
 
FJ
. et al.  
Mutation effects predicted from sequence co-variation
.
Nat Biotechnol
 
2017
;
35
:
128
135
.

77.

Cheng
 
J
,
Novati
 
G
,
Pan
 
J
. et al.  
Accurate proteome-wide missense variant effect prediction with AlphaMissense
.
Science
 
2023
;
381
:
eadg7492
.

78.

Frazer
 
J
,
Notin
 
P
,
Dias
 
M
. et al.  
Disease variant prediction with deep generative models of evolutionary data
.
Nature
 
2021
;
599
:
91
95
.

79.

Ioannidis
 
NM
,
Rothstein
 
JH
,
Pejaver
 
V
. et al.  
REVEL: an ensemble method for predicting the pathogenicity of rare missense variants
.
Am J Hum Genet
 
2016
;
99
:
877
885
.

80.

Lester
 
RA
,
Jahr
 
CE
.
NMDA channel behavior depends on agonist affinity
.
J Neurosci
 
1992
;
12
:
635
643
.

81.

Vance
 
KM
,
Simorowski
 
N
,
Traynelis
 
SF
. et al.  
Ligand-specific deactivation time course of GluN1/GluN2D NMDA receptors
.
Nat Commun
 
2011
;
2
:
294
.

82.

Xu
 
Y
,
Song
 
R
,
Perszyk
 
RE
. et al.  
De novo GRIN variants in M3 helix associated with neurological disorders control channel gating of NMDA receptor
.
Cell Mol Life Sci
 
2024
;
81
:
153
.

83.

Ashley
 
EA
.
Towards precision medicine
.
Nat Rev Genet
 
2016
;
17
:
507
522
.

84.

Adzhubei
 
IA
,
Schmidt
 
S
,
Peshkin
 
L
. et al.  
A method and server for predicting damaging missense mutations
.
Nat Methods
 
2010
;
7
:
248
249
.

85.

Rentzsch
 
P
,
Witten
 
D
,
Cooper
 
GM
. et al.  
CADD: predicting the deleteriousness of variants throughout the human genome
.
Nucleic Acids Res
 
2019
;
47
:
D886
D894
.

86.

Sim
 
NL
,
Kumar
 
P
,
Hu
 
J
. et al.  
SIFT web server: predicting effects of amino acid substitutions on proteins
.
Nucleic Acids Res
 
2012
;
40
:
W452
W457
.

87.

Stead
 
LF
,
Wood
 
IC
,
Westhead
 
DR
.
KvSNP: accurately predicting the effect of genetic variants in voltage-gated potassium channels
.
Bioinformatics
 
2011
;
27
:
2181
2186
.

88.

Li
 
B
,
Mendenhall
 
JL
,
Kroncke
 
BM
. et al.  
Predicting the functional impact of KCNQ1 variants of unknown significance
.
Circ Cardiovasc Genet
 
2017
;
10
:
e001754
.

89.

Amador
 
A
,
Bostick
 
CD
,
Olson
 
H
. et al.  
Modelling and treating GRIN2A developmental and epileptic encephalopathy in mice
.
Brain
 
2020
;
143
:
2039
2057
.

90.

Li
 
J
,
Zhang
 
J
,
Tang
 
W
. et al.  
De novo GRIN variants in NMDA receptor M2 channel pore-forming loop are associated with neurological diseases
.
Hum Mutat
 
2019
;
40
:
2393
2413
.

Author notes

Ludovica Montanucci and Tobias Brünger joint first authors.

This is an Open Access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0/), which permits unrestricted reuse, distribution, and reproduction in any medium, provided the original work is properly cited.