Multitargeted Virtual Screening and Molecular Simulation of Natural Product-like Compounds against GSK3β, NMDA-Receptor, and BACE-1 for the Management of Alzheimer’s Disease

The complexity of Alzheimer’s disease (AD) and several side effects of currently available medication inclined us to search for a novel natural cure by targeting multiple key regulatory proteins. We initially virtually screened the natural product-like compounds against GSK3β, NMDA receptor, and BACE-1 and thereafter validated the best hit through molecular dynamics simulation (MDS). The results demonstrated that out of 2029 compounds, only 51 compounds exhibited better binding interactions than native ligands, with all three protein targets (NMDA, GSK3β, and BACE) considered multitarget inhibitors. Among them, F1094-0201 is the most potent inhibitor against multiple targets with binding energy −11.7, −10.6, and −12 kcal/mol, respectively. ADME-T analysis results showed that F1094-0201 was found to be suitable for CNS drug-likeness in addition to their other drug-likeness properties. The MDS results of RMSD, RMSF, Rg, SASA, SSE and residue interactions indicated the formation of a strong and stable association in the complex of ligands (F1094-0201) and proteins. These findings confirm the F1094-0201’s ability to remain inside target proteins’ binding pockets while forming a stable complex of protein-ligand. The free energies (MM/GBSA) of BACE-F1094-0201, GSK3β-F1094-0201, and NMDA-F1094-0201 complex formation were −73.78 ± 4.31 kcal mol−1, −72.77 ± 3.43 kcal mol−1, and −52.51 ± 2.85 kcal mol−1, respectively. Amongst the target proteins, F1094-0201 have a more stable association with BACE, followed by NMDA and GSK3β. These attributes of F1094-0201 indicate it as a possible option for the management of pathophysiological pathways associated with AD.


Virtual Screening of Natural Product-like Compounds
Computational screening of a vast number of small organic molecules for their inhibitory potential against the target proteins is widely recognized and accepted, which can significantly reduce the time, cost and efforts of wet lab high-throughput screening [47][48][49][50]. In this study, we applied virtual screening through molecular docking and validation of best hits through molecular dynamics simulation (MDS) to identify the novel inhibitor for multiple targets (NMDA, GSK3β, and BACE) to manage neurological disorders. Primarily, the molecular docking protocol was validated by redocking the native ligands and found that it bounds to almost similar residues. RMSD value (≤2) between docked and native ligands was in an acceptable limit. Then, the natural product-like compounds (NPLC) library and reference inhibitors/substrate (Native ligand of GSK3β: Adenosine-5 -Diphosphate, Native ligand of NMDA receptor: 5,7-Dichlorokynurenic acid, and Native ligand of BACE: non-peptidic inhibitor) were individually docked on the active site of target proteins NMDA, GSK3β, and BACE, respectively. The results demonstrated that out of 2029 compounds total of 288 compounds have binding energy scores between −10 and −12.3 Kcal/mol for NMDA, 135 compounds exhibited binding energy scores between −9 and −11.2 Kcal/mol for GSK3β, and 213 compounds showed binding energy scores between −10 and −12 Kcal/mol for BACE. Of these best active 288, 135, and 213 compounds against NMDA, GSK3β, and BACE, respectively, only 51 compounds exhibited binding interactions with all three targets and were considered multitarget inhibitors. The top 10 best hits of compounds against all three targets were chosen (Table 1) for further analysis.

