Inhibition of GSK_3β by Iridoid Glycosides of Snowberry (Symphoricarpos albus) Effective in the Treatment of Alzheimer’s Disease Using Computational Drug Design Methods

Introduction

Alzheimer’s disease (AD) is the most prevalent form of dementia (Lin, Jones et al., 2020). This health problem results in malignant neurological disorder along with cognitive, functional, and behavioral changes progressively and irreversibly. According to statistics, by 2050, this age-associated disease will affect 1 out of 85 people globally (Brookmeyer et al., 2007), posing a social and economic burden on the future community (Association 2019).

AD has a multifactorial etiology, with the most prominent ones associated with extracellular deposition of amyloid-β (Aβ) plaques, aggregation of tau protein, inflammation, oxidative stress, and declined levels of acetylcholine (De Strooper and Karran 2016; Hall et al., 2019). Nowadays, some evidence suggests that Aβ deposition is related to aging abnormalities. Also, tau (τ) protein is probably a higher target than Aβ because several therapeutic techniques concerning Aβ have not been highly effective in the disease treatment (Kametani and Hasegawa 2018). Tau protein is a member of the microtubule-related protein family that promotes microtubule stability and cytoskeleton formation (Kolarova et al., 2012). Intracellular tau protein aggregation in nerve fiber nodes is a neuropathological feature of AD. In the brain of people with AD, tau protein undergoes hyperphosphorylation because of some factors such as, formation of Aβ and environmental factors etc (Gong and Iqbal 2008). Thus, it is separated from microtubules (Mandelkow et al., 2007). Consequently, it results in neurofibrillary tangles (NFTs), ultimately leading to synapse loss and nerve cell destruction (Lauretti et al., 2020). Different kinases and phosphatases regulate the degree of phosphorylation of tau protein (Kolarova et al., 2012). In people with AD, tau protein is hyperphosphorylated, and it causes tauopathy due to degranulation and imbalance between them (Martin et al., 2013). One of the major kinases involved in tau protein phosphorylation is the glycogen synthase kinase_3 (GSK-3) enzyme (Jouanne et al., 2017). This serine–threonine kinase has two isoforms of alpha and beta (GSK3-α and GSK3-β) (Bagyinszky et al., 2014), which are encoded by two different genes (Lee et al., 2006).

In the central nervous system of people with AD, the GSK3-β enzyme is hyperactive, and evidence supports its role in AD pathology (Llorens-Marítin et al., 2014). GSK3-β phosphorylates many tau protein sites among people with AD and results in hyperphosphorylation of tau protein. Moreover, GSK3-β hyperactivity has been associated with impaired neurogenesis, microglia activation, amyloid-beta (Aβ) deposition, and inhibition of long-term potentiation (LTP) hippocampus and memory impairment (Llorens-Marítin et al., 2014).

This protein has two lobes, i.e., n-terminal and c-terminal, with residues 25 to 134 located in the N-terminal region and 35 to 380 in the c-lobe. Furthermore, these two lobes are held together by the hinge area such that residues 133 to 139 are considered a part of this area. Activation LOOP (A-Loop) is one of the most important loops of this enzyme where residues of A-LOOP are phosphorylated and have a key role in enzyme regulation. Among the residues of this loop, the DFG motif (ASP200-GLY202) is vital in enzyme regulation (ter Haar, Coll et al., 2001; Elangovan, Dhanabalan et al., 2020).

GSK3-β maintains an ATP-binding site between the two lobes. An important region of the ATP-binding site is the hinge region containing the residues Asp133, TYR134, and Val135. Some studies revealed that VAL135 and AP133 are very important residues for the binding of small molecules to the ATP-binding site (Gyurak et al., 2012; Pandey and DeGrado 2016; Davies 2019).

Hence, given the various roles of GSK3-β in AD pathology, the therapeutic potentials for its inhibition have been extensively examined in different studies.

