Introduction

Alzheimer’s disease (AD) stands as a formidable and increasingly prevalent neurodegenerative disorder, posing substantial challenges to healthcare systems worldwide. Characterized by progressive cognitive decline, memory impairment, and behavioral changes, AD significantly impacts the quality of life for affected individuals and their families1. Approximately 10% of individuals aged 65 years experience Alzheimer’s, while the prevalence significantly escalates to affect nearly one-third of individuals aged 85 years and older2. As the most prevalent cause of dementia among the elderly, Alzheimer’s disease poses a substantial burden on healthcare systems worldwide3.

Globally, an estimated 50 million people were living with dementia in 2018, a number projected to triple by 2050, reaching 152 million individuals4. Alzheimer’s disease accounts for the majority of these cases, constituting a significant portion of the global dementia burden. With an aging global population, the incidence and prevalence of Alzheimer’s are expected to rise dramatically in the coming decades5.

In 2016, Alzheimer’s and other dementias were the fifth leading cause of death worldwide, responsible for 2.4 million deaths6. The impact of Alzheimer’s extends beyond mortality, profoundly affecting the quality of life for individuals and their families. The economic burden associated with Alzheimer’s is substantial, with global costs exceeding US$1 trillion in 2018, a figure expected to escalate as the prevalence of the disease increases7. These statistics underscore the urgent need for continued research, public awareness, and healthcare strategies to address the growing challenges posed by Alzheimer’s disease on a global scale.

As the global population ages, the incidence of AD is expected to rise, intensifying the urgency to explore diverse avenues for both understanding its pathogenesis and identifying potential therapeutic interventions8.

The etiology of Alzheimer’s disease is complex, involving a multitude of genetic, environmental, and lifestyle factors. The accumulation of beta-amyloid plaques and tau tangles in the brain is a hallmark pathological feature, leading to synaptic dysfunction and neuronal loss9. Additionally, dysregulation of neurotransmitter systems, particularly the glutamatergic system, has been implicated in the cognitive decline observed in AD10. The N-methyl-D-aspartate (NMDA) receptor Fig. 1, a key player in glutamate-mediated neurotransmission, has garnered significant attention in the context of AD pathogenesis11. In recent research, attention has expanded beyond traditional targets like the NMDA receptor to explore novel avenues for Alzheimer’s disease intervention. Emerging studies are investigating targets such as monoamine oxidase (MAO), traditionally associated with conditions like depression and Parkinson’s disease12,13. This evolving understanding underscores the intricate nature of Alzheimer’s etiology and the ongoing efforts to identify diverse therapeutic targets for this complex neurodegenerative disorder14,15.

Fig. 1
figure 1

The three-dimensional structural representation of the NMDA receptor providing a detailed examination from various perspectives, each designated with specific colors for chain identification. In the front view, the A-chain stands out distinctly in a striking Red color, offering a clear visualization of its positioning. Transitioning to the top view, the B-chain takes center stage with its prominent Blue color, facilitating a comprehensive overhead perspective. The bottom view highlights the C-chain, uniquely designated in Green, aiding in a focused examination of its spatial arrangement. Completing the structural overview, the Yellow color specifically designates the D-chain, ensuring clarity in identifying its position within the overall structure.

Full size image

One intriguing avenue for potential therapeutic intervention in Alzheimer’s disease lies in the realm of traditional Indian Ayurvedic herbs16. Among these, Withania somnifera, commonly known as Ashwagandha, has emerged as a subject of interest due to its reputed neuroprotective properties. Ayurveda, the ancient Indian system of medicine, has long recognized the diverse therapeutic benefits of Ashwagandha, considering it an adaptogen that enhances resilience to stress and promotes overall well-being17.

The neuroprotective potential of Ashwagandha is attributed to its rich composition of bioactive compounds, including alkaloids, steroidal lactones (withanolides), and flavonoids18. These constituents are believed to exert antioxidant, anti-inflammatory, and anti-apoptotic effects, which could be particularly relevant in the context of neurodegenerative disorders like Alzheimer’s disease19,20,21. Moreover, Ashwagandha has been reported to modulate various neurotransmitter systems, including the cholinergic and GABAergic systems, suggesting a multifaceted impact on brain function22.

Advancements in computational techniques have revolutionized drug discovery, offering an efficient and cost-effective approach to identify potential therapeutic agents23. In the context of Alzheimer’s disease (AD), where the need for effective treatments is urgent, computational methods such as virtual screening, molecular docking, and molecular dynamics simulations have played a pivotal role in accelerating the drug discovery process24,25,26.

In the realm of Alzheimer’s disease (AD) treatment, various therapeutic strategies have been developed with the objective of enhancing cholinergic brain function27. One prominent approach involves the introduction of drugs designed to stimulate cholinergic receptors, thereby facilitating the development of acetylcholine, a neurotransmitter crucial for cognitive function. Simultaneously, these drugs work to impede the degradation of acetylcholine by inhibiting the action of cholinesterase, an enzyme responsible for its breakdown. Additionally, these medications exhibit anti-inflammatory properties by modulating enzymes involved in inflammatory processes28.

The current focus of drug interventions for AD centers around specific agents known as acetylcholinesterase inhibitors. This class of drugs includes donepezil, galantamine, and rivastigmine29. Their primary mechanism of action involves blocking the activity of acetylcholinesterase, thereby preserving higher levels of acetylcholine in the brain. This approach aims to mitigate the cognitive decline associated with AD by enhancing cholinergic neurotransmission30.