Drug-Likeness, Pharmacokinetics, and Physicochemical Properties
This study used the SwissADME tool to analyze the physicochemical, pharmacokinetics, and drug-likeness properties of the 10 best hits (multiple targets inhibitors; Table 2). The significance of adhering to these parameters was well established. It has been said that most medications failed throughout the drug development process because they were ineffective at adhering to these parameters [51,52]. The physiological properties of the majority of orally active drugs, such as molecular weight (MW), hydrogen bond donors (HBD), hydrogen bond acceptors (HBA), and XlogP, were found to be in a specified range (MW: 160-500 g/mol, HBD: 5, HBA: 10, and XlogP: −1-6). Poor oral bioavailability is represented by chemical structures with more than 10 rotatable bonds [53]. In addition, the molar refractivity (MR) range was regarded as between 40 and 130 for more excellent intestine absorption [54]. Therefore, according to Lipinski et al. (1997), the compounds must adhere to the acceptable range for at least three physiochemical attributes out of five to be considered drug-like [52]. Our results (Table 2) illustrated that out of the 10 best hits, only five (F1094-0201, F1217-0041, F1094-0205, F1094-0196, and F1094-0198) compounds could penetrate the blood-brain barrier (BBB) which is the most important property for any central nervous system medication [55]. Furthermore, all the compounds except F3161-0307 showed high GI absorption, whereas five (F0870-0001, F1217-0041, F1094-0205, F1094-0196, and F1094-0206) compounds were found to be Pgp-substrate. Those compounds exhibited as good Pgp-substrate are unsuitable for CNS drug candidates [56]. These results illustrate that all 10 best hits except one compound (F3161-0307) followed the Lipinski rule of five, having an acceptable range of MW, HBD, HBA, XlogP, and MR. Therefore, after the analysis of these results, only two compounds (F1094-0201 and F1094-0198) were found to be suitable for CNS drug-likeness in which BOILED-Egg image for F1094-0201 to predict gastrointestinal absorption (HIA) and brain penetration (BBB) has been shown due to its higher inhibitory potential against multiple targets ( Figure 1). Our results also illustrated that all the best hits follow Veber's rule as they have fewer than 10 rotatable bonds, and the TPSA is under the acceptable limit for drug-likeness, which is less than 140 Å 2 which showed that our best hits have good oral bioavailability and better penetration possibility [53].

Molecular Docking Study
The virtual screening of two thousand twenty-nine natural product-like compounds and blood-brain barrier penetration capacity of the top 10 best hits illustrate that the F1094-0201 is the most effective inhibitor against multiple targets (NMDA, GSK3β, and BACE) with binding energy −11.7, −10.6, and −12 kcal/mol, respectively. Further analysis of the ligand-protein complex revealed the interacting amino acid residues and the type of interactions involved in stabilizing the complexes. The interactions of F1094-0201 with the active site residues of NMDA, GSK3β, and BACE are shown in Figures 2-4, respectively.

