Alzheimer’s disease, a devastating neurodegenerative disorder, has emerged as a global healthcare challenge. At present, dementia affects more than 55 million people globally, with over 60% residing in low- and middle-income nations. Each year, nearly 10 million new cases emerge. Alzheimer’s disease constitutes the most prevalent form of dementia, accounting for 60–70% of diagnoses [1]. The complex pathophysiology of Alzheimer’s disease is marked by the accumulation of abnormal protein aggregates in the brain, including beta-amyloid plaques [2] and tau fibrillary tangles [3, 4] which ultimately results in cognitive decline and memory impairment [5]. Alzheimer’s disease exhibits two neuropathies: Positive Lesions and Negative Lesions [6]. Positive Lesions, caused due to aggregation, include pathological changes like the formation of amyloid plaques and neurofibril tangles. Negative Lesions, caused by deterioration, are characterized by large atrophy (shrinking of the brain) due to neural and synaptic losses [7, 8].
Plaques are the accumulation of beta-amyloid (Aβ) proteins [9] around the neurons. These Aβ deposits are formed by the catalytic action of β-secretase (BACE-1) [10] and γ-secretase enzymes [11], causing the cleaving of transmembrane amyloid precursor proteins (APP) [12]. APPs on cleavage form multiple insoluble Aβ monomers or peptides, which are sticky and adhere to form plaques between the neurons [13]. This leads to an interruption of neuron-to-neuron signaling, amyloid angiopathy [14] and inflammation in the brain producing neurotoxic [15] effects and hence, cognitive impairments [16]. Neuro-fibrillary tangles (NFTs) are the intracellular aggregation of abnormal, hyper-phosphorylated tau proteins [17]. Phosphate kinase gets activated due to a cascade of reactions initiated by beta-amyloid plaques and transfers the phosphate group to tau proteins causing their misfolding. Structurally altered tau proteins can no longer support microtubules and start forming tangles in the neuro-cytoplasm, axons and dendrites. This alteration in neurons sends the signal for programmed cell death and nerve cells start undergoing apoptosis leading to the shrinkage of the brain [18]. The atrophy of brain cells causes memory loss, language impairment, loss of motor skills, and disorientation and ultimately makes the person bed-ridden leading to death [16, 19]. Despite significant research endeavors, the quest for effective therapies for Alzheimer’s disease continues to be challenging.
Inspiration from nature has prompted exploration into potential therapies for Alzheimer’s disease [20], leading to a focus on the neuroprotective medicinal plants Withania somnifera, Bacopa monnieri, Centella asiatica, and Crocus sativus [21]. These plants have a rich history of traditional use, with various parts of the plants believed to possess medicinal properties. Withania somnifera (also known as Ashwagandha), a member of Solanaceae family is known to possess neuroprotective [22, 23], hepatoprotective [24], anti-inflammatory, antioxidant [25] and anti-depressant [26] properties. Bacopa monnieri (family: Scrophulariaceae) also called Brahmi, has been used for years as a brain-tonic in Ayurveda [27], known for treating various neurological disorders such as improving memory, thinking skills, insomnia, seizures, and anxiety [28]. It has also been studied for its therapeutic potential to treat or prevent neurodegenerative diseases like Alzheimer’s [29] and Parkinson’s disease [30]. Centella asiatica from Apiaceae family, is another medicinally important [31] plant with antioxidant [32], antimicrobial [33], anticancer [34], neuroceutical and cogniceutical properties [35, 36]. It is known to reduce oxidative stress, Aβ levels, and apoptosis, promotes dendritic growth and mitochondrial health, improve mood and memory [21]. Crocus sativus (Saffron), a member of the Iridaceae family has long been known to possess pharmacological properties [37] like antimicrobial, antioxidant [38], cytotoxic, neuroprotective [39] and antidepressant [40] properties.
Computational biology approaches serve as invaluable resources in screening potential therapeutic agents for a particular disease target. Interactions between bioactive compounds found in the neuroprotective plants and key enzymes involved in Alzheimer’s disease pathology can be simulated and predicted virtually. Molecular docking [41] is an in-silico technique used to predict the preferred orientation and conformation of a small molecule (ligand) when bound to a target protein (receptor). By performing molecular docking studies, the interaction of natural bioactive compounds with key enzymes involved in Alzheimer’s disease pathology, such as beta-amyloid protein and beta-secretase (BACE1) could be understood. Molecular dynamics (MD) [42] is a computational method used to simulate the movement and behavior of atoms and molecules over time. These simulations provide insights into the dynamic nature of protein-ligand interactions at the atomic level. By integrating ADMET prediction and drug-likeness assessment into the drug discovery process, we can prioritize compounds with favorable pharmacokinetic and pharmacodynamic properties while minimizing the risk of toxicity and optimizing the likelihood of clinical success. This investigation aims to bridge the gap between traditional herbal knowledge and modern computational techniques by exploring the potential of natural bioactive compounds as novel candidates for Alzheimer’s disease therapy.