Another significant player in AD pharmacotherapy is memantine, which operates through a distinct mechanism. Memantine acts as a non-competitive inhibitor of N-methyl-D-aspartate (NMDA) receptors31. NMDA receptors play a crucial role in synaptic transmission and are involved in processes such as learning and memory32. By modulating the activity of NMDA receptors, memantine helps regulate glutamate, another neurotransmitter implicated in AD pathology33. The diversity in the mechanisms of action between acetylcholinesterase inhibitors and memantine highlights the multifaceted nature of Alzheimer’s disease and the need for a comprehensive treatment approach. While acetylcholinesterase inhibitors address the cholinergic deficits, memantine provides an alternative avenue by targeting glutamatergic dysfunction, offering a more holistic strategy for managing the complex pathophysiology of AD34.

In the context of Alzheimer’s disease treatment, the pharmacological landscape extends beyond drugs directly addressing cholinergic deficits. A noteworthy control compound in this ___domain is Ifenprodil, a drug known for its role as a non-competitive antagonist of the N-methyl-D-aspartate (NMDA) receptor35. Unlike memantine, which is a non-competitive inhibitor, Ifenprodil exerts its influence by blocking a specific subunit of the NMDA receptor, namely the NR2B subunit36.

Ifenprodil’s mechanism involves selectively modulating the NMDA receptor function, particularly targeting the NR2B subunit, which has implications for synaptic plasticity and memory processes37. While Ifenprodil is not a primary drug used in AD treatment, its inclusion as a control compound is instrumental in experimental settings. Researchers often leverage Ifenprodil to compare and validate the efficacy of novel compounds targeting the NMDA receptor, contributing to a deeper understanding of the intricate processes involved in glutamatergic dysfunction associated with AD38.

Keeping in view of the complex landscape of Alzheimer’s disease treatment, our present research endeavors are intricately directed towards the meticulous identification of the most promising active compound derived from Withania somnifera. This pursuit is rooted in the explicit goal of discerning a compound with the inherent potential to intricately modulate the N-methyl-D-aspartate (NMDA) receptor function. The selected compound is expected to pave ways for further investigations as a prospective candidate for Alzheimer’s disease therapeutics. This strategic focus aligns with the imperative to unravel novel avenues in AD treatment, with the ultimate aim of contributing substantially to the ongoing discourse in the field.

Methods

Protein structure retrieval, preparation and visualization

The preparation for molecular docking commenced with the acquisition of the three-dimensional structure of the N-methyl-D-aspartate (NMDA) receptor in PDB format. Given the crystallographic availability of our protein of interest in the RCSB Protein Data Bank39, we retrieved the NMDA receptor structure with the accession code 4PE5. The retrieved structure underwent meticulous scrutiny and rectification of any structural irregularities using the BioVia Discovery Studio Visualizer version 17.240. This step ensured precision by addressing missing atoms or residues.

Hydrogen atoms were then added to the protein utilizing MGLTools and PyMOL, employing PyMOL’s “AddH” tool. This critical step is pivotal for accurate docking calculations, ensuring precise hydrogen placement. Subsequently, Gasteiger charges were assigned to the protein atoms within MGLTools using the “Assign Charges” tool, laying the groundwork for the computation of electrostatic interactions during the docking process. The final step involved the transformation of the protein structure into the Autodock-compatible PDBQT format. This crucial conversion, executed with MGLTools’ “Prepare Protein” tool, incorporated essential details such as atom types and charges, prerequisites for the ensuing successful docking simulations, as elaborated elsewhere41,42.

Concurrently, the three-dimensional structures of the ligands were retrieved from the PubChem database43 in SDF format. These structures were subject to energy minimization using Biovia Discovery Studio and subsequently converted to PDB format. Following this, the “Prepare Ligand” tool within MGLTools facilitated the conversion of the ligand structures into the Autodock-compatible PDBQT format. This format encapsulates vital information, including atom types and charges, ensuring the ligands are suitably equipped for precise Autodock docking calculations. This comprehensive methodology guarantees the accuracy and reliability of subsequent docking simulations, contributing to the robustness of our investigative approach.

Withania somnifera active compounds database and ligand preparation

