Ethanol-induced conformational fluctuations of NMDA receptors

ABSTRACT Alcohol addiction ranks among the leading global causes of preventable death and disabilities in human population. Understanding the sites of ethanol action that mediate its acute and chronic neural and behavioural effects is critical to develop appropriate treatment options for this disorder. The N-methyl-d-asparate (NMDA) receptors are ligand-gated heterotetrameric ion channels, which are known to directly interact with alcohol in a concentration-dependent manner. Yet, the exact molecular mechanisms and conformational dynamics of this interaction are not well understood. Here, we conducted a series of molecular dynamics simulations of the interaction of moderate ethanol concentrations with rat's wild-type GluN1–GluN2B NMDA Receptor under physiological conditions. The simulations suggest that glutamate or glycine alone induce an intermediate conformational state and point towards the transmembrane domain (TMD) as the site of action of ethanol molecules. Ethanol interacts by double hydrogen bonds with Trp635 and Phe638 at the transmembrane M3 helix of GluN2B. Alcohol not only reduces the pore radius of the ion channel within the TMD but also decreases accessibility of glutamate and glycine to the ligand-binding sites by altering the structure of the ligand-binding domain and significantly widening the receptor in that area. GRAPHICAL ABSTRACT


Introduction
Alcohol intake induces a multi-scale spectrum of spatiotemporal effects, the underlying neurobiology and physiology of which are not yet completely understood. concentrations of 5-50 mM that also produce intoxication, ethanol has been shown to inhibit the NMDAactivated ion current in a non-competitive manner [5][6][7][8][9][10] but also the NMDA-induced calcium influx, longterm potentiation and transmitter release [11][12][13][14]. The NMDA receptor is a ligand-and voltage-gated ion channel that is in general an assembly of GluN1, GluN2 (A-D) and GluN3 subunits forming heterodimers with the twofold symmetry axis running through the entire molecule composed of an amino terminal domain (ATD), a ligandbinding domain (LBD) and a transmembrane domain (TMD). The ion channel activates upon concurrent binding of glycine or d-serine and glutamate at the LBDs of GluN1 and GluN2 subunits, respectively [15] and relief of magnesium blockade of the channel pore by membrane depolarisation. A number of studies have focused on the identification of the molecular locus of alcohol action on NMDA receptors, particularly within the membraneassociated (M) domains of GluN1 and GluN2 subunits. Thereby, substitutions of phenylalanine (F637 and F639) at M3 of GluN1 subunit [16][17][18] as well as tyrosine T822, methionine M823 and alanine A825 of the GluN2 subunit [19,20] have been shown to influence the sensitivity of NMDA receptor to ethanol. Due to the high conservation level of the M-domains between GluN1 and GluN2 subunits, it is safe to assume that these sites correspond to ethanol sensitive amino acids in the other subunit. Despite the value of these investigations, their search strategies have been hypothesis-driven and site directed as a precise description of the complete structure of NMDA receptor including the M-domains was not available until recently [21,22]. Furthermore, the relationship between the ethanol interaction and the conformational dynamics of the receptor [23,24] is still missing.
Conformational dynamics is a multi-scale process and plays a critical role in the activation, deactivation and open-close activities of ion channels in living cells. In the present in silico study, we utilise molecular dynamics (MD) simulations to analyse the global binding properties of alcohol and conformational dynamics of the NMDA receptor based on the recently identified crystal structure of the intact heterotetrameric GluN1-GluN2B ion channel of rat at 4 Å [21]. While a number of recent studies have investigated the general binding properties of mostly homological models of this receptor [25][26][27][28], to our knowledge, this is the first investigation of ethanol-induced effects on a native NMDA receptor; thereby overcoming the intrinsic bias of site directed hypothesis-driven approaches and homology modelling. In total, eight configurations were simulated (100 ns/configuration), which represent all combinatorial ensembles of physiological concentrations of ethanol, glycine and glutamate in a close extracellular environment of the receptor. Here, we report a two-fold mechanism of NMDA inhibition by alcohol through a double hydrogen bond with tryptophan 635 and phenylalanine 638 located at the transmembrane M3 helix of the GluN2B subunit: (1) a widening of the pore radius at LBD and thereby decrease of accessibility of glutamate and glycine to their binding sites and (2) by tightening the TMD.

