Isoform-specific Inhibition of N-methyl-D-aspartate Receptors by Bile Salts

The N-methyl-D-aspartate subfamily of ionotropic glutamate receptors (NMDARs) is well known for its important roles in the central nervous system (CNS), e.g. learning and memory formation. Besides the CNS, NMDARs are also expressed in numerous peripheral tissues including the pancreas, kidney, stomach, and blood cells, where an understanding of their physiological and pathophysiological roles is only evolving. Whereas subunit composition increases functional diversity of NMDARs, a great number of endogenous cues tune receptor signaling. Here, we characterized the effects of the steroid bile salts cholate and chenodeoxycholate (CDC) on recombinantly expressed NMDARs of defined molecular composition. CDC inhibited NMDARs in an isoform-dependent manner, preferring GluN2D and GluN3B over GluN2A and GluN2B receptors. Determined IC50 values were in the range of bile salt serum concentrations in severe cholestatic disease states, pointing at a putative pathophysiological significance of the identified receptor modulation. Both pharmacological and molecular simulation analyses indicate that CDC acts allosterically on GluN2D, whereas it competes with agonist binding on GluN3B receptors. Such differential modes of inhibition may allow isoform-specific targeted interference with the NMDAR/bile salt interaction. In summary, our study provides further molecular insight into the modulation of NMDARs by endogenous steroids and points at a putative pathophysiological role of the receptors in cholestatic disease.


Results
CDC inhibits preferentially the GluN subtypes 2D and 3B. The effects of the bile salts cholate and CDC on GluNs were probed on recombinant receptors of defined subunit composition expressed in Xenopus laevis oocytes. For initial experiments, agonists were applied at saturating concentration to ensure maximum activation of GluN2A, GluN2B, GluN2D, and GluN3B, respectively 4,19 . Cholate inhibited steady-state currents only moderately irrespective of receptor subtype, while CDC was more potent and reduced preferentially GluN2D and GluN3B currents (Fig. 1, Table S1). Inhibition by 100 µM CDC ranged between 7 ± 2% of GluN2A and 13 ± 1% of GluN2B steady-state currents, whereas it amounted to 37 ± 2% of GluN2D and even to 58 ± 5% of GluN3B currents. The bile salt effects were specific for GluNs, as no inhibition of glutamate-induced steady-state kainate receptor GluK2 current (n = 3, after block of desensitization with ConA) and only a minor reduction of α-amino-3-hydroxy-5-methyl-4-isoxazolepropionic acid (AMPA) receptor GluA1/GluA2 currents (7 ± 1%, n = 6) were observed after application of 100 µM CDC. Somewhat unexpected, taurine or glycine conjugation of bile salts did neither change their overall potencies in inhibiting GluNs nor their subunit preferences.

