Analysis of small molecules against the NMDA receptor: an insight from virtual screening and molecular dynamics simulation based ndings

Alzheimer’s disease (AD) is a chronic intensifying neurodegenerative disorder and accounts for three-fourths of dementia cases. To date, there is no effective treatment available which can completely cure AD. The available medications can slower AD progression and can provide symptomatic relaxation. The N-methyl-d-aspartate receptor (NMDAR) plays a paramount role in the survival of neurons and synaptic plasticity. Although, excessive function of NMDAR cause excitotoxicity. Due to this the cell death process activated resulting into neurodegeneration and promotes AD. Hence in this study, we have screened 98,072 natural compounds by using Smina and idock. After that 154 compounds were selected and ADMET is predicted by using the pkCSM web-based server. From the ADMET analysis, 18 compounds were chosen and employed for the re-docking studies by using Autodock Vina. Then from the docking result, we have selected top three complexes (NMDAR-ZINC4258884, NMDAR-ZINC8635472, and NMDAR-ZINC15675934) and employed them for the 100 ns MDS studies. Based on MDS result analysis we have concluded that NMDAR-ZINC4258884 and NMDAR-ZINC15675934 are the best stable complex and can function as a lead compound against the NMDAR. Although this is a theoretical study while we have shortlisted only two compounds out of 98072 compounds and proposed them to the scientic community worldwide for further experimental validations.


Introduction
Alzheimer's disease (AD) has been outlined as a chronic intensifying neurodegenerative disease. Its onset is slow, but it worsens gradually over time [1]. It has been believed to be a cause of 60-70% of dementia. The foremost common early symptom is claimed to be, "di culty in remembering recent events". With the advancement in illness, symptoms include Aphasia, self-negligence, and performance problems. In line with a report, in the last ten years, casualties from cardiovascular diseases have shown a 7.8 percent decrease whereas casualties from AD are hyperbolic145 percent [2]. As per an estimate, 6.2 million Americans aged 65 + are diagnosed with AD in the year 2020 [2]. In an estimate, it is predicted that in 30 years from now, dementia (showing 60-70% contribution for AD) will be affecting approximately 152 million people worldwide. An estimated four million Indians are affected by one or another form of dementia [3]. The AD is mainly characterized by the β-amyloid plaques and neuro brillary tangles (NFTs).
One major cause of AD is an accumulation of β-amyloid protein which is (in normal condition) the metabolic waste product present in the uid between brain cells [4]. The NFTs are made by the tau protein hyperphosphorylation via the Glycogen Synthase Kinase 3β (GSK3β) and Cyclin dependent kinase 5 (CDK5) enzymes. The NFTs are present inside the neuronal cell and promotes the cell death. We have proposed various compounds to reduce the tau phosphorylation induced by GSK3β [5,6] and CDK5 [7,8] recently. In AD, β-amyloid comes to create amyloid plaques which are thought to instigate neuroin ammation and disrupt the information exchange among neurons [9]. Now, in the neuron, NMDAR which are a type of ionotropic glutamate receptors (iGluR) whose function involves the mediation of the excitatory transmission in the brain is present [10]. The protagonist molecule NMDA (N-methyl-Daspartate) binds selectively to the NMDA receptors only, hence it is named as NMDA receptor. NMDA receptors areactivatedupon the binding of theglutamate and glycine (or D serine) with these receptors; upon activation they allow the ow of the positively charged ions throughthem. They have an important role in forming memory, learning, and synaptic plasticity ( Fig. 1) [11].
NMDAR is thought to be different from other ionotropic glutamate receptors as it has voltage-dependent activation through Mg 2+ blockade removal, high Ca 2+ permeability and comparably slow ligand-gated kinetics. In normal condition, the resting membrane potential is -70mV, at this potential the NMDAR Ca 2+ channel is thought to be blocked by Mg 2+ [12]. In the Long-Term Potential (LTP), there is a strong and long-lasting release of glutamate from the presynaptic terminal, which activates α-amino-3-hydroxy-5methyl-4-isoxazolepropionic acid receptor (AMPAR) and subsequent withdrawals of glutamate, which ultimately removes the Mg 2+ inhibition of NMDAR and thus allows the entry of Ca 2+ , which ultimately leads to improved synaptic strength [13]. Studies suggest that β-amyloid proteins which accumulate in the brains of AD patients can cause abnormal increase in synaptic glutamate levels by blocking glutamate uptake or making glutamate free from glial cells [14]. The binding of glutamate to the NMDA receptor triggers the extracellular Ca 2+ ight that regulates membrane permeability and synaptic transmission [14]. When glutamate levels increase abnormally, in addition to reactivation of NMDA receptors leading to an overgrowth of Ca 2+ substances that remove Mg 2+ inhibition, ultimately leading to cell division and death [15]. More in ux of Ca 2+ occurs since NMDARs increased the penetration of Ca 2+ compared to other iGluR [16]. Therefore, blocking this receptors can be a potential treatment for AD, thus stopping the excessive in ux of Ca 2+ [17][18][19]. The memantine is an FDA approved antagonist which blocks the NMDAR activity and gives the symptomatic relax in the case of AD. Although it has several side effects, hence we are proposing new compounds through this study and analysis.
We have used the structure-based virtual screening approach and 98,072 compounds were screened against the NMDAR. After that 154 compounds were chosen for the ADMET (Absorption, Distribution, Metabolism, Excretion, and Toxicity) analysis and then selected 18 compounds were used for the redocking studies by using the AutoDock Vina. Finally, 3 compounds were selected and used for the 100 ns simulation. Lastly, based on MDS results we have proposed that NMDAR-ZINC4258884 and NMDAR-ZINC15675934complexes are showing stability and these compounds can function as novel as well as potential compounds against the NMDAR. A comprehensive methodology is shown in Fig. 2. 2. Methodology