Methods
For MD simulations, wild-type crystal structure of GluN1-GluN2B NMDA receptor obtained by X-ray diffraction method was represented by Protein Data Bank ID 4PE5, in the presence of a GluN1 agonist, glycine, a GluN2 agonist, l-glutamate and an ATDbinding allosteric inhibitor, ifenprodil at 4 Å resolution [21]. The ligands were removed from the structure and the ROSETTA refinement program [29] was used to relax the structure, with protonation identical to wildtype. The protein was inserted into a pure 200 × 200 Å 2 palmitoyloleoylphosphatidylcholine (POPC) bilayer by deleting the overlapping lipid molecules, keeping 1070 POPC lipids, after which roughly 201,290 TIP3 water molecules were added into the simulation box of the size 7600 nm 3 . To neutralise the net charge and achieve physiological ion concentrations critical for the ion blockade of the NMDA receptors, 118 Na + and 50 Mg 2+ replaced the water molecules. The system preparation was done using VMD [30] and the Solvate and Ionize plugins. Fifty ethanol molecules were distributed randomly in the extracellular space using Packmol [31], reflecting a concentration of approximately 20 mM. Assuming instantaneous release of glutamate from vesicles and the glutamate concentration in the extracellular space as a function of time t and distance r from release, the glutamate concentration was calculated using the analytical solution of the two-dimensional diffusion equation c(r, t) = (N glu /4πhtD glu )e −r 2 /4tD glu . Hereby, h represents the width of the synaptic cleft ( ∼ 20 nm; [32]), N glu ( ∼ 7000; [33]) is the number of glutamate molecules in a vesicle and D glu ( ∼ 75 Å 2 /ns) is the diffusion coefficient of glutamate. The simulation of the diffusion equation suggests that within the active release zone (r < 2 nm) and after 1 ns, the simulation box will contain approximately 50 glutamate molecules represent a concentration of 160 μM. In order to achieve glycine binding at saturation level, the same number of molecules was distributed randomly in the extracellular space using Packmol. All simulations were performed using NAMD 2.9 [34] with CHARMM27 force field [35], a time-step of 2 fs, periodic boundary conditions and under constant NPT with pressure at 1 atm, temperature at 310 K. The Langevin dynamics and Nose-Hoover Langevin piston methods were used for temperature and pressure coupling, respectively. Particle mesh Ewald electrostatics was used with a 9 Å cutoff for non-bond interactions and neighbour lists updated every 10 steps. The system was energy minimised for 15,000 steps with conjugate gradient. To equilibrate the receptor with the membrane and the solvent, a 10-ns simulation was carried out. Simulation of the receptor in the closed state (without any constraints) was continued for another 100 ns. Thereafter seven independent targeted MD simulations (one ligand: ethanol, glycine and glutamate; two ligands: ethanol and glycine, ethanol and glutamate, and glycine and glutamate; three ligands: glycine, glutamate and ethanol) were carried out, each lasting 100 ns.