Retrieval of ligands

The plants selected for the present study Withania somnifera, Bacopa monnieri, Centella asiatica, and Crocus sativus were chosen based on their traditional use as neuroprotective agents and documents efficacy in treating neurological disorders. An extensive literature review guided the selection of 39 phytochemicals, focusing on phytocompounds with diverse chemical classes, known pharmacological activities, and drug-like properties. The drug Donepezil [43] was chosen as a standard control. The 3D-structures of selected ligands and drugs were retrieved from PubChem [44] database (https://pubchem.ncbi.nlm.nih.gov/) in sdf file format. The compound structures were energy minimized using Chem-3D before docking studies. (Table 1)

Retrieval of target proteins and preparation

This investigation targets two major proteins related to Alzheimer’s disease, such as acetylcholinesterase (AChE) (PDB ID: 1B41) and β-secretase (BACE-1) (PDB ID: 1TQF). The 3D structure of selected proteins were downloaded from RCSB Protein Data Bank (https://www.rcsb.org/pdb/) [45] in pdb format (Figure 1). The enzyme AChE is responsible for hydrolysis of acetylcholine to acetic acid and choline. Acetylcholine is a crucial neurotransmitter in the central nervous system playing a significant role in brain signaling and cognitive functioning. A study in the 1970s revealed that patients of Alzheimer’s disease show a deficiency of Acetylcholine in the brain [46]. This leads to a cholinergic deficit, an inability to transmit neurological impulses across cholinergic synapses, and a consequent cognitive decline [47]. AChE has additionally been implicated in the formation of plaque. AChE inhibitors serve as the most effective symptomatic treatment for AD by slowing down the biochemical breakdown of acetylcholine, leading to prolonged neurotransmission. It may also modify the deposition of β-amyloid [48].

BACE-1 enzyme is accountable for the aberrant degradation of APP [10], which play a role in neuronal growth and repair. The BACE-1 cleaves APP, resulting in the production of insoluble peptides that aggregate to form beta-amyloid plaques. The present investigation aims to examine the potential inhibitory effects of diverse bioactive compounds on the BACE-1 enzyme. By specifically targeting BACE-1, the abnormal degradation of amyloid precursor protein (APP) and subsequent formation of plaque can be mitigated.
Both the proteins were prepared before docking by the removal of water molecules, heteroatoms, and extra ligands. Then checking and repairing missing atoms, the addition of polar hydrogens, the addition of Kollman charges, and spreading the charge deficit equally. Grid maps were also created for the proteins covering their active sites. For 1B41, the grid size was set to 74×74×74 (XYZ) points, and the coordinates for the grid center were set to x= 126.361, y= 102.556, and z= -121.778. For 1TQF, the grid size was set at 50×50×50 (XYZ) and the grid-centre coordinates were x= 28.250, y= 45.805, and z= 2.86.

Molecular docking

Molecular docking between ligands and target proteins is carried out to analyze the fit of ligands in protein’s active site [41]. The binding energy of each ligand is calculated for different poses/conformations of the ligand with the protein. The ligands and conformations showing the least binding energy indicate a favorable and stable fit with the protein. Autodock vina software [49] was used to perform the molecular docking between selected ligands and target proteins. Discovery studio (DS) was used to assess the optimal configuration for each ligand-protein interaction according to binding energies. The ligands with the highest binding energy in the interaction profile were chosen for further examination through molecular dynamics simulations.

Molecular dynamics

The present investigation employed molecular dynamics (MD) simulations using the academic version of the Desmond program (version 2.0) to assess the structural stability of receptor-ligand complexes) [50]. The system in this program was constructed using the TIP3P water model [51] with a cubic periodic box containing simple point charge (SPC) (10Å×10Å×10Å) and optimised potentials for liquid simulations (OPLS) all-atom force field 2005 [52]. The system was subsequently neutralized by introducing the requisite amount of sodium ions. A receptor-ligand complex was provided for the energy minimization step and pre-equilibration for various confined steps. The molecular dynamic simulations were conducted using the OPLS 2005 force field parameters and periodic boundary conditions in the NPT ensemble system [53, 54]. The system was maintained at a constant temperature of 300 K with a relaxation time of 1 ps. The volume was kept constant, and the smooth particle mesh ewald (PME) method was employed with a tolerance limit of 10-9 mm. A cut-off distance of 9.0Å was used. The investigation of protein structures was conducted at a production time of 100 ns, with evaluations performed every 1 ns. To determine the stability, a model structure from the molecular dynamics (MD) simulation during the production phase was chosen. In addition, the root mean square deviation (RMSD) and root mean square fluctuation (RMSF) to examine any structural changes that occurred during the dynamic interaction between the receptor and ligand complexes were also investigated [55].

Binding free energy calculations

ADME and toxicity analysis

Various characteristics were examined to evaluate the drug-like properties of the selected phytochemicals. Lipinski’s rule, executed via molinspiration, was employed to analyze characteristics like the number of hydrogen acceptors (less than 10), the number of hydrogen donors (less than 5), molecular weight (more than 500 Daltons), and partition coefficient log P (not less than 5) [59]. The selected phytocompounds were subjected to their toxicity attributes. The 3D structures of the best-docked phytocompounds and Donepezil were converted to SMILES format and evaluated for toxicity using admetSAR [60] and PROTOX-III webserver [61].

Molecular docking

The results with the least (most negative) binding energy values are considered strong interactions and reflect a stable and favorable fit. The binding energies (kcal/mol) of the selected compounds with the two target proteins (1B41 and 1TQF) are tabulated in Table 2. Retinoic acid showed a strong binding energy of -9.2 kcal/mol with 1B41 protein, while Somniferine showed the highest binding energy value of -8.8 kcal/mol with 1TQF protein. Donepezil showed binding energy values of -8.3 kcal/mol and -8.4 kcal/mol with 1B41 and 1TQF, respectively (Table 2).

Discovery Studio Visualizer was used to analyze the interactions of retinoic acid with 1B41, and Somniferine with 1TQF (Figure 2). Retinoic acid formed a hydrogen bond with HIS (A):287, and two pi-sigma bonds with TYR (A):341 and TRP (A):286 residues of 1B41 protein. Other non-covalent interactions were formed with TYR (A):124 and VAL (A):294 residues of 1B41. Somniferine was found to form two hydrogen bonds with LYS (A):321 and ARG (A):307 residues of 1TQF protein. Other non-covalent interactions were formed with ILE (A):110 and ARG (A):307 residues of 1TQF. Interactions of Donepezil with both the target proteins are also shown in Table 3.

Molecular dynamics

Molecular dynamics simulation is used to study the dynamic behavior of protein-ligand complexes over time. It provides insights into the motions, interactions, and structural changes of atoms and molecules at the atomic or molecular level within a protein-ligand complex. Retinoic acid and Somniferine demonstrated strong binding energy when interacting with 1B41 and 1TQF, respectively. Subsequently, the protein-ligand complexes of retinoic acid-1B41 and somniferine-1TQF were screened through molecular dynamics simulations for 100 ns.

RMSD analysis of protein-ligand complexes

The RMSD reflects the average deviation in the positions of selected atoms of a simulated complex from a reference structure over a trajectory. The deviation for each frame in the trajectory is calculated. Figure 3 displays the RMSD evolution of the protein (left y-axis) and the ligand (right y-axis). The starting point is defined by the docked arrangement of the ligand and the protein in the complex. The deviation from this reference position during the MD simulation is then gauged by aligning all protein-ligand frames obtained throughout the trajectories. The term "lig fit prot" (line colored in magenta) implies the RMSD of the ligand from the backbone of the protein. For retinoic acid-1B41 complex, the protein backbone hovers within the range of 2Å-3.4 Å and fluctuation initially up to 30 ns. The complex was then stable between 30 to 80 ns in between 2-2.8 Å, then showed fluctuation between 80-90 ns between 2.0-3.2 Å (Figure 3A). The RMSD plot of somniferine-1TQF complex (Figure 3B) shows a more fluctuating curve. The protein backbone hovers in the range from 1.2 Å to 2.8 Å from 0-100 ns. The lig-fit-prot curve remained well below the protein backbone curve in both protein-ligand complexes, indicating that the orientation of the ligand remains the same (Figure 3).

RMSF analysis of protein-ligand complexes

The RMSF measures the localized structural changes along the protein chain during molecular dynamics simulations. Peaks on the RMSF plot indicate regions of the protein that experience the maximum fluctuations throughout the simulation. Secondary structure elements, such as α-helices and β-strands, highlighted with red and blue backgrounds, respectively, are characterized by lower fluctuations, signifying their inherent rigidity. In contrast, loop regions, often unstructured, tend to exhibit higher fluctuations, emphasizing their dynamic nature. The green-colored vertical bars mark the protein residues that encounter the ligand. These interactions play a crucial role in understanding the binding dynamics and stability of the complex.
To complement the RMSF analysis, the plot also correlates the fluctuation data with experimental B-factors. While the RMSF and B-factor definitions differ, a parallel between simulation results and crystallographic data is anticipated. This correlation enhances the reliability of the MD simulation outcomes, offering a valuable perspective on the agreement between computational predictions and experimental observations. The trajectory of 1B41-retinoic acid complex shows a plot with minimal residual fluctuations in the range of 0.5-1.2 Å and almost parallel RMSF and B-Factor curves (Figure 4A). The plots of Somniferine with 1TQF initially displayed few differences in the RMSF and B-Factor values (Figure 4B) till residue index 100-120, and after that they show overlapping RMSF and B-Factor curves. 1TQF-somniferine displayed residual fluctuation in the range of 1.0-2.0 Å.

Ligand properties

In the MD simulation, the behaviour of the two best performing phytocompounds, Retinoic acid and Somniferine was assessed over a 100 ns timeframe, with key properties evaluated to understand their stability and interaction with the surrounding environment (depicted in Figure 5). For the retinoic acid, The Root Mean Square Deviation (RMSD) fluctuated between 0.5 Å and 1.5 Å while for Somniferine, the range remained below 0.6 Å, indicating that the ligands maintained a stable structure with minor deviations from their initial conformations (time t = 0). The Radius of Gyration (rGyr), which measures the ‘extendedness’ of the ligand, fluctuated slightly between 4.8 Å to 5.1 Å for Retinoic acid and between 5.12 Å and 5.28 Å for Somniferine, showing that the ligands retained compact shape throughout the simulation. Notably, no intramolecular hydrogen bonds (intraHB) were detected in both simulations, suggesting that the ligands did not rely on internal hydrogen bonds to maintain their structural stability. The molecular surface area (MolSA), which represents the molecular surface calculated with a 1.4 Å probe radius (comparable to the van der Waals surface), ranged from 332.5 Ų to 340 Ų in Retinoic acid and between 474 Ų - 486 Ų in Somniferine, indicating minimal changes in the ligand’s surface. The Solvent Accessible Surface Area (SASA), which measures the extent to which the ligand is exposed to solvent (water), fluctuated between 50 Ų and 200 Ų for Retinoic acid and between 120 Ų and 300 Ų in Somniferine. Higher SASA values could imply that the ligand is more exposed and possibly at risk of being pulled away from the active site of the protein by the solvent, while lower values indicate stronger interaction with the protein. The polar surface area (PSA), which accounts for the solvent-accessible areas contributed by polar atoms (oxygen and nitrogen), ranged from 87 Ų to 93 Ų in retinoic acid and 150 Ų to 168 Ų in somniferine.

Protein-ligand contacts histogram

The interactions between the protein and the ligand are studied throughout the simulation and depicted through the histogram. Four types of interactions are summarized in the plots (Figure 6): hydrogen bonds (green), hydrophobic (grey), ionic (pink), and water bridges (blue). Retinoic acid complexed with 1B41 showed one ionic bond with Asp (A): 74; water bridges with Asn (A): 283, His (A): 287, Leu (A): 289, Pro (A): 290, and Arg (A): 296; polar contacts with Leu (A): 76, Tyr (A): 77, Trp (A): 86, Tyr (A): 124, Phe (A): 295, Phe (A): 297, Tyr (A): 337, Phe (A): 338 and Tyr (A): 341 (Figure 6A). Somniferine formed the polar contacts interactions with Phe (A):148, Met (A):149, Arg(A):156, Phe (A):161, Leu (A):166, Pro 178, Met (A):194, Phe (A): 213 residues of 1TQF protein (Figure 6B). Figure 6C and 6D depict the interactions of crucial amino acids of selected target proteins with retinoic acid and Somniferine. Several residues establish several explicit interactions with the ligand, as seen by a more diffuse shade of orange, as indicated by the scale on one side of the figure. These plots are extremely important because they show that retinoid acid and Somniferine interact with the amino acids of target proteins throughout the simulation and do not dissociate from their interacting site. On the other hand, small changes in the ligand’s RMSD and RMSF values, as shown in Figure 2 and 3, respectively, imply that the ligands may be reorienting themselves during the simulation.

Preservation of intermolecular contacts in molecular dynamics simulations

In the 1B41-Retinoic acid complex, the ligand interacts with the TYR341 residue of the protein. The interaction is primarily hydrophobic, indicated by the green label, reflecting non-polar side chain interactions between the tyrosine residue and the hydrophobic components of the ligand. The ligand’s structure, which includes extensive conjugation represented by alternating single and double bonds in the E-configuration along with bulky cyclic structures, enhances this interaction by maximizing hydrophobic contact. In the 1TQF-Somniferine complex, several key interactions were observed during the molecular dynamics simulation. The GLU (A):141 residue forms a hydrogen bond with the ligand, preserved for 83% of the simulation time, and contributes to stabilizing the interaction through electrostatic attraction due to its negative charge. The PHE (A):148 residue engages in a Pi-Pi stacking interaction with the ligand’s aromatic rings, with a 41% occurrence during the simulation, enhancing hydrophobic interactions. Additionally, the LYS (A):145 residue forms a Pi-cation interaction with the ligand, maintained for 99% of the simulation time, reflecting a highly stable contact between the positively charged lysine residue and the ligand’s aromatic ring. These interactions, including hydrogen bonds, Pi-Pi stacking, and Pi-cation contacts, play a significant role in maintaining the structural stability of the ligand-protein complex. Figure 7 presents the 2D interaction maps of the best-docked compounds, illustrating the maintenance of contacts throughout the simulation trajectory.

Binding-free energy calculations

Post-simulation analysis of both protein-ligand complexes was carried out by capturing snapshots of the trajectory profiles generated during MD simulations, as shown in Table 4. It was observed that both protein-ligand complexes exhibited negative ∆G_binding energies, suggesting the stability of these complexes during the MD simulation. The results showed that the Van der Waals interactions (∆G_vdW) were -44.13 ± 5.22 kcal/mol for the retinoic acid-1B41 complex and -56.97 ± 2.76 kcal/mol for Somniferine-1TQF complexes, indicating that both the ligands remain close to the interacting amino acids of the target proteins. The Coulomb energy displayed a negative value in all complexes, indicating consistently low potential energy for ligands when bound to their respective target proteins. The finding indicates that protein-ligand complexes generally exhibit improved stability. Table 4 provides an extensive overview of the contributions made by different components, such as hydrogen bonding and covalent bonding, to the total energy, along with the associated total energy values.

Drug-likeness and toxicity analysis

Table 5 displays the results of Lipinski’s rule of five for the best-docked phytocompounds and selected drugs that were chosen. Both the selected phytocompounds and selected drugs were found to follow Lipinski’s rule of five. The results from toxicity prediction showed that retinoic acid and Somniferine were non-neurotoxic, non-carcinogenic, and non-cytotoxic in nature, whereas, the drug, Donepezil was positive for all these parameters. Retinoic acid was positive for hepatotoxicity, while Somniferine was immunogenic in nature. The predicted LD₅₀ (mg/kg) was found to be 1100 (Class-4) for both retinoic acid and Somniferine, whereas the predicted LD₅₀ for Donepezil was found to be 550 mg/kg (Class-4), indicating the non-toxicity of both the phytocompounds.

Computational drug discovery methods have emerged as efficient pillars in the lifecycle of drug-development [62]. Taking the ‘experimental set-up’ to ‘virtual environment’ has allowed the screening of even billions of molecules for their drug-likeness properties against many disease targets [63]. In-silico analyses have significantly reduced the time and cost of identifying suitable molecules as drugs, expected to reduce clinical trial failures at later stages. The chemical space of compounds that can be used as drugs is huge, with a large part of them coming from traditional herbal knowledge. Natural compounds have shown promising therapeutic potential and reduced toxicity values [64] and provide cost-effective measures in terms of synthesis and collection of raw material. Building upon traditional knowledge of such natural herbs, this investigation aimed to find potential therapeutic compounds to treat Alzheimer’s disease.
Several plants have long been used traditionally as brain tonics and neuroprotective agents [65]. Our study aimed to explore the binding interactions between natural compounds from such neuroprotective herbs including Withania somnifera, Bacopa monnieri, Centella asiatica, and Crocus sativus, and enzymes associated with Alzheimer’s disease. Among the 39 compounds screened, retinoic acid (PubChem ID: 444795) and somniferine (PubChem ID: 14106343) emerged as promising candidates for Alzheimer’s disease therapy as they exhibited strong binding affinities (-9.2 kcal/mol and -8.8 kcal/mol, respectively) and stable interactions with AChE (PDB ID: 1B41) and BACE-1 (PDB ID: 1TQF), respectively. Prior research emphasizes the significance of AChE and BACE-1 inhibition in mitigating cognitive symptoms linked to Alzheimer’s [66-68]. Retinoic acid was established as a potent AChE inhibitor in this study, a key target for Alzheimer’s treatment [69, 70]. BACE-1 inhibitors are another class of compounds studied against Alzheimer’s [11, 47]. Somniferine has been estimated to be a promising BACE-1 inhibitor in this investigation.
Phytocompounds like Withaferine [71] and Withanolides [72] from Withania somnifera have been studied by researchers for their potential against AD. Our study established Somniferine, an alkaloid found in Withania as another novel candidate for Alzheimer’s disease treatment. Another compound, retinoic acid, found in many plants including Crocus sativus was established as a lead compound to manage Alzheimer’s disease via AChE inhibition. Retinoic acid has also been previously studied as a therapeutic option for AD through cholinergic restoration [73]. Additionally, computational toxicity prediction plays a crucial role in drug discovery by facilitating early compound screening and reducing clinical trial failures [74]. The favorable ADMET properties of retinoic acid and Somniferine underscore their potential as viable Alzheimer’s treatment candidates. However, in silico analyses are preliminary, and have their limitations. Understanding how drugs interact with the body and predicting their real-world outcomes can only be confirmed through experimental validation, given the dynamic and complex nature of biological systems. It’s imperative to conduct comprehensive in vitro and in vivo studies to thoroughly assess the efficacy and toxicity profile of lead compounds.

In conclusion, this research aims to connect traditional herbal wisdom with modern scientific methods, investigating how bioactive compounds from traditional herbs could be used therapeutically for Alzheimer’s disease. The study focused on 39 key compounds found in neuroprotective herbs and employed a cheminformatics approach, combining computational biology techniques such as molecular docking and molecular dynamics simulations with pharmacological and ADMET analyses. The selected compounds were thoroughly analyzed for drug-likeness, pharmacokinetic properties, and toxicity, revealing valuable insights into their potential as drug candidates for Alzheimer’s disease therapy. Among 39 selected phytocompounds, retinoic acid and Somniferine exhibited particularly promising binding interactions, demonstrating adherence to the Lipinski rule, and showing low toxicity. Both the phytocompounds outperformed the standard AD drug, Donepezil, in terms of binding affinity and toxicity properties. The molecular dynamics simulations provided a dynamic perspective on the stability and fluctuations of the protein-ligand complexes over time. Retinoic acid and somniferine displayed stable interactions with respective target enzymes (1B41 and 1TQF). However, further experimental validation, including in vitro and in vivo studies, will be crucial to substantiate the observed interactions and therapeutic potential. This research serves as a promising point for unraveling the hidden pharmacological treasures within nature. By merging traditional knowledge with cutting-edge computational techniques, the study contributes to the ongoing quest for innovative treatments for Alzheimer’s disease, a global healthcare challenge.