Molecular Docking Study
The virtual screening of two thousand twenty-nine natural product-like com and blood-brain barrier penetration capacity of the top 10 best hits illustrate F1094-0201 is the most effective inhibitor against multiple targets (NMDA, GSK BACE) with binding energy −11.7, −10.6, and −12 kcal/mol, respectively. Further of the ligand-protein complex revealed the interacting amino acid residues and t of interactions involved in stabilizing the complexes. The interactions of F1094-02 the active site residues of NMDA, GSK3β, and BACE are shown in Figures 2-4, tively.
GSK3β protein has two active sites: ATP-binding and substrate-binding sites, where key residues (VAL135 and ASP133) are available at the ATP-binding site, known as the activation loop. Whereas LYS85 and GLU97 also have a prominent role in the catalytic process [64]. The previous report suggests that ARG141 is one of the essential residues for specific ATP/ADP recognition by TPK I/GSK3 beta, and several other residues are of key importance in ATP-binding sites, including ILE62, VAL70, ALA83, LYS85, VAL110, LEU132, GLN185, LEU188, and ASP200 [65]. Our results demonstrate that both the reference ligand (Adenosine-5′-Diphosphate) and F1094-0201 compound commonly interacted with most of the key residues (ILE62, ALA83, LYS85, VAL110, LEU132, GLN185, LEU188, and ASP200). Moreover, F1094-0201 was also found to interact with ASP133. Our results correspond with previously published reports [64,66].
Elevated beta-secretase (BACE) activity may have harmful effects on CNS due to its involvement in the production of Aβ 42 , known for forming AB plaque by protein aggregation. Reduction of Aβ 42 toxicity in SH-SY5Y human neuroblastoma cells was achieved by polyphenols which modulate autophagy against neurodegeneration [40]. Moreover, some neurological disorders are accompanied by α-synuclein (α-syn) misfolding, and reports suggested that fungal extracts known to reduce α-syn toxicity through diverse processes, such as the reduction of protein aggregation, a decrease of the ROS level, and α-syn membrane delocalization supporting their anti-aggregation properties. Hence the reduction of protein aggregates can be a better approach to managing neurological disorders [42]. Therefore, BACE is considered a prime target for delaying amyloid pathology and managing AD [3,57,58]. We used BACE enzyme co-crystallized with hydroxyethyl amine inhibitor at 2.55 Å resolution [59]. The two aspartate (ASP32 and ASP228) residues are involved in the catalytic process of the enzyme [60,61]. Our results showed that a critical catalytic residue ASP32 interacted with the F1094-0201. Our results correspond with the previous report, stating that apart from ASP 32 and ASP228, the BACE inhibitors also interacted with GLY34, TYR71, LYS107, and PHE108 [62].
Previous in-vitro studies suggested that the production of toxic Aβ peptide can also be regulated by GSK3β via affecting presenilin-1 function. Studies conducted in-vitro and in transgenic AD animal models showed that Aβ stimulates GSK3β signaling; further, the brains of AD patients showed a comparable rise in GSK3β activity. Hyperactivation of GSK3β increases the phosphorylation of tau proteins and potentiates the formation of neurofibrillary tangles. However, through the implementation of an NF-kB signalingmediated approach, GSK3β inhibition decreases BACE1-mediated breakage of APP. This finding consequently implies that inhibiting GSK3β lessens the disease associated with Aβ pathology [63].
GSK3β protein has two active sites: ATP-binding and substrate-binding sites, where key residues (VAL135 and ASP133) are available at the ATP-binding site, known as the activation loop. Whereas LYS85 and GLU97 also have a prominent role in the catalytic process [64]. The previous report suggests that ARG141 is one of the essential residues for specific ATP/ADP recognition by TPK I/GSK3 beta, and several other residues are of key importance in ATP-binding sites, including ILE62, VAL70, ALA83, LYS85, VAL110, LEU132, GLN185, LEU188, and ASP200 [65]. Our results demonstrate that both the reference ligand (Adenosine-5 -Diphosphate) and F1094-0201 compound commonly interacted with most of the key residues (ILE62, ALA83, LYS85, VAL110, LEU132, GLN185, LEU188, and ASP200). Moreover, F1094-0201 was also found to interact with ASP133. Our results correspond with previously published reports [64,66].
Neuronal survival is reliant on synaptic NMDA receptor signaling. The overflow of glutamate produced by astrocytes or presynaptic terminals plays a crucial role in antagonizing the synaptic pro-survival signaling pathway and tipping the scales in favor of excitotoxicity and eventual neurodegeneration. Memantine, an FDA-approved NMDA receptor antagonist, has been shown to have positive therapeutic benefits in moderate-tosevere AD patients. It may do this by reducing extra-synaptic NMDA receptor signaling. Therefore it is beneficial to target NMDA receptors for the management of AD [67].
NMDA receptors require both glycine and glutamate for activation, with NR1 and NR2 forming glycine and glutamate sites, respectively. Here we used the high-resolution (1.90 Å) co-crystal structures of the NR1 S1S2 ligand-binding core with the antagonist 5,7-dichloro kynurenic acid (DCKA). The NR1 site has been considered for therapeutic potential [68,69]. Ugale and Bari reported that amino acid residues Arg131, Pro124, and Thr126 are essential for inhibiting the Gly/NMDA receptor [70,71]. These interactions were also noticed by Devid et al. [72] with some additional interactions (GLN13, ASP224, and Trp223 residues). Our results are in correspondence with these previous reports.

Root Mean Square Deviation (RMSD)
A protein-ligand complex's RMSD is often calculated to assess the stability and dynamic properties of the complex. The RMSD is determined by deviating from the proteinligand starting pose's structure as a function of simulation duration [73]. In this experiment, we measured the RMSD of three target proteins (BACE, GSK3, and NMDA) in the presence of the chemical F1094-0201 throughout a simulation time of 100 ns ( Figure 5). The RMSD of BACE increased during the first few seconds and remained consistent throughout the simulation. However The MDS results indicated the establishment of a stable association between proteins and ligands based on the steady-state behavior of RMSD of all target proteins in the presence of F1094-0201. Furthermore, our results suggested that the overall structure of the protein-ligand complex had not undergone any significant conformational changes.

