Computational Investigation of BMAA and Its Carbamate Adducts as Potential GluR2 Modulators

Beta-N-methylamino-l-alanine (BMAA) is a potential neurotoxic nonprotein amino acid, which can reach the human body through the food chain. When BMAA interacts with bicarbonate in the human body, carbamate adducts are produced, which share a high structural similarity with the neurotransmitter glutamate. It is believed that BMAA and its l-carbamate adducts bind in the glutamate binding site of ionotropic glutamate receptor 2 (GluR2). Chronic exposure to BMAA and its adducts could cause neurological illness such as neurodegenerative diseases. However, the mechanism of BMAA action and its carbamate adducts bound to GluR2 has not yet been elucidated. Here, we investigate the binding modes and the affinity of BMAA and its carbamate adducts to GluR2 in comparison to the natural agonist, glutamate, to understand whether these can act as GluR2 modulators. Initially, we perform molecular dynamics simulations of BMAA and its carbamate adducts bound to GluR2 to examine the stability of the ligands in the S1/S2 ligand-binding core of the receptor. In addition, we utilize alchemical free energy calculations to compute the difference in the free energy of binding of the beta-carbamate adduct of BMAA to GluR2 compared to that of glutamate. Our findings indicate that carbamate adducts of BMAA and glutamate remain stable in the binding site of the GluR2 compared to BMAA. Additionally, alchemical free energy results reveal that glutamate and the beta-carbamate adduct of BMAA have comparable binding affinity to the GluR2. These results provide a rationale that BMAA carbamate adducts may be, in fact, the modulators of GluR2 and not BMAA itself.


Calculation of BMAA analog similarity with glutamate
To calculate whether alpha-, beta-carbamate and BMAA are structurally similar with respect to the native agonist glutamate, we used the Tanimoto coefficient 2 (Tc), a similarity metric for comparing chemical structures using chemical fingerprints.A chemical fingerprint is a series of bits that represent the presence or absence of chemical substructures, patterns, physicochemical characteristics, and other properties in a molecule.If the property is present in the molecule, then the fingerprint is 1 (ON) and if the property is absent, then the fingerprint is 0 (OFF).The Tc to compare the similarity of two molecules is given by Equation 1: where Nc is the number of common ON bits (1) in both molecules, NA is the number of ON bits in molecule 1 and NB is the number of ON bits in molecule 2. The Tc was calculated using ChemBioServer 3 .The results are presented in Table S1.In addition, to compute the 3D similarities among alpha-, beta-carbamate and BMAA with respect to the native agonist glutamate, we used Shape Screening module (Maestro 2020.3,Schrödinger, Inc).3D similarity screening can be applied to match the overall shape of different compounds 4 by performing initial alignments of atom triplets and then conducting refinement and overlap scoring. 5To compute the volume overlap scoring between the molecules, we used the Pharmacophore type approach, in which every compound is treated as a collection of pharmacophore sites defined by the Phase 6,7 pharmacophore feature and then calculates the overlapping volumes between same feature sites. 8A shape similarity score of 0 denotes the maximum dissimilarity between the binding poses of two compounds, while a shape similarity score of 1 means that two compounds have the maximum 3D similarity.The results are presented in Table S2.

Glutamate-GluR2 model construction
The atomistic model of the GluR2 complex was created based on the crystal structure with PDB ID 1FTJ 1 , which comprises the GluR2 ligand binding domain (LBD) in complex with glutamate.
GluR2 consists of four subunits each one containing the natural agonist glutamate bound to LBD.
However, the crystal structure contains only three chains, where chain A includes one zinc ion and chains B and C include two zinc ions each.Zinc ions are not located inside the binding pocket and are at a distance of ca. 17 Å from the center of mass of glutamate.In this study, only chain A was used for the MD simulations/alchemical free energy calculations to reduce the computational effort.Chain A in complex with glutamate was prepared using the Protein Preparation Wizard 9 (Schrödinger 2020-3).Missing residues of chain A were added and refined using Prime 10,11 , tautomer/ionization states were assigned and the resulting structure was minimized.

Relative binding free energy calculations
The relative binding free energy for the glutamate/beta-carbamate perturbation was estimated using alchemical free energy calculations 12 .The free energy of difference between two end states in both complex and solvent phases could, in theory, be calculated by sampling configurations, computing the energies of each microstate, and evaluating the partition function Z as written in Equation 2: where k B is the Boltzmann constant, T is the temperature of the system and  & ° is the binding constant.The computation of the binding free energy depends on solving the different partition functions of the system, for example the  ()*+,-.is given by Equation 3: where  ()*+,-.is the potential energy of conformation q of the two end states in the complex phase.However, direct computation of these partition functions is computationally intractable due to the amount of sampling needed to compute absolute free energies.Instead, it is much easier to compute the differences between the binding affinities of two structural similar compounds using a relative free energy calculations formula derived from Equation 2: This approach depends only on the transformation of system A to system B in the complex and solvent phases and thus it resolves errors associated with computing absolute free energies.This transformation is typically performed by connecting the two systems using a coupling parameter λ, such that the potential energy function U(λ) interpolates between the initial state A (λ=0.0) and the final state B (λ=1.0).
To obtain accurate relative binding free energy predictions, a sufficient degree of overlap of the phase space between the initial (A) and the final (B) state is required.This is feasible by introducing a series of N alchemical intermediate λ windows along the perturbation path.When sufficient overlap is not attained, one can either rerun the calculations using more λ windows or insert intermediate compounds to bridge the conformational space between the two states.
During a relative binding free energy calculation involving multiple compounds, one ligand from the congeneric series acts as a reference to which all other molecules are aligned.As a result, the conserved binding mode as well as the improved overlap between the windows are secured.The reference compound is usually selected as the one that is representative of the ligand set and at the same time with the higher certainty of the predicted binding mode.Then, based on the reference ligand, the compound set forms a network, where the "edges" are the transformations performed using alchemical free energy simulations.The relative binding affinities calculated from the network may contain systematic errors due to inability of the force-field to describe the molecular interactions or/and motions or due to unconverged simulations.Here, an estimation of these errors was obtained using the cycle closure error method as the relative free energy should be zero in a closed cycle.

Phase space overlap between neighboring λ windows from NAMD/OPLS-AA/FEP and AMBER/ff14SB/TI protocols
The results obtained from the NAMD/OPLS-AA/FEP and AMBER/ff14SB/TI protocols concerning the phase space overlap between the neighboring λ windows for each perturbation are listed in Tables S3 -S7: Table S3.Phase space overlap between neighboring λ windows of perturbations glutamate→intermediate-        replica1 of the AMBER20/ff14SB/TI protocol, C) replica2 of the AMBER20/ff14SB/TI protocol, and D) replica3 of the AMBER20/ff14SB/TI protocol.

Figure S5 .
Figure S5.Perturbation network used for the initial RBFE calculations.Glutamate and beta-

Figure S6 .
Figure S6.Histograms of the probability distribution of the forward, Pfwd(ΔU) (black solid line)