Currently, there are no definitive treatments for AD (Eftekharzadeh et al., 2018), and the FDA-approved medicines for the disease are only for symptomatic treatment to improve behavioral changes and delay performance decline. These medicines include acetylcholinesterase inhibitors (e.g., donepezil, rivastigmine, and galantamine) and the N-methyl-D-aspartate receptor (NMDAR) (memantine). To date, several attempts have been made to find new inhibitors for AD, such as finding GSK3-β inhibitors.

The common snowberry plant (with the scientific name Symphoricarpos albus) has a fruit rich in phenolic acids, carbohydrates, and iridoids. The most important iridoid glycosides of this plant are loganin, loganetin, and secologanin (Gilbert 1995; Makarevich et al., 2009). Many biological activities like antioxidant, anti-inflammatory, anticancer, antidiabetic, antimicrobial, antiviral, anti-prosthetic, and neuroprotective activities have been reported for iridoid glycosides (Dinda et al., 2007; Dinda et al., 2019). In this respect, a study showed that iridoid glycosides effectively improve memory and cause nerve cell survival by increasing the expression of synaptophysin and neurotrophic factors in Alzheimer’s mice with cholinergic defects. Additionally, these compounds enhance cognitive impairment and inhibit hyperphosphorylation of tau protein in the hippocampus and striatum of SAMP8 mice (Ma et al., 2016).

Another study indicated that iridoid glycosides could increase pp2a activity, decrease GSK3-β activity, and finally inhibit tau hyperphosphorylation. Therefore, these treatment options can effectively improve cognitive and behavioral disorders in AD patients (Wang et al., 2020; Zhang et al., 2020). In this respect, we have used computational drug design methods, which are one of the fastest and most cost-effective ways available for drug design and finding new inhibitors of the GSK3- β enzyme. Moreover, we examined the inhibitory potential of the GSK3-β enzyme by loganin, loganetin, and secologanin with the help of molecular docking methods and molecular dynamic (MD) simulation.

Materials and Methods

Molecular Docking

Validation

Before molecular docking, we need to validate the software processing. Validation is of vital significance as it will eliminate false results. Since GSK3-β inhibition can be a treatment for AD, many inhibitors have been examined; many of them are competitive inhibitors with the ATP site and target its binding site in GSK3-β. The co-crystal structures of these inhibitors can be obtained through the Protein Data Bank (PDB) site. Codes related to the GSK3-β protein as a complex with the inhibitor are 1Q41 (Bertrand et al., 2003), 1Q3D (Bertrand et al., 2003), 1Q4L (Bertrand et al., 2003), 2OW3 (Zhang et al., 2007), 1R0E (Allard. 2004), 1Q3W (Bertrand et al., 2003), 3GB2 (Saitoh et al., 2009), 3F88 (Saitoh et al., 2009), 3F7Z (Saitoh et al., 2009), 3I4B (Aronov, Tang et al., 2009), 1UV5 (Meijer et al., 2003), 3Q3B (Coffman et al., 2011), 3L1S (Arnost et al., 2010), 1Q5K (Bhat et al., 2003), and 2O5K (Shin et al., 2007).