For the identification of phytochemicals present in Withania somnifera, a curated database, namely the Indian Medicinal Plants, Photochemistry, and Therapeutics (IMPPAT) (https://cb.imsc.res.in/imppat/home)44, was systematically employed. The 2D structures of both the identified phytochemicals and the GluN2B NMDAR inhibitor, Ifenprodil, were retrieved from the National Center for Biotechnology Information (NCBI) PubChem in .sdf format.

Subsequently, all phytochemicals, as detailed in Supplementary Material 1, underwent meticulous preparation using the open babel software. This preparation involved the conversion of the 2D structures into the AutoDock Vina-supported .pdbqt format, followed by a crucial step of universal force field45 energy minimization. This rigorous preparation ensures the compatibility of the ligands for subsequent molecular docking studies, facilitating a robust and accurate exploration of their interactions with the GluN2B N-methyl-D-aspartate receptors (NMDARs).

Virtual screening and molecular docking

Virtual screening of the 80 compounds targeting the N-methyl-D-aspartate (NMDA) receptor protein was meticulously conducted using PyRx AutoDock Vina (Version 1.1.2)46. This initial step aimed to filter compounds exhibiting favorable binding energy values from the pool. Subsequent molecular docking and interaction analyses were performed using AutoDock version 1.5.6 default protocol, extensively detailed elsewhere47.

AutoDock version 1.5.6 provided a comprehensive platform for identifying binding energy, pIC50 values, and binding poses of the compounds. The process involved two main steps: grid generation and docking. Grid generation established a spatial framework directing ligand binding, with dimensions of X-axis = 126, Y-axis = 126, and Z-axis = 126, along with a spacing of 0.375. Docking employed Lamarckian genetic algorithms, encompassing 2,500,000 evolutions and generating 27,000 docking poses. Before starting the docking with the compounds, we validated the docking protocol by redocking the co-crystallized compound and got an approximate RMSD of 0.9 angstroms. The top 10 ligands displaying the most favorable interactions with the protein were subsequently identified.

AutoDock provides binding affinity scores, usually in the form of Gibbs free energy of binding (ΔG), which reflects how energetically favorable the binding interaction is. In AutoDock, the inhibition constant (Ki) is not directly calculated by the docking tool itself but can be inferred from the docking results. AutoDock primarily focuses on predicting the binding affinity of a ligand to a target protein. Here’s a outline of how you the inhibition constant from AutoDock results was derived:

To convert the binding affinity (ΔG) to an inhibition constant (Ki), first the binding free energy was used to calculate the dissociation constant (Kd). The Kd is related to ΔG through the following equation:

$$\Delta G\, = - RT\,ln\,{K_d}$$

where ΔG is the binding free energy, RRR is the gas constant (1.987 cal/(mol·K)), and T is the temperature in Kelvin (typically 298 K for physiological conditions). Rearranging this equation to solve for Kd gives:

$$Kd\, = \,{e^{ - \Delta G/\left( {RT} \right)}}$$

In many cases, Ki can be approximated as the Kd for the ligand if you’re assuming a simple competitive inhibition scenario. To ensure the accuracy and reliability of the molecular docking simulations, both proteins and ligands underwent energy minimization using BioVia Discovery Studio Visualizer version 17.248. This software not only facilitated the energy minimization process but also enabled a detailed analysis of molecular-level interactions within the docked protein-ligand complexes.

The energy-minimized proteins and ligands were seamlessly integrated into the Autodock software through the mgltools interface. Before proceeding with the blind docking procedure, the structures were converted into pdbqt files, adhering to standard protocols elucidated elsewhere49. This meticulous approach ensured the robustness of the blind docking procedure, offering insights into potential ligand binding sites and interaction patterns at the atomic level.

Molecular dynamic simulations

Molecular dynamics simulations were systematically conducted to validate the binding interactions of the most potent binder identified from virtual screening, a candidate validated through molecular docking. The atomic-level analysis of these interactions was performed using the default protocol of the Desmond v3.6 module, extensively detailed elsewhere50,51.

In essence, the simulations employed the OPLS 2005 force field52 parameters, coupled with TIP3P water models53 under neutral pH conditions. Periodic boundary conditions were employed to establish the specific size and shape of the water box, buffered at 10 Å of the simulation box. The equilibration process involved treating Van der Waals and short-range electrostatic interactions with a cutoff of 9 Å, while long-range electrostatic interactions were computed using the Particle Mesh Ewald method54. A Reversible Reference System Propagation Algorithm (RESPA) integrator55 was utilized with a time step of 2 fs, and long-range electrostatics were computed every 6 fs.

The simulation system, which included approximately 604,798 atoms for the NMDA-17alpha-hydroxywithanolide-D complex along with a POPC membrane at a temperature of 300 K, underwent equilibration in the NPT ensemble at 300 K and 1 bar. The Nose-Hoover chain relaxation thermostat method and Martyna-Tobias-Klein relaxation Barostat method, with isotropic coupling style at 1 ps and 2 ps timescales, respectively, ensured consistent temperature and pressure conditions. Rigorous checks were undertaken to confirm uniform temperature, pressure, and volume conditions throughout the simulated timescale.

As a component of the simulation quality analysis, it was ascertained that the average total potential energy of the simulated systems remained approximately − 1,820,000 Kcal/mol for the NMDA-17alpha-hydroxywithanolide-D complex. This meticulous simulation protocol and subsequent analysis guarantee the reliability and accuracy of the acquired molecular dynamics data, enabling a comprehensive understanding of the dynamic behavior and stability of the studied complex.

Results & discussion

Virtual screening

In this study, virtual screening was meticulously conducted for all 80 ligands sourced from Withania somnifera, utilizing a curated database and literature from the IMPAAT database. Employing PyRx AutoDock Vina, binding affinities of these ligands towards the NMDA receptor protein were evaluated. The grid box was intelligently designed to envelop the entire protein, and subsequent docking scores or binding energies were subjected to thorough analysis.

The outcomes of the virtual screening revealed compelling insights. Across the spectrum of 80 ligands, all compounds demonstrated a noteworthy potential to bind to diverse domains of the NMDA receptor, showcasing binding energies within the range of -4.1 to -11.9 Kcal/mol. For comparative context, the FDA-approved Ifenprodil exhibited a binding energy of -7.8 Kcal/mol.

Further scrutiny of the results unveiled distinct binding patterns among the investigated compounds. Seven compounds exhibited binding affinity to the A-chain, 37 to the B-chain, 7 to the C-chain, and 29 to the D-chain of the NMDA receptor. Intriguingly, visual inspection of the bound compounds elucidated a predominant localization at the core binding site of the Amino Terminal Domain (ATD). Additionally, some compounds occupied the interface between ATD and the Ligand Binding Domain (LBD), while only a select few were positioned at the Transmembrane Domain (TMD) (Fig. 2). Figure 2 snapshots provide a glimpse into the molecular landscape, highlighting the specific points of interaction between the screened compound and different regions of the NMDA receptor. The visual representation offers a valuable insight into the spatial distribution and potential binding sites within the receptor’s complex architecture. This snapshot serves as a concise yet informative illustration of the compound’s binding dynamics across the diverse domains of the NMDA receptor chains.

A standout candidate, 17alpha-hydroxywithanolide D, demonstrated exceptional promise with a striking binding energy of -11.9 Kcal/mol. Given its notable performance, this compound was earmarked for further in-depth investigations through molecular docking and dynamic simulation studies. This comprehensive virtual screening underscores the potential of Withania somnifera-derived compounds as valuable candidates for targeted intervention in the modulation of NMDA receptor function, presenting avenues for future research in Alzheimer’s disease therapeutics.

Fig. 2
figure 2

The visual snapshot captures a moment in the dynamic interaction between the NMDA receptor chains and the virtually screened compound. Each receptor chain—A, B, C, and D—is depicted distinctly, showcasing the binding of the screened compound at various domains.

Full size image

Structure-activity relationship (SAR) study: establishing the proposed compound as a partial agonist of NMDA receptor

We have investigated the structure-activity relationship (SAR) of various NMDA receptor modulators, with a focus on determining the role of the potential compound found through virtual screening studies. Through comparative docking studies, we analyzed inactive, moderate, and active NMDA receptor ligands, including a co-crystallized ligand (QEL) compared to our proposed compound 17 alpha-hydroxy withanolide D, to assess binding affinities and their potential mechanisms of action (Table 1).

Inactive compounds

  • AP5 (D-APV), a selective NMDA receptor antagonist, but considered less effective in vivo due to its inability to cross the blood-brain barrier efficiently56 - displayed poor binding affinity with docking scores between − 1.6 kcal/mol and − 3.27 kcal/mol and Ki values ranging from 4.04 mM to 67.48 mM. AP5’s low binding affinity and high Ki values indicate minimal interaction with the glycine site, confirming its role as an inactive antagonist. The structure of AP5, characterized by its bulky, hydrophilic groups, may limit its ability to engage effectively with the receptor’s binding pocket.

  • L-689,560, A competitive glycine site antagonist at the NMDA receptor. Although it can inhibit NMDA receptor function, it is often categorized as having low efficacy compared to more potent modulators57 - exhibited better binding than AP5 but remained within the inactive category, with docking scores between − 5.85 kcal/mol and − 6.98 kcal/mol and Ki values ranging from 7.61 µM to 51.34 µM. Its moderately improved affinity can be attributed to a more favorable alignment within the binding site. However, the steric hindrance imposed by its molecular structure likely reduces its efficacy, contributing to its inactivity at the receptor.

Moderate compounds

  • Ketamine, A non-competitive antagonist of the NMDA receptor, often used as an anesthetic and also for its rapid antidepressant effects. While potent, its effects can vary based on dosage, hence categorized as moderate in certain contexts58 - showed moderate binding affinity with docking scores between − 7.93 kcal/mol and − 8.00 kcal/mol and Ki values ranging from 1.36 µM to 1.55 µM. Ketamine’s cyclic structure and lipophilic groups enable better interaction within the NMDA receptor’s glycine site. Its moderate binding suggests partial inhibition of receptor activity, aligning with its pharmacological role as a moderate NMDA antagonist.

  • D-cycloserine, A partial agonist at the NMDA receptor glycine site. It has been studied for enhancing cognitive function and as an adjunct to exposure therapy in anxiety disorders59 - demonstrated docking scores between − 6.65 kcal/mol and − 6.74 kcal/mol and Ki values between 11.48 µM and 13.30 µM. Its moderate binding affinity is indicative of partial activation, consistent with its structure that allows for both receptor engagement and partial agonistic action. The presence of smaller functional groups facilitates limited receptor activation, which is a hallmark of partial agonism.

Active compounds

  • MK-801 (Dizocilpine), A non-competitive NMDA receptor antagonist widely used in research as a potent and selective inhibitor. It is used to block NMDA receptor activity and induce neuroprotection in models of excitotoxicity60 - displayed high binding affinity with docking scores between − 8.59 kcal/mol and − 8.60 kcal/mol and Ki values in the nanomolar range (500.46 nM to 507.53 nM). The rigid bicyclic structure of MK-801 allows for optimal interaction with the receptor, leading to full antagonism. This strong binding correlates with its established role as a full blocker of NMDA receptor activity.

  • Memantine, A moderate-affinity NMDA receptor antagonist used clinically to treat Alzheimer’s disease. It blocks excessive calcium influx without affecting normal synaptic activity61 - showed even stronger binding with docking scores between − 8.67 kcal/mol and − 8.96 kcal/mol and Ki values ranging from 269.61 nM to 445.00 nM. Memantine’s symmetrical structure and well-defined hydrophobic interactions within the glycine site enable high-affinity binding, supporting its clinical use as a neuroprotective agent. The compact, rigid structure facilitates tight receptor blockade, indicating full antagonistic activity.

Co-crystallized Ligand (QEL)

The co-crystallized ligand QEL of NMDA receptor with PDB ID 4PE562 displayed docking scores between − 7.78 kcal/mol and − 9.42 kcal/mol, with Ki values ranging from 125.24 nM to 1.97 µM. These values suggest strong binding at the glycine site, serving as a reference for the structural alignment and binding potential of other compounds. QEL’s structure, containing both hydrophobic and hydrophilic regions, demonstrates balanced interactions that enhance receptor engagement, validating the docking methodology.

Proposed compound 17alpha-hydroxywithanolide D

The proposed compound 17 alpha-hydroxy withanolide D showed docking scores ranging from − 7.47 kcal/mol to -10.03 kcal/mol and Ki values between 44.24 nM and 3.33 µM, indicating significant binding affinity to the NMDA receptor glycine site. Notably, the variation in docking scores suggests a partial agonistic behavior, similar to known partial agonists like D-cycloserine.

Structurally, the proposed compound possesses functional groups that enable flexible interactions with the glycine site, allowing for partial receptor activation. The relatively small size and balanced hydrophobic/hydrophilic nature of the compound suggest it can engage the receptor without inducing full activation, consistent with partial agonism. The higher docking scores in certain poses (-10.03 kcal/mol) suggest that the compound can achieve stable binding. However, the moderate scores in other poses (-7.47 kcal/mol) reflect partial receptor engagement rather than full agonism.

The variability in Ki values (ranging from 44.24 nM to 3.33 µM) further supports this hypothesis, as partial agonists often display a range of binding affinities based on receptor conformation and the degree of receptor activation.

Structural insights and mechanism

The proposed compound exhibits structural features conducive to partial agonism. The flexible moieties within its molecular structure allow it to partially occupy the binding site, resulting in moderate receptor activation. The interaction profile is similar to D-cycloserine, where the compound induces submaximal receptor activity. Unlike full antagonists like MK-801 and Memantine, the proposed compound’s structural flexibility likely prevents it from fully blocking the receptor, allowing for controlled modulation of receptor activity.

Based on the SAR analysis, the proposed compound demonstrates key structural characteristics and binding profiles that suggest it may act as a partial agonist at the NMDA receptor glycine site. The variation in docking scores and Ki values, combined with the structural flexibility of the compound, align with the pharmacological behavior of partial agonists. These findings warrant further experimental validation through functional assays to confirm its potential as a therapeutic modulator of NMDA receptor activity.

Table 1 Structure-activity relationship (SAR) of various NMDA receptor modulators: through comparative docking studies, we analyzed inactive, moderate, and active NMDA receptor ligands, including a co-crystallized ligand (QEL) to our proposed compound 17 alpha-hydroxy withanolide D:
Full size table

Molecular docking

Molecular docking, recognized for its robustness in predicting binding poses across a spectrum of entities including ligands, peptides, proteins, and nanomaterials, stands as a cornerstone in structural biology42,43. However, its limitation in capturing dynamic binding effects has spurred the integration of molecular dynamics, offering a more comprehensive approach that incorporates both static and dynamic considerations44.

In this study, we seamlessly combined molecular docking with molecular dynamics simulations to elucidate the molecular-level interactions within the NMDA-17alpha-hydroxywithanolide D complex. Molecular docking, the initial phase, demonstrated the impressive binding capability of the investigated 17alpha-hydroxywithanolide D compound. The interaction energies reached a substantial magnitude of -10.03 Kcal/mol, underscoring its robust binding affinity. Furthermore, the predicted IC50 value of 44.24 nano molar attested to its potential as a potent inhibitor.

A distinctive feature unveiled by the docking studies was a hydrogen bond formation between 17alpha-hydroxywithanolide D and the Lys466 residue of the D-chain, a pivotal interaction elucidated in Fig. 3. This interaction not only substantiates the binding capability of the compound but also provides valuable insights into the molecular-level events that govern its interaction with the NMDA receptor.

The subsequent integration of molecular dynamics simulations ensures a comprehensive exploration of the dynamic behavior and stability of the NMDA-17alpha-hydroxywithanolide D complex. By amalgamating these two methodologies, our study aims to unveil a nuanced understanding of the intricate interplay between the compound and the receptor, paving the way for potential therapeutic applications and further advancements in NMDA receptor target specific drug discovery.

Fig. 3
figure 3

Molecular level interactions between NMDA-17alpha-hydroxywithanolide D complex. Left panel shows the interactions in 2D, whereas the right panel shows the interactions in three dimensional representation.

Full size image

MD simulations of NMDA-17alpha-hydroxywithanolide D complex

In our pursuit to unravel the dynamics of the NMDA-17alpha-hydroxywithanolide D complex, we conducted extensive molecular dynamics simulations spanning 100 nanoseconds (100,000 picoseconds). The resulting simulation data, pivotal for understanding the impact of 17alpha-hydroxywithanolide D binding on the dynamics of the NMDA receptor protein, was subjected to thorough analysis.

Key parameters, including the Root Mean Square Deviation (RMSD) of the protein backbone, energy fluctuations, and the percentage of Secondary Structural Elements (SSE), were systematically evaluated across the simulated timescale. The RMSD of the protein backbone exhibited fluctuations within the range of 4.0 to 8.0 Å, with an average of 7.5 Å, as illustrated in Fig. 4.

Detailed analysis of the RMSD graph unveiled noteworthy insights. The NMDA-17alpha-hydroxywithanolide D complex demonstrated substantial stability in its overall structure, particularly beyond the initial 30 nanoseconds of the simulated timescale. Notably, a progressive expansion towards higher RMSD values was observed as the ligand RMSD shifted from 3.0 Å to 6.0 Å. This intriguing observation suggests a modulatory impact exerted by the 17alpha-hydroxywithanolide D compound on the structural dynamics of the NMDA receptor, indicating a dynamic interplay between the compound and the receptor over the course of the simulation.

Fig. 4
figure 4

Root mean square deviation (RMSD) of protein and ligand backbone of NMDA-17alpha-hydroxywithanolide D complex. X-axis represents the timeline of the simulation in terms of nanoseconds, whereas the y-axis represents the calculated RMSD in angstrom units.

Full size image

In order to validate and reinforce the overall stability of the NMDA-17alpha-hydroxywithanolide D complex, a comprehensive analysis of the energy parameters governing the simulated system was undertaken. The potential energy of the simulated system emerged as a critical focus, providing crucial insights into the system’s behavior. The analysis revealed that the potential energy of the NMDA-17alpha-hydroxywithanolide D complex approximated − 1,820,000 Kcal/mol. This substantial negative potential energy underscores the robust stability of the complex. The pronounced negative binding energies, exemplified by this considerable value, suggest a formidable inhibition of the NMDA-17alpha-hydroxywithanolide D complex. This observation is further substantiated by the enormous negative binding energies, as depicted in Fig. 5.

In essence, the formidable negative binding energies affirm the strong inhibitory potential of the NMDA-17alpha-hydroxywithanolide D complex, providing a quantitative measure of the complex’s stability and reinforcing its potential as a potent inhibitor within the context of the studied molecular dynamics simulations.

Fig. 5
figure 5

Calculated potential energies of the simulated systems of NMDA-17alpha-hydroxywithanolide D complex. X-axis represents the timeline of the simulation in terms of nanoseconds, whereas the y-axis represents the calculated potential energy in kcal/mol units.

Full size image

In our pursuit to comprehensively understand the influence of the 17alpha-hydroxywithanolide D compound binding on the structural characteristics of the NMDA protein, a meticulous investigation into the Secondary Structural Elements (SSE) was undertaken throughout the simulated timescale. The comprehensive analysis of the Secondary Structural Elements (SSE) within the NMDA-17alpha-hydroxywithanolide D compound complex offers a detailed portrayal of the distribution and behavior of key protein structural elements, particularly alpha-helices (depicted in orange) and beta-strands (highlighted in blue). Monitored throughout the entirety of the simulation, this analysis provides insights into the SSE distribution by residue index across the protein structure.

As illustrated in Fig. 6, the NMDA-17alpha-hydroxywithanolide D compound complex consistently maintained an average of 39.90% total SSE. Within this composition, alpha helices accounted for 30.14%, beta helices constituted 9.77%, while the remainder was comprised of loops. This visualization enhances our understanding of how the NMDA-17alpha-hydroxywithanolide D complex influences the distribution and maintenance of alpha-helices and beta-strands over the simulated timescale. Such detailed structural insights contribute to unraveling the complex interplay between the compound and the NMDA receptor, shedding light on the mechanisms through which 17alpha-hydroxywithanolide D may modulate the structural dynamics of the protein.

This analysis offers valuable insights into the preservation of structural integrity within the NMDA-17alpha-hydroxywithanolide D complex. The maintenance of a substantial percentage of alpha helices and beta helices, crucial elements of protein structure, underscores the structural stability conferred by the binding interaction with the 17alpha-hydroxywithanolide D compound.

The preservation of these secondary structural elements throughout the simulation provides further evidence of the compound’s influence on maintaining the structural robustness of the NMDA receptor. This information contributes to a comprehensive understanding of the complex’s behavior, emphasizing its potential as a modulator of the structural dynamics of the NMDA protein.

Fig. 6
figure 6

Analysis of Secondary Structural Elements (SSE) of NMDA-17alpha-hydroxywithanolide D compound complex. Protein secondary structure elements (SSE) like alpha-helices (orange colored) and beta-strands (blue colored) are monitored throughout the simulation. As delineated in the plot above, the X-axis corresponds to the amino acid residues index, spanning the entire protein structure. On the other hand, the y-axis represents the calculated SSE percentage concerning the simulated timescale. The distinctive coloring of alpha-helices and beta-strands aids in visually discerning their respective contributions to the overall secondary structure.

Full size image

Molecular interactions observed during molecular dynamic simulations

A thorough examination of the interactions within the NMDA-17alpha-hydroxywithanolide D complex revealed a total of 23 contacts, showcasing the complexity of molecular associations. Among these, hydrogen bonds emerged as crucial contributors, involving six specific residues. Concurrently, hydrophobic interactions, known for enhancing structural stability, engaged six residues, further fortifying the complex. Additionally, water binding interactions, a noteworthy mediator in biomolecular dynamics, encompassed 21 residues. Notably, direct hydrogen bonding interactions with D-chain residues ASP165, ARG431, THR433, LYS466, and TYR476 emerged as pivotal, underscoring their critical role in orchestrating the formation and stability of the complex. This comprehensive interaction analysis provides a nuanced understanding of the diverse molecular contacts within the NMDA-17alpha-hydroxywithanolide D complex, shedding light on the intricacies of its binding mechanism and emphasizing the significance of specific residues in driving these interactions (Fig. 7). The c-panel from Fig. 6 illustrates the total number of specific contacts established by the protein with the ligand throughout the trajectory. Simultaneously, the bottom panel delineates the residues involved in interactions with the ligand in each trajectory frame. Notably, darker shades of orange signify residues making multiple specific contacts, providing a nuanced depiction of the dynamic and multifaceted nature of the interactions. This timeline representation serves as a valuable tool for discerning the temporal evolution and diversity of contacts between the NMDA receptor and 17alpha-hydroxywithanolide D, contributing to a comprehensive understanding of the complex’s behavior.

Fig. 7
figure 7

Comprehensive Analysis of Molecular Interactions in NMDA-17alpha-hydroxywithanolide D Complex: (a) Molecular Dynamics Simulations Interactions: During the extensive molecular dynamics simulations of the NMDA-17alpha-hydroxywithanolide D complex, intricate molecular-level interactions were observed. (b) Interaction Bar Chart Panel: The bar chart panel above provides a graphical representation of ligand-interacting amino acid residues on the x-axis and their respective interaction percentages concerning the simulated timescale on the y-axis. (c) Timeline Representation of Interactions and Contacts: The timeline representation encapsulates a comprehensive summary of various contacts, including H-bonds, hydrophobic interactions, ionic interactions, and water bridges.

Full size image

Our investigation has identified specific residues within the D-chain of the NMDA receptor—ASP165, ARG431, THR433, LYS466, and TYR476—as pivotal players in the binding interactions with Withania somnifera-derived 17alpha-hydroxywithanolide D.

Understanding the importance of targeting these residues provides a nuanced perspective on potential therapeutic strategies for Alzheimer’s disease (AD). The identified residues within the NMDA receptor, namely ASP165, ARG431, THR433, LYS466, and TYR476, play distinctive roles in ligand-receptor interactions. ASP165 forms critical hydrogen bonds, contributing significantly to stabilizing the complex. ARG431, as a positively charged residue, is integral to electrostatic interactions, influencing receptor modulation and signaling pathways. THR433, known for hydrogen bonding, plays a crucial role in maintaining structural stability, impacting the receptor’s conformation and downstream signaling. LYS466 catalyzes essential electrostatic interactions vital for ligand binding, with mutations in lysine residues potentially altering receptor kinetics. Finally, TYR476, with its aromatic nature, engages in hydrophobic interactions crucial for stabilizing ligand binding. Understanding the specific contributions of these residues provides valuable insights into the molecular mechanisms governing NMDA receptor functionality and opens avenues for targeted interventions in modulating receptor behavior.

The outcomes of our investigation highlight the pivotal role of direct hydrogen bonding interactions, coupled with water bridging interactions, in facilitating the binding of the 17alpha-hydroxywithanolide D compound to the NMDA receptor protein. Notably, the specific hydrogen bonding was observed with the hydroxyl group of the 17alpha-hydroxywithanolide D compound, underscoring the precision and specificity of this interaction.

Moreover, the interactions within the NMDA-17alpha-hydroxywithanolide D complex were found to be influenced by various factors, including molecular weight, hydrophobicity, and the amino acid composition and sequence of the compound. The distinctive nature of these interactions emphasizes the nuanced response of NMDA receptor proteins to compounds, providing insights into the complex interplay between the ligand and the receptor.

Specifically, the binding mechanism involved the hydrophobic regions of the protein receptor interacting with the non-polar aromatic ring of the 17alpha-hydroxywithanolide D compound. This hydrophobic interaction contributes to the overall stability of the complex, showcasing the importance of structural complementarity in the binding process.

Additionally, electrostatic interactions were identified as contributors to the binding mechanism, facilitated by the hydroxyl groups of the 17alpha-hydroxywithanolide D compound and charged groups on the proteins. These electrostatic interactions further underscore the multifaceted nature of the molecular interactions governing the binding of 17alpha-hydroxywithanolide D to the NMDA receptor protein.

Discussion

Prior investigations have established through 3D classification that the GluN1b-GluN2B heterodimers exhibit three principal structures denoted as ‘Non-active1,’ ‘Non-active2,’ and ‘Active.’63 These structural classifications are characterized by specific conformations, and their distinctions are crucial in understanding the modulations within the receptor complex. The key divergences among these 3D classes revolve around three major aspects:

  1. 1.

    Subunit Arrangement within GluN1b-GluN2B ATD Heterodimers: The configuration of subunits within the ATD heterodimers is delineated by the distance between GluN1b Lys178 and GluN2B Asn184 in the R2 lobes.

  2. 2.

    Relative Orientation of GluN1b-GluN2B LBD Heterodimers: The positioning of LBD heterodimers is defined by the distance between GluN1b Arg510 in loop1 (L10) and GluN2B Leu425 in loop2 (L2) in the D1 lobes, known as the L10-L2 distance.

  3. 3.

    Relative Tension of LBD-M30 Linker of GluN2B Subunits: The tension within the LBD-M30 linker of GluN2B subunits, a critical factor influencing the channel gate, is quantified by the distance between the two GluN2B Gln662 residues from the respective subunits.

To gain insights into the plausible mode of action of the compounds, three representative snapshots from the simulation, captured at the first frame, 50th frame, and the last frame, were meticulously analyzed. This approach allowed for a comprehensive examination of the compounds’ potential impact on the intricate dynamics of the GluN1b-GluN2B receptor complex. The assessment, based on the specified criteria, provides valuable perspectives on the compounds’ interactions within the receptor environment, contributing to a nuanced understanding of their modulatory effects.

Subunit Arrangement within GluN1b-GluN2B ATD Heterodimers: In the earlier study, the identification of distinct conformational states within the GluN1b-GluN2B receptor complex was elucidated by monitoring the dynamic interplay between GluN1b Lys178 and GluN2B Asn184 in the R2 lobes of ATDs. This dynamic interaction is governed by ‘rolling’ motions intrinsic to GluN2B ATDs. The delineation of these states was based on the measured distances: 17.4 angstroms characterized as non-active1, 17.8 angstroms as non-active2, 12.6 angstroms as active, and the shortest at 12.0 angstroms as active-SS. Notably, a shorter distance implies heightened activity, while a longer or relaxed distance indicates a non-active state.

Our current investigation has shed light on the impact of our proposed compound on the NMDA receptor. In the presence of our compound, the receptor displayed a distance of 8.3 angstroms in the first frame, 12.29 angstroms at the 50ns mark, and 10.36 angstroms in the last frame of the simulation. These findings distinctly demonstrate the capability of the proposed compound to modulate and inhibit the activity of the NMDA receptor, providing crucial insights into its potential therapeutic role in regulating receptor dynamics (Fig. 8).

Fig. 8
figure 8

Visualization of NMDA subunit re-arrangement within GluN1b-GluN2B ATD Heterodimers during the simulated timescale from molecular dynamic simulations trajectory analysis. The configuration of subunits within the ATD heterodimers is delineated by the distance between GluN1b Lys178 and GluN2B Asn184 in the R2 lobes at (a) first frame, (b) 50ns time scale frame and (c) last frame of the simulation.

Full size image

Relative orientation of GluN1b-GluN2B LBD heterodimers: In the previous investigation, the dynamic positioning of LBD heterodimers within the GluN1b-GluN2B receptor complex was defined through the L10-L2 distance, measured as the distance between GluN1b Arg510 in loop1 (L10) and GluN2B Leu425 in loop2 (L2) in the D1 lobes. The categorization of states was based on this measured distance, with 18.8 angstroms representing non-active1, 19.0 angstroms as non-active2, 13.0 angstroms as active, and the shortest at 9.1 angstroms as active-SS. A shorter L10-L2 distance indicated heightened activity, while a longer or relaxed distance signified a non-active state.

Our current investigation has unveiled the impact of our proposed compound on the NMDA receptor’s L10-L2 distance. In the presence of our compound, the receptor exhibited a distance of 37.32 angstroms in the first frame, 45.97 angstroms at the 50ns mark, and 41.54 angstroms in the last frame of the simulation. These findings convincingly establish the compound’s efficacy in inhibiting the activity of the NMDA receptor, providing valuable insights into its potential therapeutic role in regulating the dynamic conformational states of the receptor (Fig. 9).

Fig. 9
figure 9

Visualization of NMDA relative orientation of GluN1b-GluN2B LBD heterodimers during the simulated timescale from molecular dynamic simulations trajectory analysis: The positioning of LBD heterodimers is defined by the distance between GluN1b Arg510 in loop1 (L10) and GluN2B Leu425 in loop2 (L2) in the D1 lobes, known as the L10-L2 distance at (a) first frame, (b) 50ns time scale frame and (c) last frame of the simulation.

Full size image

Relative tension of LBD-M30 Linker of GluN2B subunits: In a previous investigation, the study highlighted the pivotal role of the LBD-M30 linker tension within the GluN2B subunits, a determining factor influencing the channel gate of the NMDA receptor. This tension was quantified by measuring the distance between the two GluN2B Gln662 residues from their respective subunits. The classification of states was established based on this measured distance, with 50.3 angstroms indicative of non-active1, 51.0 angstroms for non-active2, 52.8 angstroms denoting an active state, and a stretched 61.0 angstroms representing active-SS. A longer distance between the GluN2B Gln662 residues signified increased activity, whereas a shorter distance maintained tension in the LBD-M30 linker, characterizing a non-active state.

In our current investigation, the influence of our proposed compound on the NMDA receptor’s LBD-M30 linker tension has been elucidated. In the presence of our compound, the receptor exhibited a distance of 49.60 angstroms in the first frame, 37.95 angstroms at the 50ns mark, and 35.19 angstroms in the last frame of the simulation. These results provide compelling evidence for the proposed compound’s efficacy in inhibiting the activity of the NMDA receptor, specifically through its impact on maintaining the tension within the LBD-M30 linker (Fig. 10).

Fig. 10
figure 10

Visualization of relative tension of LBD-M30 linker of GluN2B subunits during the simulated timescale from molecular dynamic simulations trajectory analysis: The tension within the LBD-M30 linker of GluN2B subunits, a critical factor influencing the channel gate, is quantified by the distance between the two GluN2B Gln662 residues from the respective subunits at (a) first frame, (b) 50ns time scale frame and (c) last frame of the simulation.

Full size image

Conclusion

In this exhaustive exploration, we have intricately probed the molecular nuances dictating the interaction between Withania somnifera-derived 17alpha-hydroxywithanolide D and the NMDA receptor, holding profound implications for Alzheimer’s disease (AD) therapeutics. Employing advanced computational methodologies such as virtual screening, molecular docking, and dynamic simulations, our study has unveiled promising compounds, particularly highlighting the potential of 17alpha-hydroxywithanolide D in the context of AD treatment. Virtual screening of 80 Withania somnifera compounds uncovered candidates with substantial binding affinity across diverse domains of the NMDA receptor. Noteworthy among them, 17alpha-hydroxywithanolide D exhibited a remarkable binding pose, garnering attention for its potential therapeutic application.

Molecular docking simulations provided detailed insights into the binding orientation, while dynamic simulations shed light on the stability, energy parameters, and secondary structural elements of the NMDA-17alpha-hydroxywithanolide D complex over a simulated timescale. Our study underscored the critical role of specific residues, such as ASP165, ARG431, THR433, LYS466, and TYR476, through direct hydrogen bonding interactions in orchestrating the binding of 17alpha-hydroxywithanolide D to the NMDA receptor. The color-coded 3D structural representation and snapshots not only visually delineated the interaction dynamics but also accentuated the importance of these residues in the intricate molecular interplay among receptor chains.

Our investigation further delved into the dynamic conformational states within the GluN1b-GluN2B receptor complex, elucidating the intricate interplay between GluN1b Lys178 and GluN2B Asn184 in the R2 lobes of ATDs. The observed modulation of this interaction, categorized into non-active1, non-active2, active, and active-SS states, was influenced by the proposed compound. The NMDA receptor, in the presence of our compound, demonstrated distinct distances in its first frame (8.3 angstroms), 50ns mark (12.29 angstroms), and last frame (10.36 angstroms), showcasing the compound’s ability to modulate and inhibit receptor activity.

Furthermore, our study examined the relative orientation of GluN1b-GluN2B LBD heterodimers through the L10-L2 distance. The compound’s impact on this distance, demonstrated by values of 37.32 angstroms, 45.97 angstroms, and 41.54 angstroms in the first frame, 50ns mark, and last frame, respectively, underscored its efficacy in inhibiting NMDA receptor activity. Lastly, the investigation into the tension within the LBD-M30 linker of GluN2B subunits provided crucial insights. The proposed compound exhibited the capability to influence the NMDA receptor’s LBD-M30 linker tension, with distances of 49.60 angstroms, 37.95 angstroms, and 35.19 angstroms in the first frame, 50ns mark, and last frame, respectively.

In conclusion, our multifaceted approach not only validates the therapeutic potential of Withania somnifera-derived compounds for AD, especially 17-alpha-hydroxy withanolide-D as potential partial agonist unraveling a nuanced understanding of the molecular intricacies. The inclusion of residue-specific information not only enhances our comprehension of the binding mechanism but also provides a roadmap for targeted drug development. This study propels us toward innovative therapeutic strategies, enriching the arsenal against Alzheimer’s disease with a more precise and informed approach through further functional assays to confirm its potential as a therapeutic modulator of NMDA receptor activity.