Structural analyses
In order to compare the affinities of the glutamate, glycine and ethanol molecules towards the NMDA receptor and identify their sites of action, average distribution densities were computed using the Volmap plug-in of VMD with a resolution of 1 Å, and averaged over the entire simulation period (100 ns) for ligand trajectories. It creates a map of the weighted atomic density at each grid-point. Thereby it replaces all atoms in the selection with normalised Gaussian distributions of width (standard deviation) equal to their atomic radius. The various Gaussians are then additively distributed on a grid. In this study, the weight of 'none' was chosen; therefore the Volmap results relate to number densities.
Following the identification of the sites of action of alcohol, the calculation of the diffusion coefficient of ethanol as well as the its interaction energy (Gibbs free energy) with the NMDA amino acids provide measures for the binding affinity with the amino acids. Furthermore, differences in the mobility of the agonists and alcohol may provide an explanation for the temporal framework of ethanol's modulating effects. In order to estimate the average diffusion coefficients, as a dynamical measure of the mobility of alcohol, glutamate and glycine molecules during the simulations, we have utilised the Diffusion Coefficient Tool Plugin for VMD (http://multiscalelab.org/utilities/ DiffusionCoefficientTool), which computes the threedimensional mean squared displacements (MSD)-based diffusion coefficients, based on the Einstein relation as D(τ ) = MSD(τ )/6τ . The mean square displacement is defined by MSD(τ ) = |r(τ ) − r(0)| 2 . To estimate the appropriate lag time τ for obtaining the diffusion coefficients, first the τ -average mean square displacement of ethanol molecules for an interval of 1000 frames (0.1 ns/frame) and step size of 10 frames was calculated as 19,110.01 Å 2 . Since the diffusion coefficient of ethanol in water is 84 Å 2 /ns, the appropriate lag time τ assumes 37.9 ns.
For an amino acid AA, the Gibbs free energy G ETOH AA was defined statically as the difference of the energy of the unbound V ub and bounded V b ethanol-AA complexes. The energy was given by the weighted sum of van der Waals, hydrogen bonds, electrostatic and desolvation potentials [36]: The parameters of this equation were obtained from the CHARMM27 force field. Although the Gibbs free energy G ETOH AA obtained with this strategy does not utilise the dynamical properties of the system, it provides estimations for the ranges of interaction energy between ethanol and NMDA amino acids.
The HOLE program [37] was used to calculate the pore radii along the pore axis for each snapshot as a measure for conformational dynamics of the closed-state and activation simulations (sampled once per 50 ps). Averages over 500 snapshots from the last 5 ns of the closed-state simulation and 1000 snapshots from the last 10 ns of activation simulations were taken as values for the closed, ligand-induced intermediate and open states, respectively. Furthermore, root mean square deviation (RMSD) of backbone amino acids were calculated, presented as heat maps and compared among different conditions to characterise the conformational dynamics at amino acid level. 3Vee [38] were used to calculate pore size (i.e. diameter and accessible volume) for all cases for temporal snapshots at 25, 50, 75 and 100 ns.

Statistical analysis
One-way ANOVA was used for analysis of the raw and normalised data. Whenever significant differences were found, post hoc Student Newman Keul's tests were performed. The chosen level of significance was p < .05.

Diffusion coefficients
The τ -average mean square displacements of glycine and glutamate were 19156.42 and 17033.83 Å 2 , respectively. Using a lag time of 37.9 ns, the respective average diffusion coefficients were determined as D gly = 84.2 Å 2 /ns and D glu = 74.9 Å 2 /ns. Glycine and ethanol did not differ effectively in their diffusion properties during the simulation, while the mobility of both molecules was higher to a level of 112% in comparison to glutamate.

Distribution density of alcohol
In agreement with previous investigations [39], the simulations suggest a high distribution density of glutamate and glycine molecules at LBDs of GluN2 and GluN1 subunits, respectively. In particular, glutamate tend towards a close proximity of residues [478; 485] (such as Pro478, Thr480, Arg485), as well as Ser654 and Thr655 at the GluN2B subunit. In contrast, glycine molecules are rather attracted toward the residues located at [516; 523] (such as Pro516, Thr518, Arg523) and Ser688 and Val689 at GluN1 subunit.
The analysis of the average distribution density of ethanol molecules during the 100 ns simulations suggests that residues that participate in alcohol interaction are located within the TMD (Figure 1(A)) at M3 positions Trp635 and Phe638 of the GluN2B subunit ( Figure  1(B)), which may interact directly with the second binding site for glutamate molecules. Thereby, the Gibbs free energy calculated as the difference of interaction energies of the unbounded and bounded ethanol-Trp635 and ethanol-Phe638 takes values −2.9 and 2.94 kcal/mol (Figure 1(C)), respectively.