Root Mean Square Fluctuation (RMSF)
The monitoring of RMSF is often used to identify the local conformational changes in a protein's side chains caused by the binding of a ligand [73]. Here, the RMSFs of BACE, GSK3β and NMDA were accessed in the presence or absence of the ligand F1094-0201 ( Figure 6). Fluctuations at proteins' N-and C-terminals are due to their higher flexibilities. The results indicate that the RMSF plots of target proteins in the presence of F1094-0201 were almost overlapping with the RMSF plot of target proteins alone. The average RMSF values of BACE, GSK3β, and NMDA in the presence of F1094-0201 were 0.88 ± 0.03 Å, 0.94 ± 0.06 Å, and 0.79 ± 0.05 Å, respectively, while the average RMSF values of BACE, GSK3β, and NMDA in the presence of F1094-0201 were 0.94 ± 0.06 Å, 1.14 ± 0.08 Å, and 2.82 ± 0.09 Å, respectively. These findings indicate that the protein-ligand combination is stable in nature and that the binding of the F1094-0201 compound did not significantly alter the overall structure of the target proteins.

The Radius of Gyration (Rg) and Solvent-Accessible Surface Area (SASA)
The radius of gyration measures the protein-ligand complex's compactness, and the solvent-accessible surface area determines its exposure to solvent molecules. Each of these characteristics sheds light on the protein-ligand complex's stability throughout the simulation [73]. In this study, the Rg of BACE, GSK3β, and NMDA was determined in the presence of F1094-0201 during 30-100 ns, as shown in Figure 7A  The radius of gyration measures the protein-ligand complex's com solvent-accessible surface area determines its exposure to solvent molec  Monitoring the variations in a protein's secondary structure elements (SSE) in tein-ligand complex is significant for evaluating any structural changes in a prote to ligand binding [74]. Here, we have evaluated the total SSE (α-helix + β-sheet) of GSK3β, and NMDA in the presence of F1094-0201 as a function of simulation time 8). The average SSE of BACE, GSK3β, and NMDA in complex with F1094-0201 wa ± 1.84% (α-helix = 5.57 ± 0.91% and β-sheet = 33.00 ± 2.22%), 38.11 ± 1.54% (α-helix ± 1.08% and β-sheet = 16.02 ± 1.87%), and 41.92 ± 2.56% (α-helix = 22.09 ± 1.37% sheet = 19.84 ± 1.98%), respectively. The results indicate that the total SSE of all the ta proteins in the presence of F1094-0201 did not undergo significant changes, there plying a stable protein conformation in the protein-ligand complex.

Contacts between Protein and Ligand
The total number of contacts formed between protein and ligand was also evaluated over the simulation (Figure 9). The total number of interactions between ligand and protein in BACE-F1094-0201, GSK3β-F1094-0201, and NMDA-F1094-0201 complexes varied in the range of 0-11, 0-11, and 0-12, respectively. During the 30-100 ns simulation, the average contacts between F1094-0201 and BACE, GSK3β, and NMDA were 7, 5, and 6, respectively.

Analysis of Prime/MM-GBSA Free Energy
Free energy calculation by Prime/MM-GBSA is a versatile and accurate method to evaluate the stability of a protein-ligand complex in the presence of solvent molecules [75]. Here, free energy and its constituents were determined using the Prime/MM-GBSA approach, and the results are presented in Table 3. Amongst the target proteins, F1094-0201 formed the most stable complex with BACE, followed by NMDA and GSK3β. The free energies of BACE-F1094-0201, GSK3β-F1094-0201, and NMDA-F1094-0201 complex formation were −73.78 ± 4.31 kcal mol −1 , −72.77 ± 3.43 kcal mol −1 , and −52.51 ± 2.85 kcal mol −1 respectively. In all the protein-ligand complexes, van der Waal forces, Coulombic forces, packing interactions, and lipophilic interactions favored the formation of a stable complex. Conversely, covalent interactions and solvation-free energy opposed the formation of a stable protein-ligand complex.