Preparation and retrieval of Ligand and Protein
The RCSB protein data bank (PDB) was used to download the structure of the NMDA receptor protein (PDB ID:1PBQ, 1.90 Å, X-ray) [20]. The DK1 (5,7-Dichloro-4-Hydroxyquinoline-2-Carboxylic Acid) is an experimentally proved inhibitor and co-crystallized with this structure. The preparation of structure is done by using Chimera 1.13.1 software, San Francisco, CA [21]. Then, by using Amber ff99SB force eld [22] the structure was subjected to minimization using Chimera 1.13.1. After the minimization step, further protein structure is used in AutoDock Tools to convert from PDB to .pdbqt. The ZINC database [23] is used to download the structure of 98,072 compoundsin .mol2 format. These ligands were converted from mol2 to .pdbqt le format for the screening by using a Python script (http://autodock.scripps.edu/faqs-help/how-to/how-to-prepare-a-ligand-le-for-autodock4).

Analysis of binding site for 1PBQ
The NMDA receptor protein is co-crystallized with its known inhibitor DK1. Therefore, the binding cavity related to DK1 is selected for virtual screening. We check for the DK1 associated cavity by analysing the DK1 interactions with NMDAR. For the grid preparation of 1PBQ, numerous residues like Gln13, Asp224, Phe92, and Ser180, etc. were selected.

Virtual Screening
Virtual screening is an e cient process to nd a potent inhibitor using computational approaches [24,25].
It is an e cient and less time-consuming process where we can screen a large compound dataset within few days. Here we have used Smina -A fork of AutoDock Vina and idock [26] for the virtual screening. They both are inspired from the Autodock Vina [27] method but they are less time-consuming than Autodock Vina. The idock is a multithreaded and improved algorithm for virtual screening and very fast compared to all other available tools. It can screen a compound within seconds. The user can de ne the computer threads and based on a de ned core it can screen thousands of compounds within hours. Hence we have used these software for the virtual screening and then 154 compounds which were showing the >= -11.5 Kcal/mol binding a nity were adopted selected for the ADMET analysis.

ADMET prediction
The ADMET prediction was done of the 154 ZINC compounds by using the pkCSM tool (http://biosig.unimelb.edu.au/pkcsm/prediction), which is an ADMET predicting web-based server [28]. This software is freely available and to perform the prediction user can provide the ligand with the SMILES string or draw the chemical structure on the screen. This web server contains a large dataset of compounds approved by the FDA. Also, the dataset is available from scienti c literature present in PubMed (https://pubmed.ncbi.nlm.nih.gov/) and Google Scholar (https://scholar.google.com/). The ADMET descriptors in absorption are CaCo2 permeability, Intestinal absorption (human), water-solubility, p-glycoprotein substrate, p-glycoprotein 1 and 2 inhibitors, skin permeability; in Distribution are Volume of distribution (VDss) Human, Blood brain barrier (BBB) permeability, Fraction unbound (Human), CNS permeability; in Metabolism are Cytochrome P450 inhibitors, CYP2D6/CYP3A4 substrate; in Excretion are Renal OCT2 substrate, Total Clearance; in Toxicity are Rat LD50, T. Pyriformis toxicity, AMES Toxicity, Minnow Toxicity, Maximum Tolerated Dose, Hepatotoxicity, Oral Rat Chronic Toxicity, Skin Sensitization, hERG 1 and 2 inhibitors etc. We enumerated all the above stipulations for all these 154 compounds using the ADMET tool and choose the best out of these, which are employed for molecular docking analysis.

Molecular Docking Simulation
By employing ADMET prediction, 18 compounds were picked and further re-docking was employed to them through the software AutoDock Vina [27] with the control ligand DK1. Prior to the docking the preparation of the protein and ligands were performed by using the AutDock Tools [29]. The protein was prepared by adding all the hydrogens, Kollman charges and by assigning the AD4 type radii. After that the protein is saved in to the pdbqt le format. In the ligand preparation, all the hydrogens and Gasteiger charges were added. After this the ligands were saved as a .pdbqt le. The setting of the grid box was: Center_X = 3.909, Center_Y = 38.484, Center_Z= -18.612, numbers of points were set as, X_dimension = 40, Y_dimension = 40, Z_dimension = 40 with 0.419 Å spacing. After this step the docking were carried out to get the binding a nity and binding pose using Autodock Vina. For each ligand, 8 binding poses were produced and the analysis is done on the basis of least binding energy and the best binding pose.
AutoDock Vina is more e cient when compared to AutoDock because of its rapid docking of ligands and more accurate results [27].
2.6 Molecular dynamics simulation GROMACS 2018.2 suite [30] was used to perform molecular dynamics (MD) simulations on the apo-NMDAR, NMDAR-DK1, NMDAR-ZINC4258884, NMDAR-ZINC8635472 and NMDAR-ZINC15675934 complexes. A rhombic dodecahedral box having a distance of 10 Å from the closest edge was used to solvate the complexes, thereby, adding more than 19,734 water molecules. Also, for the purpose of neutralising the system, 4 Cl − counter ions were included. For the ligand GAFF force eld was used; for protein, Amber ff99SB-ILDN force eld was used and TIP3P model was used for water. Using an isobaricisothermal ensemble, production runs were performed for 1 ns, after undergoing a standard preparation routine including energy minimization, annealing, and equilibration. Other simulation conditions like the parameters for the barostat/thermostat, the use of constraints, the modelling of electrostatic and nonelectrostatic interactions was also decided beforehand. Undergoing subsequent analysis, after the elimination of the protein rot translation by least-squares t with respect to the Cα atoms. Various other analyses have also been carried out to predict the protein-ligand stability.

Virtual screening analysis
The virtual screening has been carried out by using Smina-A fork of AutoDock Vina and idock software.
We have screened 98,072 compounds against the NMDAR and binding a nity of all these compounds were displayed in Supplementary Table S1. The ZINC04277685 was the highest binding energy compound in Smina virtual screening result with − 13.1 Kcal/mol binding a nity while idock showed ZINC04258868 as the top compound having the binding a nity value as -13.29 Kcal/mol. The difference between these top compounds prediction could be the result of different tools while we have taken the consensus result from both the software. The ZINC39410302 is predicted as the lowest binding a nity compound by Smina virtual screening with − 3.9 Kcal/mol binding a nity while idock predicted ZINC01690436 as the lowest binding a nity compound with − 3.99 Kcal/mol. After that, we have taken 154 compounds which are showing >= -11.5 Kcal/mol binding a nity in both the tools and used for the ADMET prediction.

ADMET Descriptors Analysis
ADMET stands for Absorption, Distribution, Metabolism, Excretion, and Toxicity, these features are essential in support of drug designing studies and the drug cannot be validated to enter into the market until it ful ls the ADMET requirements. Therefore, this step is of utmost important. We employed the selected 154 compounds in the ADMET analysis via the pkCSM server. As we are targeting the Central Nervous System, the BBB is the most important parameter here, 31 compounds were rejected based on BBB. The second priority was given to CNS permeability and 12 compounds were rejected according to this. The next priority was given to the absorption of the drug, as it is an important factor for oral drug discovery. In the absorption category, HIA is evaluated, as it tells us if the drug is absorbed in the intestine or not, based on this, no compound was rejected. Caco-2 cell permeability is also checked to see the drug assimilation in the large intestine, 94 compounds were rejected based on this. Next, in absorption, it is checked that whether a given compound is likely to be a substrate of p-glycoprotein (Pgp) or not. Pgp usually functions as a biological barrier by excluding xenobiotics and toxins out from the cells, 2 compounds were rejected in this step (Supplementary Table S2 and S3).
Next, we check for the Cytochrome P450 barriers (CYP450). Here it is assessed that for a given isoform, the potential drug molecule is likely to be a CYP450 inhibitor or not. This enzyme is present in the liver and is responsible for the detoxi cation of the xenobiotics. Hence, it can deactivate many drugs. Here we look for the CYP2D6/CYP3A4 substrate, which lets us know that whether the drug can be metabolized so that it can be removed from the body. In the case of CYP2D6, we got only 11 positive results whereas in the case of CYP3A4 we got 151 positive results and only3 compounds were rejected (Supplementary   Table S4).
Another important criterion for drug discovery is to assess the toxicity of the potential drug molecules. Even if the e cacy of the drug is high, it cannot be launched in the market if it is toxic, we checked for the AMES toxicity which tells us about the possibility of carcinogenicity, based on this, 66 compounds were rejected. hERG is an ion channel known for the electrical activity of the heart. hERG codes for the potassium channels. Inhibition of these channels leads to long QT-syndrome which leads to the fatal ventricular arrhythmia. Therefore, we have to look for the drug candidates that do not inhibit hERG channels, keeping this as the basis, further 2 compounds were rejected (Supplementary Table S5). Considering all these parameters, a total of 18 potential drug targets were selected from the 154 compounds.

Molecular Docking Simulation Analysis
The 18 compounds chosen via ADMET analysis along with the control compound DK1 were submitted to molecular docking through AutoDock Vina software. The control compound is redocked using the AutoDock Vina, giving the binding a nity to be as -8.6 Kcal.mol − 1 . The compound ZINC705167 illustrated the highest binding a nity of -12.7 Kcal.mol − 1 and the compound ZINC705168 illustrated the lowest binding energy of -8.7 Kcal.mol − 1 , followed by the least binding a nity of the control compound of -8.6 Kcal.mol − 1 . The binding energy, hydrogen bonds, interacting residues from AutoDock Vina, Smina and idock are tabulated in the Supplementary Table S6 for all the chosen selected ligands along with the DK1. The respective ZINC IDs, binding a nities, 2D chemical structure for the top 3 selected compounds with DK1 are shown in Table 1.

DK1
The DK1 exhibited the binding a nity of -8.6, -8.6 and − 8.71 Kcal.mol − 1 by using Autodock Vina, Smina, and idock respectively. The DK1 showed less binding a nity in comparison to all the other ligands which represents that predicted hits are good as compare to DK1. We have seen 2, 4 and 5 hydrogen bonds from Autodock Vina, Smina, and idock respectively. We have seen many common residues such as Pro124 and Thr126 etc. from all three software's. The residue interaction diagram is shown in Fig. 3A

ZINC4258884
The ZINC4258884 exhibited the binding a nity of -12.5, -12.5, and − 12.77 Kcal.mol − 1 by using Autodock Vina, Smina, and idock respectively. We have analysed the Autodock Vina docking complex which is showing one hydrogen bond with the Ser180 and various other residues are involved in the hydrophobic interaction. Various common residues are also found in all three docking software. The residue interaction diagram is shown in Fig. 3B. The detailed residues are shown in the Supplementary Table S6.

ZINC8635472
The ZINC8635472 exhibited the binding a nity of -12.1, -12.1, and − 12.18 Kcal.mol − 1 by using Autodock Vina, Smina, and idock. Using AutoVina software, the complex forms 1 hydrogen bond with Thr126. The complex is also stabilized by many hydrophobic interactions. Common residues are also found in all three docking software. The residue interaction diagram is shown in Fig. 3C

ZINC15675934
The ZINC15675934 is also docked by all the software against the NMDAR. It showed − 11.9, -12.2, and − 12.3 Kcal.mol − 1 from Autodock Vina, Smina, and idock software respectively. We have not seen any hydrogen bonds between the NMDAR and ZINC15675934 through the AutodockVina while the complex is seen to be stabilized by many other interactions like hydrophobic interactions etc. The residue interaction diagram is shown in Fig. 3D. The detailed residue interaction is shown in the Supplementary Table S6.

Molecular dynamics simulation
The Molecular dynamics simulation is a widely used technique to evaluate the docked stability of proteinligand complex. Therefore this study comprises of using MDS for the validation of docking complexes. We have selected 3 top complexes with the apo-NMDAR and NMDAR-DK1 for the 100 ns simulation. Various structural parameters such as root mean square deviation (RMSD), root mean square uctuation (RMSF), Radius of gyration (Rg), Solvent accessible surface area (SASA), number of hydrogen bonds and Principal component analysis (PCA) were carried out.

Stability analysis
For the purpose of prediction of the stability of simulation, the RMSD has been carried out. We have computed the RMSD value for 100 ns and plotted it in Fig. 4A. It represents the deviation from the initial structure to the next structure. The Fig. 4A represents that after the 50 ns all the simulations got the stability and can be used for the analysis. The average RMSD for apo-NMDAR, NMDAR-DK1, NMDAR-ZINC4258884, NMDAR-ZINC8635472 and NMDAR-ZINC15675934 was 0.54, 0.53, 0.27, 0.48 and 0.45 nm respectively. The average RMSD that after ligand binding all the complexes got the stability as compare to apo-NMDAR and control ligand DK1. The predicted hits also showed less RMSD value in comparison to the DK1 and representing a well-stable complex. Hence from here, we have predicted that all our complexes got the equilibrations and the simulation trajectories are producing accurate results. Hence we considered last 50 ns trajectories for further analysis.

Residue mobility analysis
The residue level mobility analysis by using the RMSF analysis was carried out. This is a very important analysis in the case of NMDAR because it is an ion channel and very exible in nature. Hence it will show the higher uctuation in several residues. We have plotted the residue mobility in the Fig. 4B. In the Fig. 4B, we have seen very high uctuation between the residues 94-106. The plot showed higher uctuation in the maximum region of the protein as compared to the general uctuation due to the ion channel nature of the NMDAR. The average RMSF for apo-NMDAR, NMDAR-DK1, NMDAR-ZINC4258884, NMDAR-ZINC8635472 and NMDAR-ZINC15675934 was 0.17, 0.13, 0.10, 0.18 and 0.12 nm respectively. The highest uctuation is showed by NMDAR-ZINC8635472 whilst in comparison to the apo-NMDAR and NMDAR-DK1, NMDAR-ZINC4258884 and NMDAR-ZINC15675934 were seen to show a lesser uctuation. It displayed the stability of NMDAR-ZINC4258884 and NMDAR-ZINC15675934 and showed their potential use as a lead compound against NMDAR.

Compactness analysis
The compactness of the complexes has been analyzed by using the Rg analysis. The Rg has been measured from the radius of the protein and represents the compactness of the protein. Hence here also we have calculated the Rg for the last 50 ns equilibrated trajectory. It is plotted in Fig. 5A. In Fig. 5A we can see that NMDAR-ZINC4258884 showed very least value as compare to all other complexes. Other complexes and apo-NMDAR showed a similar type of Rg value. The average Rg value for apo-NMDAR, NMDAR-DK1, NMDAR-ZINC4258884, NMDAR-ZINC8635472, and NMDAR-ZINC15675934 was 2.16, 2.16, 2.06, 2.19, and 2.18 nm respectively. Here also the NMDAR-ZINC8635472 showed the highest value as compare to all others while NMDAR-ZINC4258884 and NMDAR-ZINC15675934 showed the good Rg value and showing less uctuation. Hence from here also we can say that NMDAR-ZINC4258884 and NMDAR-ZINC15675934 are stable complexes in the respect of Rg analysis.

Solvent accessible surface area analysis
The SASA has been carried out for the prediction ofthe ligand-induced solvent-accessible area changes.
We worked upon usingthe last 50 ns trajectory and calculated the SASA value and plotted it in Fig. 5B. The Fig. 5B

Hydrogen bonds analysis
After that, the number of hydrogen bonds for the last 50 ns was computed and was plotted in Fig. 5C. The hydrogen bonds are crucial for ligand-protein stability. From the Fig. 5C, we can see that NMDAR-DK1 exhibited highest number of hydrogen bonds in comparison to all the estimated complexes. The NMDAR-DK1 showed the average 4-6 hydrogen bonds throughout the simulation. After that NMDAR-ZINC8635472 showed more hydrogens bonds such as 1-4 for all simulations. The NMDAR-ZINC4258884 and NMDAR-ZINC15675934 showed 0-2 hydrogen bonds throughout the simulation. The hydrogen bond analysis represents the stability of all the complexes and more number of hydrogen bonds is observed in the DK1.

Principal Component analysis
The PCA analysis has been carried out for predicting the correlated motions ligand binding. Here, the calculation of the eigenvalue vs. and eigenvector was performed. Due to the clear depiction of the result only the rst 40 eigenvectors were selected for the analysis and were plotted in Fig. 6A.
In the Fig. 6A  From here we have seen that the overall dynamics of the protein is accountable for the rst few eigenvectors. Hence, the rst two eigenvectors were selected and plotted against each other in phase space and shown in Fig. 6B. The Fig. 6B represents that NMDAR-ZINC15675934 and NMDAR-ZINC4258884 is the most stable cluster in comparison to all other complexes. These complexes showed very stable cluster which represents the complex stability. After these two the NMDAR-DK1 is showing the stable cluster. We have seen a very wide and dispersed cluster for the NMDAR-ZINC8635472 and apo-NMDAR. The 2D projection report also represents that ZINC15675934 and NMDAR-ZINC4258884 are the stable complex as compare to the NMDAR-ZINC8635472.
Lastly, we have predicted the residue-wise correlated motions from the last 50 ns trajectory using only one eigenvector and plotted in Fig. 6C. Here we have seen a similar pattern as like RMSF analysis. We have seen a high peak from residues 87-109. The overall eigRMSF is also high for all the residues. The apo-NMDAR, NMDAR-DK1, NMDAR-ZINC4258884, NMDAR-ZINC8635472 and NMDAR-ZINC15675934 showing average eigRMSF is 0.11, 0.07, 0.04, 0.11 and 0.06 nm respectively. Here it was observed that average eigRMSF for apo-NMDAR and NMDAR-ZINC8635472 is high as compare to other complexes. The NMDAR-DK1 is also showing less eigRMSF value while our predicted hits NMDAR-ZINC4258884 and NMDAR-ZINC15675934 are showing very less eigRMSF value and representing the well stable complexes.
The overall PCA result also agrees with the above analysis namely RMSD, RMSF, hydrogen bond analysis and SASA. These all result indicates that out of three complexes the NMDAR-ZINC4258884 and NMDAR-ZINC15675934 are stable and has a potential to act as lead compounds.

Conclusion
AD is a progressive neurological disorder mainly affecting old age people. There is no medication available to cure AD. It is a multifactorial disease and several proteins simultaneously activate and spread the disease. The NMDAR is a key target in the case of AD hence we have targeted the NMDAR with the natural compounds. We started our study with 98,072 compounds and we have picked out 154 compounds from the virtual screening and performed the ADMET analysis. From the ADMET we have selected the best 18 compounds and these compounds were employed for the docking studies by using AutoDock Vina software. Finally, three compounds were chosen and MDS study of 100 ns with the apo-NMDAR and NMDAR-DK1 was performed on them. Then we have carried out various structural parameter analyses, in particular RMSF, Rg, SASA, hydrogen bonds analysis, and PCA. From all these analyses, we have proposed that ZINC4258884 and ZINC15675934 can act as a lead compound against the NMDAR to treat AD. The worldwide scientists can test these compounds employing both in-vivo and in-vitro methods.

Figure 1
An illustrative diagram to represent the NMDAR activation and its binding site.