Allosteric inhibition of GluN2D by tauro-CDC.
To dissect the mechanism of GluN2D inhibition by bile salts, we recorded dose-response relations for CDC at both saturating and EC 50 NMDAR agonist concentrations 20,21 . We used routinely taurine-conjugated CDC in these experiments, because of higher water solubility and hence more feasible handling. As depicted in Fig. 2a, maximum inhibition of glutamate/glycine-induced currents by tauro-CDC was significantly smaller at saturating than at EC 50 agonist concentrations (n = 8-13; p = 0.01 and 0.03 for glutamate and glycine, respectively).
Conversely, dose-response curves for glutamate and glycine in the absence or presence of 180 μM tauro-CDC, which was calculated to be the mean IC 50 (179 ± 11 μM, n = 13, saturating agonist concentration), showed that the maximum agonist-induced currents were significantly reduced by the bile salt (Fig. 2b,c; n = 7-13, p = 0.0002 and 0.003 for glycine and glutamate, respectively). There were neither significant differences between the EC 50 values for glutamate or glycine at varying concentrations of tauro-CDC (n = 7-13, p = 0.3 and 0.2 for glutamate and glycine respectively) nor were there different IC 50 values for tauro-CDC at varying concentrations of glycine (n = 8-13; p = 0.1; Table S2). Only the IC 50 value for tauro-CDC at the glutamate EC 50 concentration was significantly smaller compared to inhibition at a saturating glutamate concentration (n = 11-13, p = 0.009).
Application of 100 µM tauro-CDC in the absence of NMDAR agonists showed no effects on GluN2D or GluN3B receptors (n = 4) excluding partial agonist action of bile salts on NMDA receptors (data not shown). Taken together, these results strongly suggest that tauro-CDC inhibits GluN2D receptors by an allosteric mechanism.
Competitive inhibition of GluN3B by tauro-CDC. Next, we sought to investigate the mechanism of action of tauro-CDC-mediated inhibition of GluN3B receptors. As depicted in Fig. 3a, the mean IC 50 for tauro-CDC was significantly higher in the presence of 10 μM glycine than in the presence of 5 μM glycine (83 ± 8 μM vs. 44 ± 7 μM, n = 6, p = 0.007). The maximum block by tauro-CDC, however, did not depend on the glycine concentration (95 ± 3% vs. 92 ± 3%, n = 6, p = 0.5). These data indicate that tauro-CDC inhibits GluN3B receptors -other than GluN2D receptors -in a competitive manner. We also recorded glycine dose-response curves in the presence or absence of 100 μM tauro-CDC (n = 6-11, Fig. 3b). Known receptor desensitization at higher glycine concentrations prevented calculation of precise EC 50 values 4,22 . However, as depicted in Fig. 3b, the bell-shaped concentration-response curve shifted to higher glycine concentrations in the presence of tauro-CDC.
Tauro-CDC has no effect on ion selectivity. Bile salts are known to integrate into the plasma membrane of cells 23 , which might impact ion channel pores including their selectivity filters. Therefore, we checked for potential bile salt effects on GluN receptor ion selectivity. Current-voltage (I/V) relationships of agonist-induced currents were recorded for GluN2D and GluN3B in the absence and presence of 100 µM tauro-CDC (n = 10-11). For better comparison of I/V curves, individual oocyte data were normalized to their maximum inward current at −140 mV. As shown in Fig. 4, tauro-CDC did neither change the shape of I/V relationships nor did it shift the reversal potentials of the receptor currents. These results render a change in ion selectivity of the two GluNs by bile salts unlikely.
Bile salt inhibition of GluNs when activated by alternative agonists. As the inhibition of GluNs by bile salts required concentrations that are hardly reached in the CNS, we tested whether tauro-CDC also reduced currents elicited by the alternative non-neuronal NMDAR agonists L-homocysteic acid (L-HCA) and quinolinic acid (QN) 5 .
As shown in Fig. 5, GluN2D currents induced by application of L-HCA at its reported EC 50 of 13 µM 24 were blocked by 30 ± 2% (n = 12), whereas currents elicited by QN at its reported EC 50 of 7.2 mM 24 were reduced by 26 ± 2% (n = 7) when co-applying tauro-CDC at 100 µM. Lowering the agonist concentration of QN to 1 mM did not change inhibition by tauro-CDC, which was still 29 ± 2% (n = 7). For all these experiments involving the non-neuronal NMDAR agonists L-HCA and QN on GluN2D receptors, 10 µM glycine was used as co-agonist. GluN3B could not be activated by L-HCA, even at high concentrations of up to 1 mM (n = 4). In contrast, currents activated by 7.2 mM QN were reduced by about 68 ± 4% (n = 5), with QN showing no effect on uninjected oocytes (n = 3).
Binding modes of tauro-CDC in GluN2D LBD and GluN3B LBD . The pharmacological experiments indicated that tauro-CDC acts as a non-competitive allosteric antagonist of GluN2D, but as a competitive antagonist of GluN3B. In order to assess the molecular origin of isoform-specific inhibition of NMDARs by tauro-CDC, we performed induced fit docking experiments of tauro-CDC in respective ligand binding domains (LBDs) of the NMDAR isoforms (Fig. 6a,b), which allows to consider receptor flexibility in terms of side chain rearrangements.
For both GluN2D LBD and GluN3B LBD , binding modes of tauro-CDC (Fig. 6c,d) were predicted with similar docking scores (IFDScore 25 ) of ~7.4 kcal·mol −1 . In agreement, both binding modes were of similar geometry with an in-place root-mean-square deviation (RMSD) of the ligand atoms of 2.98 Å after superimposition of the C α atoms of the protein structures (RMSD: 1.45 Å). The similarity in binding mode geometry was paralleled by a similar interaction profile in both complexes (Fig. 6c,d).
Hence, the static structures obtained from docking experiments cannot explain the observed isoform specificity of NMDAR inhibition by tauro-CDC.
www.nature.com/scientificreports www.nature.com/scientificreports/ Structural variability of tauro-CDC/GluN2D LBD and tauro-CDC/GluN3B LBD binding modes. We next asked whether the predicted binding modes of tauro-CDC in GluN2D LBD or GluN3B LBD structurally vary over time scales accessible by molecular dynamics (MD) simulations. We therefore carried out ten MD simulations (S 1 -S 10 ) of 500 ns length each for both the tauro-CDC/GluN2D LBD and tauro-CDC/GluN3B LBD systems and quantified the structural deviation from the initial tauro-CDC pose within each simulation and between independent simulations. Similarly, we quantified the conformational variability of tauro-CDC to account for potential alternative binding modes. Taking all simulations but S 9 into account, the average RMSD of the tauro-CDC pose in GluN2D LBD was 2.86 ± 0.35 Å higher than in GluN3B LBD ( Supplementary Fig. 3a). Simulation S 9 of GluN2D LBD was omitted from this and all subsequent analyses when tauro-CDC unbound after ~225 ns ( Supplementary Fig. 3b). No unbinding event was observed in the simulations of tauro-CDC/GluN3B LBD .
To conclude, these results reveal that the predicted binding mode of tauro-CDC in GluN2D LBD is structurally less stable than in GluN3B LBD . In one simulation of tauro-CDC in GluN2D LBD , tauro-CDC even unbound. Note that simulating ligand unbinding events is challenging also for ligands with moderate affinities due to low off-rates 27 , which may explain why the unbinding occurred only in one simulation. As to GluN3B LBD , the results are in line with our experimental data, according to which tauro-CDC must bind to the glycine binding site in order to competitively inhibit GluN3B-containing NMDARs.