Analysis of Prime/MM-GBSA Free Energy
Free energy calculation by Prime/MM-GBSA is a versatile and accurate method to evaluate the stability of a protein-ligand complex in the presence of solvent molecules [75]. Here, free energy and its constituents were determined using the Prime/MM-GBSA approach, and the results are presented in Table 3. Amongst the target proteins, F1094-0201 formed the most stable complex with BACE, followed by NMDA and GSK3β. The free energies of BACE-F1094-0201, GSK3β-F1094-0201, and NMDA-F1094-0201 complex formation were −73.78 ± 4.31 kcal mol −1 , −72.77 ± 3.43 kcal mol −1 , and −52.51 ± 2.85 kcal mol −1 respectively. In all the protein-ligand complexes, van der Waal forces, Coulombic forces, packing interactions, and lipophilic interactions favored the formation of a stable complex. Conversely, covalent interactions and solvation-free energy opposed the formation of a stable protein-ligand complex.

Computational Hardware and Software
Target proteins' three-dimensional (3D) crystallographic structures have been retrieved from the PDB database at http://www.rcsb.org/pdb/ (accessed on 12 December 2020). Through Autodock vina-enabled PyRx software, molecular docking was carried out, and the Lamarckian genetic method was used as a scoring function [76,77]. The visualization of molecular interactions was analyzed through the Biovia Discovery Studio visualizer [78]. The molecular dynamics simulation study was performed by Desmond (Schrodinger-2020, LLC, New York, NY, USA) software. An Intel Xenon (E3-1245-8C) workstation powered by NVIDIA Quadro P5000 with a 3.50 GHz processor and 28 GB RAM was used for the computational study.

Preparation of Ligands and Proteins
The 2029 natural product-like compounds library was retrieved in .sdf format from Life Chemicals (www.lifechemicals.com) as accessed on 2 November 2020. After that, all the compounds were energy minimized using a universal force field and converted to Autodock suitable ".pdbqt" format. Moreover, for the 3D structure of target proteins, we used GSK3β co-crystallized with a substrate (Adenosine-5 -Diphosphate) with a resolution of 2.10 Å (PDB Id: 1J1C), the high-resolution (1.90 Å) co-crystal structures of the NR1 S1S2 ligand-binding core with the antagonist 5,7-dichloro kynurenic acid (DCKA) well-known as an NMDA-receptor (PDB Id: 1PBQ), and BACE enzyme co-crystallized with hydroxyethyl amine inhibitor at 2.55 Å resolution (PDB Id: 1W51) were downloaded from the PubChem database [59,65,69]. After that, all the heteroatoms, including native ligands and water molecules, were removed from target proteins and hydrogens (polar only) were added. Then, geometric optimization and energy minimization of these edited target proteins were performed using PyRx built-in tool. The finalized proteins were later converted to the ".pdbqt" format.

Molecular Docking
Virtual screening of natural product-like compounds against target proteins was performed through PyRx-Python 0.8, a freely available tool coupled with AutoDock 4.2, as described earlier [8]. The ligands were docked individually after protocol validation by docking native ligands in the active site coordinates of the co-crystallized structure. The grid box of size 60 × 60 × 60 Å for target proteins of active site coordinates were finalized through the discovery studio visualizer at the attributes of native ligands docked in their specific proteins. The grid box was centered at 21.61 × 17.28 × −8.57 Å for GSK3β, at 5.61 × 38 × −17.16 Å for NMDA-receptor, and at 73.79 × 54.27 × 11.51 Å for BACE-1, respectively. The exhaustiveness was set to 8, and all other docking parameters were set to default values. The binding energy (∆G) of best poses was used to calculate the binding affinity (Kd) with the following equation [8]: In this equation, R stands for Boltzmann's gas constant, and T stands for temperature.