After downloading the PDB file corresponding to each of these codes, the structure of the GSK3-β protein and the cognate ligand of each code were prepared and optimized in ViewerLite 5.0 (Lite, 1998) software. All water molecules and nonpolar hydrogen of protein and nonpolar hydrogen of inhibitors were removed. In the next step, we performed redocking for each complex using both AutoDock4.2 (Forli et al., 2012) and AutoDock Vina (Trott and Olson 2010) software. We have used the two software to select the most appropriate software during the validation phase. Moreover, 100 runs were considered for each complex. Afterward, the best ligand conformation with the most negative ΔG and the lowest root mean square deviation (RMSD) was considered for each complex. Eventually, using Visual Molecular Dynamics (VMD) 1.9.3(http://www.ks.uiuc.edu) software, we obtained the RMSD between the selected conformations for the ligand and the ligand in the obtained crystalline structure from the PDB site.

In the next step, the obtained RMSDs of redock using two software were compared to each other for choosing the best software. Between these two tools, whichever showed the lower RMSD than the other, we chose it for loganin, secologanin, and loganetin docking with GSK3-β because it means that it can better predict the active site of the enzyme than the other.

Additionally, after selecting the software, the code with the lowest RMSD was used as the selected code for protein–ligand docking.

Docking Loganin, Secologanin, and Loganetin With GSK3-β Protein

Ligand Structure Preparation

First, the two-dimensional structure of each ligand was prepared in the ChemDraw Ultra 2d 8.0 (Mendelsohn 2004) tool and then was transmitted into ChemDraw Ultra 3d 8-0 (Mendelsohn 2004) software for preparing the three-dimensional structure. Then, the structures of ligands were subjected to energy minimization with the default of the MOPAC program in ChemDraw Ultra 3d 8.0 and then were saved as in the format of .mol. So, we converted them to the PDB format by ViewerLite 5.0 and used them as inputs of selected software for molecular docking.

The Gasteiger–Marsili procedure was used for calculating the partial charges of atoms (Shahlaei et al., 2011). Nonpolar hydrogens were deleted and then rotatable bonds were determined. All rotatable bonds of ligands were assumed to be flexible.

Protein Structure Preparation

We used the protein structure related to a selected complex in the validation step that it had lowest RMSD. And then, all water molecules and co-crystallized ligands were removed by the ViewerLite 5.0 tool, and just chain A of the protein was considered for molecular docking. Afterward, all hydrogens were added, and nonpolar hydrogens were removed by the selected tool. Then, Kolman charges were determined for the protein.

Docking Procedure

Molecular docking was performed by selected software in the validation stage. We used the Lamarckian genetic algorithm (LGA) method in AuotoDock4.2 based on the previous studies because they have revealed that other approaches (simulated annealing and genetic algorithm) are less effective than LGA (Morris et al., 1998). A total of 100 independent runs were considered for each ligand.

The dimensions of the center grid box for loganin, loganetin, and secologanin were considered 94.63 A° × 68.1680 A° × 9.7880 A°, and the software was allowed to search the entire volume of the protein. After finishing the docking, the two-dimensional figures of the results (e.g., the binding site and interactions between GSK3-β and inhibitors) were exhibited using LigPlot⁺ (Laskowski and Swindells 2011) software.

Molecular Dynamic Simulation

We carried out MD simulation in two steps. First, the structure of GSK3-β from a selected complex in the validation step was simulated in a water box. Second, to investigate the effects of ligand binding on GSK3-β conformation, the lowest binding energy conformation from the docking step of ligands was imported into the MD simulation.

For simulation, we used GROMAX 2019.1 (https://doi.org/10.5281/zenodo.2564761) software. At first, the protein was modified using the Swiss_PDB viewer (SPDBV) (Kaplan and Littlejohn 2001) program since the protein had the missing atom. Then, all the heteroatoms and water molecules were removed from the protein. In the next stage, topology and coordinate files were prepared. CHARMM27 was used as a force field, with the intermolecular (nonbonded) potential demonstration as a sum of Lennard–Jones (LJ) force and pairwise Coulomb interaction, and the long-range electrostatic force was defined by the particle mesh Ewald (PME) method. The velocity Verlet algorithm was applied for the numerical intermixture (Gharaghani et al., 2013).

Under periodic boundary conditions, the system was placed in a cubic water box of extended simple point-charge (SPC) water molecules with dimensions of 9.0465 × 9.06465 × 9.06465 nm³ with 1 nm away from the wall on each side.

Moreover, the system has to be examined for electrostatic neutrality, and anions or cations should be added if necessary. Here, 6 chloride ions were added to neutralize the system, and the entire system was included 5,579 atoms of GSK3-β, 6 Clˉ, and 22,560 solvent atoms. In the next stage, the energy was minimized in 50,000 steps for 2 fs. Then, the system was equilibrated at a constant temperature (NVT) of 300 k using the Berendsen thermostat. The cutoff radius was considered as 1.2 nm. Then, the system was balanced in the same way at a constant pressure (NVT) of 1 bar.

The goal of balancing the system at constant pressure and temperature for solvent molecules was to reach their best arrangement around the solute. Ultimately, the system was placed in the main run of the simulation for 100 ns at 300 K temperature and 1 bar pressure. To simulate the GSK3-β–ligand complex, first, the ligand-related topology file was generated using the SWISSPARAM server for the CHARMM27 force field. Afterward, the topology parameters and ligand coordinates were added to the topology parameters and coordinates of GSK3-β. Molecular dynamic simulation of the protein–ligand complex was performed similar to protein simulation for 100 ns. Ultimately, all the analyses related to the simulation were carried out in GROMACS 2019.1 software. These analyses include the evaluation of intermolecular hydrogen bonds, RMSD, Rg, RMSF, and free binding energy using the MM–GBSA method. The MD simulation was carried out on Ubuntu 18.04 Linux on an Intel Core 12 Quad 6800K 3.6 GHz, Gpu = gtx 1080ti NVIDIA, and 16 GB RAM.

Binding Free Energy Calculation

The binding free energy is a critical step of an in silico drug design approach which determines the binding affinity of inhibitors to the receptor. The binding free energy was calculated for the ligand receptor of loganin, loganetin, and secologanin complexes using MM/PBSA which defines the binding free energy of the protein and ligand as

Here, ΔG_complex demonstrates the total MMPBSA energy of the protein–ligand complex; ΔG_protein, and ΔG_ligand are solution free energies of the individual protein and ligand, respectively. The free energy of the individual existence can be expressed as follows:

E_MM shows the average molecular mechanical’s potential energy in the vacuum; G_solvation defines the free energy of solvation. T and S explain the temperature and entropy, respectively, and together TS exhibits the entropic contribution to the free energy in vacuum. In addition, the E_MM contains both bonded and nonbonded interactions of the molecules, including the bond angle, torsion, and electrostatic (E_elec) and van der Waals (E_vdw) interactions. Last, the free energy of solvation and G_solvation include both electrostatic and nonelectrostatic (G_polar and G_nonpolar) components. The binding free energies of the loganin_GSK3-β, loganetin_GSK3-β, and secologanin_GSK3-β complexes were calculated for 200 snapshots obtained from the last 20 ns of the trajectories (Padhi et al., 2021).

Drug Likeness Prediction and Toxicity

The OSIRIS property explorer (http://www.organic-chemistry.org/prog/peo/) was used to investigate the likeness of drugs. Pharmacokinetic properties like Log S calculation, TPSA, Clog P calculation, molecular mass, and toxicities like mutagenicity, tumorigenicity, irritation force, and hazard of the reproductive force of three compounds have been established.

Results and Discussion

Molecular Docking

In the molecular docking stage, it was seen that the obtained RMSD by AutoDock4.2 in most codes is less than the obtained results by AutoDock Vina software (Table 1). Since the AutoDock4.2 software can accurately predict the active site of GSK3-β protein, it was selected to dock our ligands. Moreover, according to the RMSD results of all codes by this software, the code 1UV5 with RMSD = 0.6679 had the lowest RMSD compared to other codes. 1UV5 contains the 6-Bromine dirubio-3-oxime inhibitor (BRW1383), which inhibits GSK3-β protein competitively with the ATP active site.

TABLE 1. Result of redock by AutoDock4.2 and AutoDock Vina software.

Among the most important amino acids with a crucial role in regulating this kinase, one can name ASP133, ASP200, and VAL135. As shown in Figure 1A, interacted BRW1383 in AutoDock 4.2 software with a −9.56 Kcal/mol binding energy could be matched on the BRW1383 from X-ray crystallography well and could establish the same bonds as those reported in the X-ray crystallographic structure (Figure 1B).

FIGURE 1. (A) Interacted BRW1383 (pink) and BRW1383 from X-ray crystallography (blue). (B) Interactions of BRW1383 with GSK3_β protein obtained from X-ray crystallography.

Thus, the docking method used in this study could effectively identify how the inhibitor binds to the enzyme. The average binding energy was, respectively, −7.15, −5.43, and −4.98 kcal/mol in the molecular docking of loganin, loganetin, and secologanin with GSK3-β protein. As can be seen, the binding energy of loganin was positive compared to that of BRW1383 but negative compared to that of other ligands. Docking results showed that loganin had a hydrogen bond with the amino acids VAL135 and ASP200 and a hydrophobic bond with ASP133. The details of the interaction between loganin and GSK3-β provided by Ligplot software are given in Figure 2B. As already stated, these three amino acids are the key amino acids in the regulation of GSK3-β protein (Figure 2A). However, BRW1383 has better binding energy vs. loganin, as shown in Table 2 which demonstrates that loganin interacts with all of the key amino acids in the active site, but the reference inhibitor interacts with just ASP133 and VAL135. As loganin is a natural compound, it can be a choice for inhibiting GSK3-β protein compared with BRW1383 which is a chemical compound. According to Figure 3, loganetin has established a hydrogen bond with VAL135 and a hydrophobic bond with ASP133. Additionally, secologanin has established two hydrogen bonds with ASN64 and ARG141 and several hydrophobic bonds outside the active site (Figure 4).

FIGURE 2. (A) Three-dimensional view of loganin at the GSK_3β protein binding site. (B) Interactions of loganin with amino acids of the GSK3_β protein binding site.

TABLE 2. Docking results of loganin, loganetin, secologanin, and BRW1383 with GSK3-β protein.

FIGURE 3. (A) Three-dimensional view of loganetin at the GSK_3β protein binding site. (B) Interactions of loganetin with amino acids of the GSK3_β protein binding site.

FIGURE 4. (A) Three-dimensional view of secologanin at the GSK_3β protein binding site. (B) Interactions of secologanin with amino acids of the GSK3_β protein binding site.

The results in Table 2 show that loganin binds more strongly to the active site of the GSK3-β protein compared to other ligands. To ensure the stability of the docked compound, the distance of hydrogen bonds between the compounds and the active site residues was investigated. According to Jeffrey’s study (Jeffrey 1997), the acceptable hydrogen bond distance between the donor and acceptor should be from 2.7 to 3.3 A°. The 2.2–2.5 Å distance is considered strong and covalent, 2.5–3.2 Å is considered medium and electrostatic bonds, and 3.2–4 Å is considered weak electrostatic bonds. Moreover, Ippolito and colleagues (1990) and Raschka and colleagues (2018) have revealed that the average distance between the protein and ligand from 2.4 to 3.5 Å could be considered hydrogen bonding.

Thus, the hydrogen bond interaction of our ligands with the protein was considered the medium hydrogen bond type. Therefore, in this study, the moderate hydrogen bond seems to be enough to inhibit the GSK3-β protein.

Molecular Dynamic Simulation

Although the docking results provide good information on the quality and interactions of the ligand and protein, they cannot predict how the protein conformations will change after ligand binding. Thus, to examine the changes and dynamics of the protein in interaction with three selected ligands, simulations were performed throughout for 100 ns

The before and after MD simulation of protein–ligand interaction residues are compared in Table 3 (Saddala and Adi 2018).

TABLE 3. The before and after simulation of protein–ligand interaction residues.

The active site amino acids of all compounds in the molecular docking were almost similar to the active site residues in the MD simulation. Thus, these compounds have stayed stable in the active site after simulation and did not change the protein structure. Although, secologanin before and after simulation could not interact with active site residues well and just could interact with ASP200. Consequently, secologanin was not considered a GSK3-β inhibitor. Also, we have prepared 3d and 2d demonstrations of all compounds after MD simulation as shown in Figures 5–7. As can be seen, the most important functional groups of loganin and loganetin that participate in hydrogen bond interactions with active site residues are related to the sugar part of these iridoid glycosides.

FIGURE 5. (A) 2d and (B) 3d view of the GSK_3β_loganin complex after MD simulation.

FIGURE 6. (A) 2d and (B) 3d view of the GSK_3β_loganetin complex after MD simulation.

FIGURE 7. (A) 2d and (B) 3d view of the GSK_3β_secologanin complex after MD simulation.

In molecular docking, only the ligand has a flexible state. However, both the protein and ligand were considered flexible in MD simulation. Therefore, during MD simulation, we could evaluate all observed interactions in molecular docking and find potential and new effects regarding the conformational changes of the ligand and the enzyme-binding site. As shown in Table 4, in hydrogen bond examinations of the GSK3-β_loganin complex, the most stable interactions have been established between O val135 and H20 atoms with 97% occupancy and HG1 THR138 with O4 loganin with 94% occupancy. Here, val135 is the active site amino acid and is involved in the regulation of GSK3-β protein. Among the hydrogen bonds that loganin has established with the protein, only the hydrogen bond with val135 is identical to the docking results. Although the protein is considered rigid in docking, in MD simulation, both the protein and ligand are considered flexible; because of that, MD results are more reliable. It is also noteworthy that only two hydrogen bonds were predicted for secologanin in docking, but several hydrogen bonds have been predicted for MD. As Table 4 shows, the occupancy listed for each of the secologanin hydrogen bonds is nonsignificant, and none of these bonds are established in the protein active site.

TABLE 4. Analysis of intermolecular hydrogen bonds.

Loganetin has two other hydrogen bonds, where its H17 atom is bonded to oxygen ASP133, and its O2 atom is bonded to HN VAL135, whereas these two amino acids are of the residues of the active site of proteins; however, the occupancy of the loganin hydrogen bonds is greater. Thus, loganin binds more strongly to the protein compared to the other two ligands.

Root Mean Square Deviation

The first critical analysis of MD simulation is the RMSD, which shows the stability and structural changes during the simulation. RMSD is a parameter that shows the deviation of the particle position from the original structure at any point in time (Czelen 2017; Sargsyan et al., 2017).

The lower RMSD value indicates less fluctuation during the simulation, suggesting that the stability of the protein is high (Elangovan et al., 2020). Therefore, the system will be balanced. Figure 6 presents the variation in RMSD values of GSK3-β protein comparing with the three complexes. The RMSD value revealed that all pathways reached equilibrium after 10 ns Thus, RMSF, Rg, SASA, and hydrogen bonds were analyzed at an equilibrium state.

As Figure 8 shows for the GSK3-β protein (green diagram), the mean RMSD has remained at about 0.19 nm with minimal fluctuation. Consequently, the protein does not undergo severe deformation and denaturation while simulating.

FIGURE 8. RMSD analysis GSK3-β protein (green), loganin-GSK3-β complex (yellow), loganetin-GSK3-β complex (blue), and secologanin-GSK3-β complex (red).

As can be seen from the GSK3-β_secologanin complex, the RMSD has higher values and fluctuations (0.23 nm) relative to the free protein and compared to loganin and loganetin about 0.19 and 0.21, respectively. Thus, the protein–secologanin complex is more unstable than other structures. Here, the protein–loganin RMSD compared to GSK3-β protein is more stable than the other two complexes. Also, the enzyme structure has not changed significantly in the presence of loganin, and the system is stable during the simulation.

Root Mean Square Fluctuation

Another used parameter to measure stability and flexibility is the RMSF during 90 ns. The RMSF estimates the average deviation of a particle (like protein residues) from a reference position (typically the average particle position) over time (Reddy et al., 2015).

Therefore, RMSF analyzes those parts of the structure with the most and least fluctuations from their dependent structure, and this parameter describes how ligand binding can result in conformational changes at the residual level (Elangovan et al., 2020; Saravanan et al., 2020).

Overall, the peaks seen in the RMSF diagram showed the most unstable target protein residues.

Figure 9 illustrates the high residue fluctuations at the start and end of the RMSF chart.

FIGURE 9. RMSF analysis GSK3-β protein (green), loganin-GSK3-β complex (yellow), loganetin-GSK3-β complex (blue), and secologanin-GSK3-β complex (red).

Moreover, from Figure 9, loganin has less fluctuation compare with loganetin and secologanin, and it has the most adaption with the RMSF GSK3-β protein diagram.

Moreover, amino acids' RMSF of the active site in each of the three ligands reveals that binding of the ligand to the protein in the active site does not cause the fluctuation of residues, and the structure remains stable. The mean RMSF values for GSK3-β protein, loganin, loganetin, and secologanin were 2.16, 2.16, 2.17, and 2.15, respectively. These values suggest that ligand binding does not change the original conformation of the residues. Also, it is inferred that the compounds do not fluctuate much, they seem consistent with protein fluctuations, and the complexes seem stable.

Radius of Gyration

The gyration radius is a variable that is studied to examine the compression changes during MD simulation. Protein compression during interactions with the ligand is affected by protein chains, in which the flexibility between the ligand and the protein depends on it (Shukla et al., 2019; Ghosh et al., 2020).

The calculated Rg value for the three complexes (Figure 10) was 2.16, which is associated with the protein size as well. The Rg of GSK3-β protein was 2.16, indicating that GSK3-β protein remains compressed during MD simulation and binding to all three ligands.

FIGURE 10. Rg analysis GSK3-β protein (green), loganin-GSK3-β complex (yellow), loganetin-GSK3-β complex (blue), and secologanin-GSK3-β complex (red).

Solvent-Accessible Surface Area

The SASA value is used to analyze the magnitude and significance of ligand binding to the receptor and the changes in protein conformation due to ligand binding (Shukla et al., 2019; Elangovan et al., 2020).

As shown in Figure 11, all complexes show a similar magnitude of protein conformation that interacts with the ligand compared to the GSK3-β protein.

FIGURE 11. SASA analysis GSK3-β protein (green), loganin-GSK3-β complex (yellow), loganetin-GSK3-β complex (blue), and secologanin-GSK3-β complex (red).

On the other hand, SASA changes are similar to the changes in the Rg diagram, confirming the accuracy of the molecular dynamic simulations. Thus, one can infer that all protein–ligand complexes are stable in the SASA analysis.

The Analysis of Hydrogen Bonds

Different interactions like hydrogen bonds, hydrophobic interactions, and ionic interactions stabilize the protein–ligand complex. Between them, the hydrogen bonds are more important and specific transient interactions than others for protein–ligand stabilization (Shukla et al., 2019). The number of hydrogen bonds was calculated for the last 90 ns trajectory. The number of hydrogen bonds vs. time is revealed for each complex in Figure 12. The average number of hydrogen bonds for loganin_GSK3-β (Figure 12A) and loganetin_GSK3-β was (Figure 12B) 3 and 2, respectively. But the hydrogen bonds of secologanin_GSK3-β are variable.

FIGURE 12. H-bond analysis (A) loganin, (B) loganetin, and (C) secologanin.

As can be seen, a partial change in the hydrogen bond formation between the ligand and protein for loganin_GSK3-β and loganetin_GSK3-β was observed. So, these complexes were stable for most parts of the simulation trajectory compared to secologanin_GSK3-β. These results are approximately similar to the results that are presented in Table 4.

Binding Free Energy Calculation

The binding free energy is a critical step of the in silico drug design approach which determines the binding affinity of inhibitors to the receptor. The binding free energy was calculated for the ligand receptor of loganin, loganetin, and secologanin complexes using the MM/PBSA method. The nonpolar and polar solvation energy was analyzed in the electrostatic interaction, van der Waals energy, and SASA energy. The major desirable portions of the stereoisomer binding were the van der Waals and electrostatic energies for all complexes. Besides, the terms polar solvation free energy and SASA energy are unfavorable for binding in three complexes. The total free binding energy (ΔG_Binding) was calculated for each compound in Table 5.

TABLE 5. Table represents the Van der Waals, electrostatic, polar solvation, SASA, and binding energy in kJ.mol_1 for secologanin, loganetin, and loganin.

Loganin and loganetin have an almost similar amounts of binding energy, respectively, −52.421 and −52.878 KJ/mol, but they are much higher than that of secologanin (−1.870 KJ/mol). We have concluded that loganin and loganetin are more energetically favorable than secologanin, and these results are adopted with docking conclusions.

Drug Likeness Prediction and Toxicity

Pharmacokinetic properties and toxicities were predicted for each compound used by the OSIRIS property explorer server, and the results are revealed within Table 6. The mutagenicity, tumorigenicity, irritation force, and hazard of the reproductive force were predicted for toxicity verification.

TABLE 6. Drug likeness prediction and toxicity of loganin, loganetin, and secologanin.

The higher cLog P (logarithm of compounds' partition coefficient between water and n-octanol) value indicates lower permeation and absorption due to lower hydrophilicity. The solubility and molecular weight (MW) affect the absorption rate; thus, high solubility and lower MW increase absorption. The topological polar surface area (TPSA) reveals the surface of polar atoms in the compounds. An increased TPSA value is relevant to the least permeability of the membrane. So, the compounds with a larger TPSA value will be a better substrate for p-glycoprotein which is amenable for the efflux of a drug from the cell, and thus, the diminished TPSA was useful for the drug-like property. Some investigations also predicted that a molecule with a better CNS penetration should have a lesser value of the TPSA (Reddy et al., 2015; Srivastava et al., 2020).

Furthermore, the drug score mixes drug likeness, clog P, TPSA, MW, and toxicity risk parameters to reveal a compound as a drug candidate. Based on Table 6, loganetin has a higher score (0.71) than loganin (0.45) and secologanin (0.27) because loganin has a lower TPSA, MW, and cLog P and higher solubility than other compounds.

Conclusion

This study aimed at finding potential and available inhibitors for GSK3-β protein using computational drug design methods like molecular docking and MD simulation.

We performed molecular analysis of the iridoid glycosides in the Common snowberry plant, including loganin, secologanin, and loganetin, which have beneficial effects on improving memory.

Molecular docking analysis revealed that loganin and loganetin bind exactly to the ATP-binding site. Loganin established hydrogen bonding with VAL135 and ASP200 and hydrophobic bonding with ASP133, and loganetin established hydrogen bonding with VAL135 and hydrophobic bonding with ASP133, but the binding affinity of loganin was more than that of loganetin.

RMSD, RMSF, Rg, SASA, hydrogen bond analyses, and MMPBSA were carried out in the molecular dynamic simulation. The results of the RMSD value indicated that the loganin–GSK3 complex has more stability than the other two complexes.

The RMSF diagram indicated that ligand binding does not increase fluctuations, and the complexes stay stable following the binding.

The SASA analysis was equal with variations of the Rg value, and equality of the Rg value of the GSK3-β protein with complexes indicates that the structure remains compact during MD simulation.

Loganin and loganetin have an almost similar amount of binding energy, and loganetin has a better drug score than others.

Consequently, both loganin and loganetin may effectively inhibit GSK3-β because these compounds established strong hydrogen bonds with the active site residues before and after MD simulation. Although the binding affinity of loganin and loganetin was slightly less than that of the BRW1383 inhibitor in molecular docking, they can be considered to prevent tauopathy since those are available and nontoxic herbal compounds. Nonetheless, this study requires more examination, and it is better to carry out laboratory studies along with molecular studies to prove the effects of iridoid glycosides on AD prevention and treatment.

Data Availability Statement

The original contributions presented in the study are included in the article/Supplementary Materials, further inquiries can be directed to the corresponding author.

Author Contributions

ME and SP analyzed data and were responsible for manuscript writing. PK was the advisor and assisted in gathering information and research. BA assisted in research and writing the primary version of the manuscript. SG and JK were the research supervisors. All authors read and approved the final manuscript.

Conflict of Interest

BA is the president and CEO at CaroGen Corporation in Farmington, CT United States of America.

The remaining authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Acknowledgments

We gratefully appreciate the support of the Faculty of Pharmacy of Lorestan University of Medical Sciences for providing the facilities to carry out this study.

References