Contributions of hydrogen bonds to tauro-CDC/GluN2D LBD and tauro-CDC/GluN3B LBD complex stability.
To assess whether a difference in structural variability between tauro-CDC/GluN2D LBD and tauro-CDC/GluN3B LBD may be attributed to differences in intermolecular interactions within the complexes, we quantified the amount of hydrogen bonds formed between tauro-CDC and an arginine residue that is crucial for binding of the agonist glutamate to GluN2D LBD (Arg543) and of glycine to GluN3B LBD (Arg538) over the course www.nature.com/scientificreports www.nature.com/scientificreports/ of all MD simulations. Additionally, we quantified the total amount of hydrogen bonds between tauro-CDC and GluN2D LBD or GluN3B LBD . We focused here on identifying hydrogen bonds as interactions that contribute to the stability of the tauro-CDC/GluN3B LBD complex in order to predict mutations for ex vivo validation of a suggested tauro-CDC binding mode.
To conclude, hydrogen bond formation between tauro-CDC and Arg543 in GluN2D LBD occurred less frequently than between tauro-CDC and Arg538 in GluN3B LBD ; the same held true when considering hydrogen bond formation to the complete LBD. The tauro-CDC/GluN3B LBD complex was particularly stabilized by Arg538 and, to a lesser extent, Asp447, whereas no particular amino acid contributed to stabilizing tauro-CDC/ GluN2D LBD .
Effective binding energies of tauro-CDC/GluN2D LBD and tauro-CDC/GluN3B LBD complexes. To complement our findings from the structural analysis of tauro-CDC complexes, we quantified the effective binding energy of tauro-CDC to GluN2D LBD or GluN3B LBD with the end-point effective energy approaches MM/PBSA and MM/GBSA.
The effective binding energy computed by MM/GBSA for tauro-CDC/GluN2D LBD (mean: −33.53 ± 0.67 kcal·mol −1 ) was significantly less negative than for tauro-CDC/GluN3B LBD (mean:  (Fig. 9d); higher dielectric constants have been used before to account for highly charged binding sites 28 , as given in the case of the LBDs.
To conclude, both approaches revealed more negative effective energies of tauro-CDC binding to GluN3B LBD than to GluN2D LBD . Thus, they indicate a stronger binding of tauro-CDC to GluN3B LBD .
Identification of an alternative binding site in the GluN1 LBD /GluN2D LBD interface. We further aimed to identify an alternative binding mode of tauro-CDC to a potentially allosteric site in the GluN1 LBD / GluN2D LBD tetramer. To do so, we first detected potential binding pockets in a conformational ensemble obtained from 500 ns of accelerated MD simulations of the tetrameric GluN1 LBD /GluN2D LBD interface (Fig. 10a) using the MDpocket software 29 . A single binding pocket at the GluN1/GluN2D LBD interface was identified, and tauro-CDC was docked into the receptor conformation in which the pocket showed the highest mean local hydrophobic density 29,30 (Fig. 10b). The identified pocket largely overlaps with the allosteric binding pocket of the negative allosteric modulator MPX-007 in the GluN1/GluN2A interface 31 .
During 500 ns of conventional MD simulations, tauro-CDC displayed an average RMSD value of 1.90 ± 0.01 Å (2.17 ± 0.01 Å without fitting, calculated on the last 250 ns of the simulation), indicating a low structural variability of the predicted binding mode. The binding mode is stabilized by interactions of the acid moiety with the positively charged residues Lys531 and Lys769 in the GluN1 subunit and Lys797 and Arg798 in the GluN2D subunit (Fig. 10d); Leu808 in the GluN2D subunit makes a hydrophobic contact with the cholane scaffold (Fig. 10d). Functional data on GluN1/GluN2A NMDARs indicates that Phe754 in GluN1 and Val783 in GluN2A form a molecular switch that mediates allosteric inhibition 31 . The corresponding Leu808 in GluN2D may assume a similar role as Val783 in GluN2A.
These tetramer simulations suggest that the interface between the GluN1 and GluN2D subunit forms a transient pocket, to which tauro-CDC can bind.

Point mutants.
In order to verify the interaction of tauro-CDC with residues of the GluN3B LBD, Arg538 would be the primary residue for mutation studies due to its pronounced involvement in hydrogen bonds with the bile salt (Fig. 8b). However, we did not mutate Arg538 because this residue is reported to be important for www.nature.com/scientificreports www.nature.com/scientificreports/ agonist binding 32 . Rather, we mutated Asp447, Tyr505, and Ser701 (Fig. 8b) that also engage in hydrogen bonds with the bile salt. Yet, as these interactions are transient and, hence, likely rather weak, a decrease in the extent of tauro-CDC block would indicate an interaction of the bile salt with the mutated residue, while no effect would not exclude such an interaction due to potential enthalpy-entropy compensation effects 33 .
Mutation of Asp447 to Ala or Arg had neither an effect on agonist-induced current amplitudes (n = 2-3 dose-response curves) nor on the extent of tauro-CDC block in comparison to wild type receptor block (100 µM tauro-CDC, n = 9-12 at 10 µM Gly, n = 9-10 at 5 µM Gly). Mutation of Tyr505 to an Ala resulted in GluN3B receptors that could not be activated by its agonist glycine anymore (n = 9). Mutation of Ser701 to Ala strongly reduced glycine-induced currents from 835 ± 85 nA (n = 30) to 106 ± 12 nA (n = 29). While we observed an increase of 15% in CDC and tauro-CDC block (both 100 µM) for this mutant (n = 11-13), which might indicate an influence on CDC and tauro-CDC binding, it may also well be due to reduced glycine binding.
To conclude, none of the mutants showed a significant influence on tauro-CDC block such that no further information about binding (or non-binding) of tauro-CDC to GluN3B LBD may be inferred. www.nature.com/scientificreports www.nature.com/scientificreports/

Discussion
Here, we have shown that the bile salt CDC effectively inhibits NMDARs with a preference on GluN2D and GluN3B containing receptors. Neither did conjugation of CDC per se nor the type of conjugated amino salttaurine or glycine -change the overall inhibitory potency of the bile salt or shift its subunit preference. Based on detailed pharmacological analyses and supported by MD experiments, we infer from our data that CDC inhibits NMDARs by different mechanisms that are allosteric inhibition of GluN2D and competitive block of GluN3B containing receptors.
The observed subtype preference of NMDAR inhibition by CDC is in good agreement with the one reported for structurally similar neurosteroids. Thus, pregnenolone sulfate has been shown to significantly potentiate GluN2A and GluN2B containing NMDARs but has much lower efficacy on GluN2C and GluN2D receptors [6][7][8][9] . Also, the apparent insignificance of amino acid conjugates in the interaction of bile salts with receptors is In panels c and e, data for MD simulation S 9 is not shown because here tauro-CDC unbinds after 225 ns (c.f. Supplementary Fig. 3b), and the corresponding data was excluded from statistical analysis. (2019) 9:10068 | https://doi.org/10.1038/s41598-019-46496-y www.nature.com/scientificreports www.nature.com/scientificreports/ paralleled in CDC-mediated activation of the bile acid receptor TGR5, which does not discriminate between conjugation substrates 34 . Intriguingly, however, our pharmacological experiments indicated different modes of inhibition for different GluN subtypes by CDC: allosteric inhibition of GluN2D and competitive block of GluN3B. In line with this suggestion, MD simulations predicted less structural stability for tauro-CDC binding in the orthosteric pocket of GluN2D LBD than in the one of GluN3B LBD . Effective binding energy computations comparing tauro-CDC/GluN2D LBD with tauro-CDC/GluN3B LBD complexes performed by two different approaches consistently revealed more negative effective binding energies of tauro-CDC binding to GluN3B LBD than to GluN2D LBD , thus supporting stronger binding of tauro-CDC to GluN3B LBD .
Previous studies have shown that inclusion of configurational entropy is crucial for calculating absolute binding free energies 35,36 . In the present study, however, we were rather interested in relative binding free energies and, thus, decided to neglect contributions due to changes in the configurational entropy of the ligand or the receptor upon complex formation in order to avoid introducing additional uncertainty in the computations 35,37,38 . Finally, MD simulations of the GluN1 LBD /GluN2D LBD tetramer suggested that the interface between the GluN1 and GluN2D subunit forms a transient pocket, to which tauro-CDC can bind, which may explain the allosteric mechanism of action in this GluN subtype. We selected this interfacial pocket based on a combined approach of geometry-based cavity detection and druggability prediction that has been shown to be able to successfully identify binding-competent pocket conformations from an MD ensemble 29 . Notably, the detected pocket largely overlaps with the binding site of the negative allosteric modulator MPX-007 in the GluN1/GluN2A interface 29 . In order to verify binding of tauro-CDC to the orthosteric pocket of GluN3B LBD but not GluN2D LBD , we quantified hydrogen bond formation between tauro-CDC and neighboring residues based on geometric criteria over the course of the MD simulations. The results show that hydrogen bonding occurred less frequently between tauro-CDC and the GluN2D LBD than between tauro-CDC and the GluN3B LBD , in agreement with the above structural and energetic analyses of MD simulations. The tauro-CDC/GluN3B LBD complex is predominantly stabilized by a hydrogen bond to Arg538 and, to a lesser extent, Asp447. Of the residues mutated to alanine in the GluN3B LBD , positions 505 and 701 resulted in GluN3B receptors that could either not or hardly be activated by glycine anymore. For Ser701, it was reported that it forms hydrogen bonds with the glycine carboxylate group 32 , so it is not surprising that mutation of this site affected channel function. Mutation of Asp447 to an Ala or Arg had neither an effect on agonist-induced current amplitudes nor on the extent of tauro-CDC block, which does not exclude that Asp447 interacts with tauro-CDC, as enthalpy-entropy compensation effects upon mutation may account for the loss of interaction in the mutant 33 . In particular, we speculate that the Asp447Ala mutation leads to a loss in binding enthalpy upon disruption of the tauro-CDC-4Asp47 hydrogen bond, but this loss in binding enthalpy could be compensated by both a gain in entropy of the ligand due to an increase in translational and rotational degrees of freedom and a gain in entropy of the receptor due to increased conformational freedom of the residues surrounding the mutation site. The latter effect might especially be relevant due to the stabilizing interaction between Asp447 and Arg474, which would break upon mutation of Asp447 to alanine. Similarly, in the Asp447Arg mutant, the interaction between tauro-CDC and receptor could be preserved if arginine acted as a hydrogen bond donor for the 3-OH group of tauro-CDC, which would in turn act as a hydrogen bond acceptor. Moreover, repulsive forces between Arg447 in the Arg447Arg mutant and Arg474 could lead to a gain in receptor entropy.
The inhibition of GluNs by bile salts required concentrations that are hardly reached in the CNS, where CDC and cholate exists in the nmol/g range 39,40 . However, NMDARs are also widely expressed in non-neural peripheral tissues. These receptors have many distinct physiological and pathophysiological roles, and there is evidence that peripheral NMDARs may use alternative agonists such as L-HCA and QN 5 . We could show that tauro-CDC also reduced currents elicited by these alternative NMDAR agonists, indicating potential modulation of peripheral NMDAR function by bile salts. Whereas under physiological conditions only low µM amounts of bile salts are present in serum and urine, in cholestatic disease states like inherited progressive cholestasis, severe forms www.nature.com/scientificreports www.nature.com/scientificreports/ of intrahepatic cholestasis of pregnancy or acute hepatitis, plasma concentrations may rise to hundred µM and more [10][11][12] . As blood cells including platelets express NMDARs, it may well be hypothesized that bile salt mediated receptor inhibition contributes to the reported platelet activation defects and impaired thrombus formation in cholestatic liver disease 41 . Despite the observed preference of bile salts for specific receptor isoforms, it remains difficult to predict the exact physiological or pathophysiological consequences of their action on NMDAR  67 . In panels a-c, data for MD simulation S 9 is not shown because here, tauro-CDC unbinds after 225 ns (c.f. Supplementary Fig. 3b), and the corresponding data was excluded from statistical analysis. (2019) 9:10068 | https://doi.org/10.1038/s41598-019-46496-y www.nature.com/scientificreports www.nature.com/scientificreports/ signaling. The net effect of bile salts on NMDARs will eventually depend on the overall expression, stoichiometric assembly, and posttranscriptional and posttranslational modifications of individual GluN isoforms in a given cell type 42 . On a pharmacological perspective, however, the distinct modes of NMDAR inhibition by bile salts, i.e. allosteric versus competitive, may be exploited for developing targeted interference strategies holding isoform specificity.

Electrophysiological recordings. Two-electrode voltage clamp recordings of oocyte current responses
were performed 3-5 days after cRNA injection at −70 mV holding potential with a Turbo Tec-03X amplifier (npi electronic, Germany) controlled by Pulse software (HEKA, Germany). Electrodes were filled with 3 M KCl and had resistances of 0.5-1.5 MΩ. Oocytes were continuously superfused with calcium-free Ba 2+ -Ringer's solution (in mM: 115 NaCl, 2.5 KCl, 1.8 BaCl 2 and 10 HEPES, pH 7.2) to prevent the activation of endogenous Ca 2+ -gated chloride channels. Agonists and bile salts were applied for 20-40 s by superfusion. Agonist and bile salt blocker potencies were determined by recording current responses induced by the application of increasing agonist concentrations or increasing bile salt concentration to the same oocyte. Agonist-induced currents were normalized to the maximum response under saturating agonist concentration. In case of inhibition by bile salts the relative reduction in agonist-induced current was determined. Data from each oocyte was fitted separately to the Hill equation using Origin 9 software. The resulting EC 50 or IC 50 values were averaged. Current-voltage relationships were determined by ramping the holding potential from −140 mV to +50 mV corrected for background conductivities. Data are reported as mean ± standard error of the mean (SEM). Statistical significance was determined using the unpaired, two-tailed Student's t-test.

Molecular modeling and molecular dynamics simulations. Preparation of tauro-CDC/GluN2D LBD
and tauro-CDC/GluN3B LBD complex structures. In order to generate full-length structures of rat GluN2D and mouse GluN3B LBDs (from now on referred to as GluN2D LBD and GluN3B LBD ), we first generated a multiple sequence alignment of all rat and mouse GluN1, GluN2, and GluN3 subunit sequences using the MAFFT 43 server with default settings applied. The alignment was further refined using the GLProbs software 44 with two passes of consistency transformation and 100 passes of iterative refinement. The resulting alignment was then used as a template to manually align the sequence portions resolved in the selected template structures for GluN2D LBD (PDB ID: 3OEK 45 48 , keeping all other settings at their default values. The conformer with the lowest potential energy that did not display intramolecular hydrogen bonding was subjected to quantum mechanical geometry optimization at the HF/6-31 G(d) level of theory in the GAUSSIAN 09 software (Revision A.02) 49 . A tight convergence criterion for the self-consistent field iteration process was set, and computation of electrostatic potential points with a density of ~1.68 pt/au 2 (6 pt/Å 2 ) was invoked. The geometry-optimized structure was visually inspected for the absence of intramolecular hydrogen bonds and subjected to the standard RESP procedure 50,51 as implemented in the Antechamber set of programs to obtain the atomic point charges.
Binding modes for tauro-CDC in GluN2D LBD and GluN3B LBD were generated using the standard Induced Fit protocol implemented in the Schrödinger Suite 25 . The docking grid was centered on the coordinates of the bound ligand in the GluN2D and GluN3B template structure, respectively. In the initial Glide docking stage, the van der Waals radii of the receptor and ligand atoms were reduced to half their initial value, and a maximum of 20 poses were allowed to be carried over to the Prime (Schrödinger Release 2018-3: Prime, Schrödinger, LLC, New York, NY, 2018) refinement stage. For each generated pose, the receptor residues within 5 Å were energy-minimized and their side chain rotamers optimized. During the final Glide (Schrödinger Release 2018-3: Glide, Schrödinger, LLC, New York, NY, 2018) 52,53 docking stage, tauro-CDC was redocked into all receptor conformers using XP precision 54 . The tauro-CDC/GluN2D LBD and tauro-CDC/GluN3B LBD complexes with the lowest IFD score 25 were considered for all subsequent calculations.
Molecular dynamics simulations of tauro-CDC/GluN2D LBD and tauro-CDC/GluN3B LBD complexes. All molecular dynamics (MD) simulations were performed using the mixed single precision/fixed precision GPU (CUDA) version of PMEMD 55 in the Amber14 suite of molecular simulation programs 56 . Hydrogen mass repartitioning 57 was employed to enable a time step of 4 fs for integration. The Langevin thermostat 58 with a collision frequency of 0.01 ps −1 and a target temperature of T = 300 K was used for temperature control. Covalent bonds involving hydrogen atoms were constrained using the SHAKE algorithm 59 . Long-range electrostatic interactions were (2019) 9:10068 | https://doi.org/10.1038/s41598-019-46496-y www.nature.com/scientificreports www.nature.com/scientificreports/ estimated using the Particle Mesh Ewald method 60 , and a cutoff of 8 Å was used for short-range electrostatics and van der Waals forces.
All following steps were performed 10 times for each of the two complexes (tauro-CDC/GluN2D LBD and tauro-CDC/GluN3B LBD ) so as to obtain 10 independent simulations per complex. The initial structures were minimized using a three-step procedure. First, the coordinates of all solute molecules were restrained by a harmonic potential with a force constant of 2.0 kcal mol −1 Å −2 while 2,000 steps of steepest descent minimization followed by 3,000 steps of conjugate gradient minimization were carried out. This step was repeated with the restraints switched from the solute to the solvent molecules. During the following final minimization step, the restraints were removed, and 3,000 steps of steepest descent minimization followed by 7,000 steps of conjugate gradient minimization were performed. 20 ps of NVT-MD (the solute was restrained with a force constant of 2.0 kcal mol −1 Å −2 ) were carried out while heating the system from 0 K to 300 K, followed by additional 5 ps of NVT-MD at 300 K. Density adaptation of the system was achieved by 75 ps of NPT-MD (solute restrained, force constant: 2.0 kcal mol −1 Å −2 ). An additional 1.7 ns of restrained NPT-MD were carried out before switching to the NVT ensemble. Here, 3.2 ns of restrained MD were performed prior to the start of the production phase, with harmonic restraints (force constant: 2 kcal mol −1 Å −2 ) applied to only those C α atoms that are closest to the center of mass of a secondary structure element (α-helix or β-sheet). The subsequent production phase consisted of 500 ns of NVT-MD (restraints as in the final NVT step of the equilibration phase), resulting in an aggregate simulation time of 5 μs per complex (10 μs in total). Coordinates for analysis and post-processing were saved every 20 ps. Post-processing and analysis of the MD trajectories was performed in CPPTRAJ 61 as implemented in AmberTools15. Trajectories were visually inspected with VMD 62 . The two-tailed Student's t-test was used to determine statistically significant differences in time-averaged quantities between GluN2D and GluN3B.
Calculation of effective binding energies. Effective binding energies (ΔG eff ) of tauro-CDC to GluN2D LBD and GluN3B LBD were calculated using the MM/PBSA and MM/GBSA approaches as implemented in the MMPBSA. py module 63 in AmberTools15 56 . All calculations were performed on 500 snapshots per trajectory extracted at regular intervals of 1 ns, corresponding to a total of 5,000 snapshots per approach for both the tauro-CDC/ GluN2D LBD complex and the tauro-CDC/GluN3B LBD complex. The MM/PBSA calculations were performed with values of 1 and 4 for the internal dielectric constant (ε int ), to test for robustness of the relative ΔG eff with respect to this parameter 28 . The polar part of the contribution to the solvation free energy was calculated by the linear Poisson-Boltzmann equation using an ionic strength of 100 mM and an external dielectric constant (ε ext ) of 80. Solutions for the linear PB equation were computed with 1,000 iterations on a cubic grid with 0.5 Å spacing between grid points. The nonpolar part of the contribution to the solvation free energy was considered proportional to the solvent accessible surface area (SASA) and computed with the equation ΔG nonpol,solv = γ(SASA) + β, where γ and β are parameterized constants 64 . The SASA was calculated with a solvent probe radius of 1.4 Å, Tan & Luo radii 65 for the protein, and mbondi radii 66 for the ligands. The MM/GBSA calculations were performed using the modified OBC model (GB OBC II) 67 , an ionic strength of 100 mM and default values for γ (0.005 kcal mol −1 Å −2 ) and β (0.0 kcal mol −1 ).
Results are expressed as ΔG eff ± standard error of the mean (SEM), where ΔG eff is the average effective binding energy over all snapshots and the overall SEM is estimated as the SEM of the single trajectories propagated over all trajectories. The two-tailed Student's t-test was used to determine statistically significant differences in ΔG eff between GluN2D and GluN3B.
Preparation of GluN1 LBD /GluN2D LBD tetramer structures. The structure of the rat GluN1 LBD (GluN1 LBD ) was prepared analogously to the structure of GluN2D LBD (see above for further details), using the crystal structure of the rat GluN1 LBD (PDB ID: 1PB7 68 ) as a template to which only missing residues and loops were added during the modeling step. To construct the LBD tetramer from the individual LBD structures, the coordinates of the respective LBD monomers were superimposed onto their corresponding monomer units in the structure of the GluN1/GluN2D NMDA receptor (PDB ID: 5FXH 69 ) (Fig. 10a).
Accelerated molecular dynamics simulation of the GluN1 LBD /GluN2D LBD tetramer. A single, 500 ns long accelerated MD (aMD) simulation 70 of the GluN1 LBD /GluN2D LBD tetramer was prepared, heated, and equilibrated analogously to the classical MD simulations of the tauro-CDC/GluN2D LBD and tauro-CDC/GluN3B LBD complex (see above). During the production phase, bias parameters of α D = 781.0 kcal mol −1 and α P = 50,035.0 kcal mol −1 , a total dihedral energy threshold of E threshD = 18,950.0 kcal mol −1 , and a total potential energy threshold of E thresh P = −341,718.0 kcal mol −1 were set on the basis of empirical formulas that were previously used to estimate these parameters 70,71 .
Preparation of tauro-CDC//GluN1 LBD /GluN2D LBD complex structures. To identify a possible allosteric binding site in the GluN1 LBD /GluN2D LBD tetramer, the MDpocket software 29 was applied on all 25,000 snapshots of the 500 ns long aMD trajectory (see above). A binding pocket at the GluN1/GluN2D LBD interface was identified into which tauro-CDC was docked, using the snapshot in which the selected pocket displayed the highest mean local hydrophobic density 29,30 (Supplementary Fig. 2).
Molecular dynamics simulations of the tauro-CDC//GluN1 LBD /GluN2D LBD complex. A single 500 ns long conventional MD simulation of the tauro-CDC//GluN1 LBD /GluN2D LBD tetramer complex was generated and analyzed analogously to the classical MD simulations of the tauro-CDC/GluN2D LBD and tauro-CDC/GluN3B LBD complex (see above).