Molecular Dynamics Simulation
As previously reported, Desmond (Schrodinger-2020, LLC, New York, NY, USA) was employed to run a molecular dynamics simulation to evaluate the stability of target proteins with the best-shortlisted molecule after molecular docking and ADME study [8,66]. Briefly, the simulation was performed inside an orthorhombic box after placing the docked protein-ligand complex at the center at a 10 Å distance from the boundaries. TIP3P water molecules were used to solvate the box, and the appropriate Na + or Cl − counterions were added to neutralize it. To simulate physiological circumstances, 150 mM of sodium chloride was added. The system's energy was minimized through the OPLS3e forcefield by performing 2000 iterations with a convergence criterion of 1 kcal/mol/Å. Finally, a 100 ns production run was conducted at 298 K temperature and 1 bar pressure, which were kept constant using a Nose-Hoover Chain thermostat and a Martyna-Tobias-Klein barostat [80,81]. The energies and structures were recorded every 10 ps, with the time-step fixed at 2 fs. The molecular dynamics simulation results have been investigated for the root mean square deviation (RMSD), root mean square fluctuation (RMSF), the radius of gyration (Rg), solvent accessible surface area (SASA), secondary structure elements (SSE), and molecular interactions established between ligand and protein during simulation. The findings of each experiment were run independently, in triplicate, and are shown as mean ± standard errors.

Free Energy Calculations Using Prime/MM-GBSA
The free energy of interaction between protein and compound was calculated using Prime (Schrodinger-2020, LLC, New York, NY, USA) with the help of the MM-GBSA approach. Here, the free energies were calculated on the final 10 ns of the equilibrated trajectories. The docked protein-ligand complex was first optimized by molecular mechanics using Prime, followed by energy minimization by OPLS-AA forcefield using GBSA (Generalized Born Surface Area) continuum solvent model. Finally, the binding free energy of a protein-ligand complex was computed using the following equation: r ∆G Bind = ∆G Coulomb + ∆G vdW + ∆G Covalent + ∆G H−bond + ∆G Sol_Lipo +∆G Solv_GB + ∆G Packing + ∆G Sel f −contact where ∆G Bind , ∆G Coulomb , ∆G vdW , ∆G Covalent , ∆G H-bond , ∆G Sol_Lipo , ∆G Solv_GB , ∆G Packing , and ∆G Self-contact represent total free energy of protein-ligand complex, Coulombic free energy, van der Waals' interactions free energy, covalent bonds free energy, hydrogen bonds free energy, the free energy of the surface area, Solvation free energy, packing free energy, and free energy of self-contact respectively.

Conclusions
Out of 2029 natural product-like compounds, F1094-0201 displayed the best druglikeness, pharmacokinetics, and physiological properties that can cross the blood-brain barrier, and high GI absorption, based on the results of the virtual screening and molecular dynamics simulation (MDS) study. The MDS results indicated the establishment of a stable association between proteins and ligands based on the steady-state behavior of RMSD of all target proteins in the presence of F1094-0201. The RMSF plots of target proteins in the presence of F1094-0201 almost overlap with the RMSF plot of target proteins alone. Fluctuations in Rg and SASA of protein of interest in the presence of F1094-0201 did not diverge significantly. The total SSE of all the target proteins in the presence of F1094-0201 did not undergo any significant changes. These observations further emphasize that the F1094-0201 remained inside the binding pocket of target proteins and formed a stable protein-ligand complex. The free energies (MM/GBSA) results indicated that in all the protein-ligand complexes, van der Waal forces, Coulombic forces, packing interactions, and lipophilic interactions favored the formation of a stable complex. Amongst the target proteins, F1094-0201 formed the most stable complex with BACE, followed by NMDA and GSK3β. These attributes of F1094-0201 indicate it as a possible option for the management of pathophysiological pathways associated with AD. Our perspective is to utilize the beneficial effects of the lead compound (F1094-0201) against neurological disorders, and further, our strategy is to analyze the therapeutic potency of this compound by in-vivo experiments; thereafter, the drug candidate can be evaluated through clinical trials and get approval for patient care. In this study, the candidate drug's neuroprotective potentials were solely examined using in-silico methods, thereby opening the door for verification of its clinical effectiveness using in-vitro and in-vivo biological systems.