Conformational dynamics
Since NMDA receptors form ion channels, the analysis of substance-induced alterations in their pore geometry may provide insights on the functioning of the receptors. Using HOLE and 3Vee, we have estimated the changes in pore radius and volume of the receptor during simulation due to the action of the co-agonists as well as ethanol ( Figure 2 and Table 1). In agreement with experimental studies, the results suggest that single binding of glutamate or glycine induces an intermediate pore size enhancement within the TMDs. While a full transformation of the receptor into an open conformational state cannot be observed in the time-scales of our simulation (see also [25][26][27]), our simulations demonstrate the affinity of both co-agonists to their known binding sites as well as enhancements in pore volume that are significantly larger than the volume changes induced by a single agonist.
The findings further suggest a two-fold mechanism of action of alcohol that may explain the inhibition of the receptor activity: (1) alcohol decreases the accessibility of glutamate and probably glycine to their ligandbinding sites by altering the structure of LBD and widening the receptor in this area. In particular, due to the higher mobility of alcohol molecules in comparison to glutamate, ethanol-induced conformational changes at GluN2B may propagate to the nearby ligand-binding site for glutamate and thereby exacerbate the binding process of glutamate to the receptor; (2) alcohol decreases the   pore radius in TMD significantly and thereby may reduce the permeability of the receptor for calcium ions. In addition, the simulations suggest that the concurrent action of both co-agonists may stabilise the ATD (Figure 3(D)) and consequently improve ion selectivity and permeation through the channel [40]. Since single exposure to ethanol (Figure 3(A)), glycine (Figure 3(B)) or glutamate (Figure 3(C)) does not appear to have a similar effect of the receptor. In fact, the interaction of a single co-agonist or ethanol increases the pore radius at ATD particularly unilateraly at α5 and α6 helices and 8 β-sheet of GluN2B subunit.
This observation may further support the hypothesis on the involvement of the ATD in optimal channel gating and regulation of agonist potency. All-residue RMSD calculation within the GluN2B segment containing α5 and α6 represented as heat maps. In agreement with the results on the ligand-induced pore-radius alterations, ethanol (A) most strongly alters this pore-forming component of the ATD. The co-agonists glycine (B) and glutamate, similarly but less dominantly influence the conformation of this segment throughout the simulation. The concurrent interaction of glutamate and glycine in turn, stabilises the ATD and induces only a temporally limited response, vanishing towards the end of simulation period.

Discussion
In this study, we used, to our knowledge for the first time, MD simulations unbiased by homological modelling to characterise the ethanol-induced conformational dynamics of the crystal structure GluN1-GluN2B NMDA receptor. The results suggest the appropriateness of MD simulations to pursue this question, as they are in agreement with previous experimental studies identifying sites of action of glutamate and glycine at these receptors. Our investigation emphasises the role of phenylalanine 638 in the ethanol-induced inhibition of the NMDA receptor. Furthermore, it suggests the potential involvement of tryptophan 635, which yet requires experimental validation.
The simulations associate the hydrogen bonding of ethanol to the above-mentioned amino acids with complex conformational changes throughout the different domains of the receptor. Firstly, ethanol as well as glutamate and glycine tend to interact with the amino acid terminal and increase the pore radius unilateraly. This effect is only controlled by the presence of both coagonists. Secondly, ethanol increases the pore radius at the LBD through a propagation of the conformational changes from its binding sites to the binding site of glutamate. We hypothesise that the widening of the receptor in this area may decrease the accessibility of glutamate to the receptor, reduces its binding probability and thereby contributes to the well-known inhibition of NMDA receptor activity. And thirdly, ethanol induces a reduction of the pore radius in the TMD, which leads to a reduction of ion permeability through the channel.
Our findings may improve current understanding of ethanol-induced conformational dynamics of NMDA receptor, which is critical for development of efficient drugs to treat alcohol-related disorders. However, despite the advantages of our approach, it bears natural and intrinsic limitations. The MD simulations are performed at time-scales that may not capture the entire physiologically relevant dynamics, particularly they are blind to dissociation kinetics of ligands. Furthermore, we are unable to simulate the extracellular matrix, as well as active and passive intracellular processes that may influence the receptor function. In addition, the finite and limited size of the simulation box and the application of periodic boundary conditions may lead to unrealistic processes as the large ligands may appear within the intracellular space and interact with the receptor components as well as the membrane without diffusing through the membrane.
Based on these limitations, the simulation results should be considered with special care and experiments such as site directed mutagenesis combined with singlechannel recordings or FRET technology [24] are necessary to validate the simulation results on ethanol-induced conformational dynamics.