Identification of the most potent acetylcholinesterase inhibitors from plants for possible treatment of Alzheimer’s disease: a computational approach

Being one of the rapidly growing dementia type diseases in the world, Alzheimer’s disease (AD) has gained much attention from researchers in the recent decades. Many hypotheses have been developed that describe different reasons for the development of AD. Among them, the cholinergic hypothesis depicts that the degradation of an important neurotransmitter, acetylcholine by the enzyme acetylcholinesterase (AChE), is responsible for the development of AD. Although, many anti-AChE drugs are already available in the market, their performance sometimes yields unexpected results. For this reason, research works are going on to find out potential anti-AChE agents both from natural and synthetic sources. In this study, 50 potential anti-AChE phytochemicals were analyzed using numerous tools of bioinformatics and in silico biology to find out the best possible anti-AChE agents among the selected 50 ligands through molecular docking, determination of the druglikeness properties, conducting the ADMET test, PASS and P450 site of metabolism prediction, and DFT calculations. The predictions of this study suggested that among the selected 50 ligands, bellidifolin, naringenin, apigenin, and coptisine were the 4 best compounds with quite similar and sound performance in most of the experiments. In this study, bellidifolin, naringenin, apigenin, and coptisine were found to be the most effective agents for treating the AD targeting AChE. However, more in vivo and in vitro analyses are required to finally confirm the outcomes of this research.


Background
First described by Alois Alzheimer in 1907, Alzheimer's Disease (AD) has become one of the most prevalent dementia type diseases in the world which is increasing its numbers rapidly [1,2]. Intellectual morbidity, psychomotor dysregulation, delusions, hallucinations etc., are some of the familiar symptoms of AD [3]. In the familial and congenital cases of AD, genetic factors play critical roles [4]. Different hypotheses have been developed by the scientists that shed light on several reasons for AD onset and development. One of these hypotheses is the amyloid cascade hypothesis, which describes that the deposition of β-amyloid plaques in the brain is mainly responsible for the development of AD which is the result of abnormal processing of the amyloid precursor protein (APP) by the β-secretase enzyme. These plaques interfere with the normal functions of the brain [5]. On the other hand, another hypothesis called the oxidative stress hypothesis describes that because of the deposition of increased amounts of ions like iron, aluminium, and mercury, free radicals and reactive oxygen species (ROS) are generated very rapidly in the brain which are responsible for increased lipid peroxidation and protein and DNA oxidation. The stresses produced by these oxidation events lead the way for the AD onset [6]. There is another hypothesis of AD development which is known as the cholinergic hypothesis. According to this hypothesis, the loss of functions of the cholinergic neurons and thus the cholinergic signaling and neurotransmission in the brain maybe responsible for AD [7]. Table 1 lists the current status of therapeutic agents that are intended to or being used to alleviate the complications related to AD. This experiment was conducted focusing on the cholinergic hypothesis of AD development.

The cholinergic hypothesis and the development of Alzheimer's disease
The cholinergic hypothesis concerns with one of the major neurotransmitters, acetylcholine (ACh) which is regulated by two enzymes, acetylcholinesterase (AChE) and choline acetyltransferase (ChAT) [13]. ACh is involved in many important functions of the brain like learning and memory generation processes. It performs its functions through binding to two types of receptors, i.e., nicotinic (α7 and α4β2) and muscarinic receptors (M1 muscarinic receptor). ACh is synthesized by the enzyme choline acetyltransferase (ChAT) which catalyzes the transfer of an acetyl group from acetyl coenzyme A (Ac-CoA) to choline (Ch) in the pre-synaptic neuron and thus synthesizes the ACh. Thereafter, the ACh is secreted by the pre-synaptic neuron into the synapse where it mediates its effects by binding to either the nicotinic receptor or the muscarinic receptor. To maintain the optimal concentrations of ACh required for proper functioning of the brain, another enzyme called acetylcholinesterase (AChE) is synthesized by a serine hydrolase enzyme which hydrolyzes ACh to acetate and Ch. Then the Ch is again taken up by the pre-synaptic neuron for recycling and reusing. Thus, the balance of ACh is maintained in the normal brain. However, there is evidence that, in the brain of AD patients, the overexpression of AChE occurs. This phenomenon decreases the amount of ACh required for proper functioning of the brain which is why the neuron cells cannot operate properly and complications like brain damage as well as memory loss occur. These complications lead to the onset of AD development (Fig. 1) [14][15][16][17][18].
Donepezil and rivastigmine are two FDA-approved drugs that are currently used for mild to moderate AD treatment targeting the AChE enzyme. However, both of them have several side effects like nausea, diarrhea, anorexia, syncope, abdominal pain, and vomiting [19,20]. Therefore, scientists are searching for more effective agents that can provide more efficacy than these drugs with much lesser side effects. Scientists have also focused on the natural resources for potential anti-AChE agents since the natural agents are generally much safer than synthetic chemicals. Galantamine is such a natural drug isolated from Galanthus woronowii which is also currently used for AD treatment alongside other approved chemical drugs. But since none of these drugs are found to be quite satisfactory to stop the progression or development of AD, research is going on to find out new compounds from natural sources with anti-AChE properties [21][22][23].
Molecular docking is a widely accepted and used technique in drug R&D which reduces both time and costs of lead discovery processes. This method is also known as computational drug design which has already been used for designing over 50 novel drugs, and many of them have also gained FDA approval for marketing. By simulating the interaction between ligand and receptor in the computer software, the docking system assigns scoring functions to the bound ligands which reflect their binding affinity. The lower docking score represents the greater binding affinity and vice-versa [24,25]. The current study was designed to predict the best ligands among 50 selected phytochemicals with potential anti-AChE activities based on the molecular docking analysis (Table 2). Thereafter, the pharmacodynamics and physicochemical characteristics of the best selected ligands were predicted by determining their druglikeness properties, conducting the ADMET test, PASS, and P450 site of metabolism prediction and DFT calculations.

Methods
Total 50 phytochemicals were selected as ligands in this study by reviewing numerous literatures along with their During protein preparation, the bond orders were assigned and hydrogen molecules were added to heavy atoms as well as all the waters were deleted and the side chains were adjusted using Prime [61]. After that, the structure was optimized and minimized using force field OPLS_2005, which was conducted setting the maximum heavy atom RMSD (root-mean-square-deviation) to 30 Å and any remaining water less than 3 H-bonds to non-water was again deleted during the minimization step.

Ligand preparation
Three-dimensional structures of the 50 selected ligand molecules as well as the control were downloaded from PubChem database (www.pubchem.ncbi.nlm.nih.gov). These structures were prepared for docking using the Lig-Prep module of Maestro Schrödinger Suite [62]. Minimized 3D structures of ligands were generated by Epik2.2 using OPLS_2005 force field and within pH 7.0 ± 2.0 [63].

Receptor grid generation
Grid confines the active site of the receptor protein to a shortened specific area for the ligand to dock specifically. In glide, a grid was generated where the default Van der Waals radius scaling factor 1.0 and charge cutoff 0.25 was used. The grid was then subjected to OPLS_2005 force field. A cubic box was generated around the active site (reference ligand active site). Then, the grid box volume was adjusted to 14 × 14 × 14 for docking test.
Glide standard precision (SP) and extra precision (XP) ligand docking, Prime MM-GBSA calculation and induced fit docking SP and XP adaptable glide dockings were conducted using the Glide module in Maestro Schrödinger Suite [64]. The Van der Waals radius scaling factor and charge cutoff were kept at 0.80 and 0.15, respectively, for all the ligand molecules. Final score was assigned by the module by analyzing the pose of docked ligand within the active site of the receptor. After SP and XP ligand docking, the docked complexes were subjected to molecular mechanics-generalized born and surface area (MM_GBSA) rescoring with the help of Prime module of Maestro Schrödinger suite for further evaluation. This technique utilizes the docked complex and uses an implicit solvent that assigns more accurate scoring function and improves the overall freebinding affinity score upon the reprocessing of the complex. It combines OPLS molecular mechanics energies (E MM ), surface generalized born solvation model for polar solvation (G SGB ), and a nonpolar salvation term (G NP ) for total free energy (ΔG bind ) calculation. The total free energy of binding was calculated by the following equation [65]: The agents with best results in the SP and XP docking were selected for the MM_GBSA and IFD studies.
Thereafter, to carry out the IFD of the selected ligand molecules from SP and XP docking, again OPLS_2005 force field was applied after generating grid around the co-crystallized ligand of the receptor, and this time, the best five ligands were docked rigidly. Receptor and ligand Van Der Waals screening was set at 0.70 and 0.50, respectively, and residues within 2 Å were refined to generate 2 best possible posses with extra precision. Four best performing ligands were selected according to their MM_GBSA score, IFD score, and XP Gscore . The 3D representations of the best pose interactions between the best four ligands and their respective receptors were obtained using Discovery Studio Visualizer [66]. At these stages, the docking parameters of the compounds under investigation were compared with the control

Ligand-based drug likeness property and ADMET prediction
The drug likeness properties of the 4 selected ligand molecules were analyzed using SWISSADME server (http://www.swissadme.ch/) [67]. After that, the ADME T for each of the ligand molecules was conducted using the online based server, ADMETlab (http://admet.scbdd. com/) to predict various pharmacokinetic and pharmacodynamic properties [68,69]. The numeric and categorical values of the results generated by the ADME Tlab server were converted into qualitative values according to the documentation and explanation described in the ADMETlab server (http://admet.scbdd.com/ home/interpretation/) for the convenience of interpretation.

PASS and SOM prediction
The PASS (Prediction of Activity Spectra for Substances) prediction of the best four selected ligands were carried out by the PASS-Way2Drug server (http://www. pharmaexpert.ru/passonline/), using the canonical SMIL ES from PubChem server (https://pubchem.ncbi.nlm.nih. gov/) [70]. While carrying out PASS prediction, the P a (probability to be active) was kept greater than 70%, since the P a > 70% threshold generates highly reliable prediction [71]. In the PASS prediction study, 15 possible biological activities were predicted. The P450 Site of Metabolism (SOM) of the four best selected ligand molecules was determined by online tool, RS-WebPredictor 1.0 (http://reccr.chem.rpi.edu/Software/ RS-WebPredictor/) [72]. Moreover, the LD50 value and toxicity class of the compounds were predicted using the ProTox-II server (http://tox.charite.de/protox_II/) [73].

DFT calculations
Minimized ligand structures obtained from LigPrep were used for DFT calculation using the Jaguar panel of Maestro Schrödinger Suite using Becke's threeparameter exchange potential and Lee-Yang-Parr correlation functional (B3LYP) theory with 6-31G* basis set [74][75][76]. Quantum chemical properties such as surface properties (MO, density, potential) and multipole moments were calculated along with HOMO (Highest Occupied Molecular Orbital) and LUMO (Lowest Unoccupied Molecular Orbital) energy. Then, the global frontier orbital was analyzed and hardness (η) and softness (S) of selected molecules were calculated using the  [77,78]. The DFT calculation was done for the 3 best ligand molecules.

Molecular docking study
All the 50 selected ligands were docked successfully with their receptor protein, AchE. The ligand molecules that had the lowest binding energy were considered the best ligand molecules because lower binding energy or docking score represents higher binding affinity [79]. In the MM-GBSA study, the most negative ΔGBind score is also considered as the best ΔGBind score [80]. IFD study was conducted to determine the accurate binding mode and accuracy of active site geometry. The lowest values of IFD score and XP GScore represent the best results [81][82][83][84]. From the initial 50 ligands, the ligands that had both SP and XP scores of less than − 8.00 Kcal/mol were selected for MM_GBSA analysis and IFD study. Thirteen ligands were found to have both SP and XP scores less than − 8.00 Kcal/mol: bellidifolin, isoharmnetin, quercetin, myricetin, imperatorin, naringenin, harmaline, harmine, isoimperatorin, apigenin, coptisine, liriodenine, and scoulerine (Table 3). From these 13 ligands, 4 ligands were finally selected based on their lower scores in both MM_GBSA (ΔG bind score) and IFD studies (XP G score and IFD score), i.e., bellidifolin, naringenin, apigenin, and coptisine (Table 4). These four ligands were selected for further analysis to predict and determine their drug potentials. When compared to the two positive controls: donepezil and galantamine, it can be declared that these four compounds generated sound results in the molecular docking study.

Binding mode of best ligands with respective targets
The 3D representations as well as the interaction of different amino acids of the binding pocket of AChE with bellidifolin, naringenin, apigenin, and coptisine are illustrated in Fig. 2

Druglikeness properties
The druglikeness property experiment was conducted for only the 4 best selected ligands: bellidifolin, naringenin, apigenin, and coptisine. All the 4 ligands were predicted to follow the Lipinski's rule of five. Moreover, bellidifolin was found to have the highest topological polar surface area (TPSA) value of 100.13 Å 2 . Again, coptisine was predicted to have the highest molar refractivity of 87.95 as well as the highest concensus Log P o/w of 2.40 (Table 5). All the found compounds generated quite sound performance in the analysis when compared with the two positive controls, donepezil and galantamine.

ADMET prediction
The results of the ADMET test are summarized in Table  6 along with their comparison with donepezil and galantamine. All the four ligands performed quite similarly in the ADMET test. In the absorption section, all the ligands were predicted to be optimal Caco-2 permeable, whereas all of them were P-gp (P-glycoprotein) noninhibitors and P-gp non-substrates. In the distribution section, apigenin showed relatively poor performance with no capacity to be blood-brain barrier (BBB) permeable. In the metabolism section, all of the ligands were predicted to be CYP450 1A2 inhibitor and only apigenin was found to be CYP450 1A2 non-substrate. However, all of them were CYP450 2C9 and CYP450 2C19 noninhibitors and bellidifolin and apigenin were substrates for CYP450 2C9. Moreover, only bellidifolin and coptisine were found to be substrates for CYP450 2C19 and CYP450 2D6. In the excretion section, coptisine was found to have the highest half-life of 1.8 h. In the toxicity section, all the four ligands were found to be drug induced liver toxicity (DILI) positive and only bellidifolin was Ames positive as well as hERG non-blocker. Moreover, both naringenin and coptisine were found to have human hepatotoxicity. The performance of the 4 finally selected ligands was not very different from that of the two positive controls.

PASS and P450 site of metabolism prediction
The predicted LD50 value for bellidifolin, naringenin, apigenin, and coptisine were 4000 mg/kg, 2000 mg/kg, 2500 mg/kg, and 1000 mg/kg, respectively. Moreover, both bellidifolin and apigenin were predicted to be in toxicity class 5, whereas naringenin and coptisine were found to be in toxicity class 4. The prediction of activity spectra for substances (PASS prediction) was conducted for the 4 ligands to identify 15 intended biological activities. The PASS prediction results of the 4 selected ligands along with the positive control donepezil are listed in Table 7. Coptisine did not show any of the activities; however, both bellidifolin and apigenin were

Analysis of frontier's orbitals
In the analysis of Frontier's orbitals, the DFT calculations were performed. The results of the DFT calculations are listed in Table 9. Coptisine showed the lowest gap energy of 0.047 eV as well as the highest dipole moment of 6.459 debye. On the other hand, apigenin generated the highest gap energy of 0.130 eV and the second highest dipole moment of 2.163 debye. The order of dipole moments of these 4 compounds were, bellidifolin < apigenin < naringenin < coptisine. Comparing with the positive controls, it can be declared that all the four compounds showed sound results in the DFT calculations (Fig. 3).

Discussion
Molecular docking is an effective strategy in computer-aided drug designing which works on some specific algorithms and assigns affinity scores depending on the poses of the ligands inside the binding pocket of a target. In molecular docking, the lowest docking score corresponds to the highest affinity which reflects that the complex remains more time in contact with good stability [85,86]. In this study, 50 compounds that were proven to have anti-AChE properties were selected by literature reviewing. Initially, from these 50 compounds, 13 compounds were selected based on their SP and XP docking scores (compounds having both SP and XP docking scores over − 8.0 Kcal/mol). When compared with the positive controls, donepezil and galantamine, all the compounds were found to have relatively good results in molecular docking analysis. The 13 ligands selected from the initial docking analysis had predicted results better than donepezil. Later, from these 13 compounds, according to the MM_GBSA and IFD scores, 4 compounds, i.e., bellidifolin, naringenin, apigenin, and coptisine, were selected as final ligands for further analysis. These final compounds were also predicted to have better results than donepezil and galantamine in the MM_GBSA and IFD studies.
AChE possesses several notable amino acids in its active-site gorge, i.e., serine 200, glutamic acid 327, histidine 440, tryptophan 279, tyrosine 121, phenylalanine 330, and tryptophan 84. All the 4 finally selected ligands were found to form potential hydrogen and hydrophobic interactions with histidine 440, phenylalanine 330, and tryptophan 84 within the active site of AChE. Since hydrogen and hydrophobic interactions are important for strengthening the receptor-ligand interactions, it can be declared that these ligands might be able to effectively inhibit the AChE enzyme at its active site [87,88].
The prediction of druglikeness properties facilitates the drug development processes. The molecular weight and TPSA are two important characteristics of a drug that affect its permeability through the biological barriers. Higher molecular weight and TPSA represent lower permeability of the drug molecule and vice-versa. Again, the lipophilicity (expressed as LogP) influences the absorption of the drug molecule in the body and the higher LogP of a drug represents its lower absorption in the body. The capability of a drug molecule to cross the cell membrane is also influenced by the number of hydrogen bond acceptors and donors beyond the acceptable range. Moreover, the number of rotatable bonds (acceptable range is less than 10) also influences the    [89][90][91][92]. All the 4 selected ligands were found to follow the five Lipinski's rules of druglikeness properties. All of them showed quite good results in the druglikeness property experiment and for this reason, and all of them were considered as potential drug candidates. ADMET prediction of a candidate drug aids in determining the pharmacological and pharmacodynamic properties of the drug molecule within a biological system. For this reason, the ADMET properties are important determinants for the success of a drug development process. BBB permeability is the most important characteristic of the drugs which primarily target the brain cells. Again, the transportation of drugs is influenced by the P-glycoprotein (P-gp) in the cell membrane; therefore, inhibition of this protein severely affects the drug transport. Moreover, if a drug is found to be permeable in the Caco-2 cell line, then that particular drug is considered to be absorbed well in the intestine. Orally absorbed drugs travel through the blood circulation and deposit in the liver where they are metabolized by group of enzymes of the Cytochrome P450 family and later excreted as bile or urine. Therefore, inhibition of any of these enzymes affects the metabolism and biodegradation of the drugs [93,94]. Furthermore, if a drug is found to be a substrate for one or more CYP450 enzyme or enzymes, then that drug is considered to be metabolized well by the respective CYP450 enzyme or enzymes [95]. A drug's pharmacodynamic properties depend on the degree of its binding with the plasma protein. A drug can pass through the cell layers if it binds to the plasma proteins less efficiently and vice versa [96]. Human intestinal absorption (HIA) is another crucial determinant for the orally administered drugs [97][98][99]. The half-life of a drug describes that the greater the half-life, the longer it would stay in the body and the greater its potentiality [100][101][102].
Furthermore, HERG is a K+ channel which is found in the heart muscle. Blocking the hERG channel or signaling may cause cardiac arrhythmia leading to death [103,104]. On the other hand, human hepatotoxicity (H-HT) refers to any type of acute liver injury which may lead to organ failure and death [105,106]. The Ames test is a mutagenicity assay that is used to detect the potential mutagenic chemicals [107]. Drug induced liver injury (DILI) reflects the injuries in the liver that are caused by administration of drugs [108]. Bellidifolin, naringenin, and coptisine showed quite similar and sound performance in the AMDET prediction. However, apigenin was found to be BBB impermeable and also a good binder to the plasma membrane, which could affect its biodistribution throughout the body as a potential inhibitor of AChE. As a consequence, the performance of apigenin in ADMET test was declared to be non-satisfactory.
Prediction of Activity Spectra for Substances or PASS prediction is a process of determining the biological activities of drug candidates or drug-like molecules. The PASS prediction estimates biological activities by predicting two potential probabilities: Pa, which represents the probability of a compound to be active and Pi, which represents the probability of a compound to be inactive. By predicting these two probabilities, the PASS prediction determines potential biological activities of a drug-like molecule. The Pa value of greater than 0.7 represents the greater probability of exhibiting the activity of a substance in an experiment. The Pa > 70% threshold corresponds to a highly reliable prediction, while Pa > 0.5 gives quite reliable prediction and the Pa > 0.3 threshold describes that a compound is less likely to have the desired activity [70,71,109]. Most of the intended activities (13 activities) were common between bellidifolin, naringenin, and apigenin, for example, antineoplastic, antimutagenic, membrane integrity agonist, HIF1A expression inhibitor, and apoptosis agonist ( Table 8).
The toxicity classes of the 4 compounds were determined by the ProTox-II server (http://tox.charite.de/ protox_II/index.php?site=compound_input), which predicts the toxicity of a compound and then classifies the compound into a toxicity class that ranges from 1 to 6. To classify the compounds into toxicity classes, the server uses the Globally Harmonized System of Classification and Labelling of Chemicals (GHS). According this system, since both bellidifolin and apigenin had toxicity class 4, they would be harmful if swallowed. Again with the toxicity class of 5, both naringenin and coptisine might be harmful if swallowed [73,110]. Both the positive controls, donepezil and galantamine were found to have the toxicity class of 4.
The cytochrome P450 or CYP450 represents a family of 57 isoforms of the heme-containing P450 enzymes. This enzyme family is one of the most important enzyme family for proper metabolism within the human body, which is responsible for catalyzing the phase-I metabolism of about 90% of all the marketed drugs. They also convert the lipophilic drugs or compounds to more polar, hydrophilic substances. Among these 57 isoforms of the CYP450 enzymes, 9 isoforms are most important for biological activities, i.e., CYPs 1A2, 2A6, 2B6, 2C19, 2C8, 2C9, 2D6, 2E1, and 3A4 [111,112]. When the P450 SOMs of the 4 best selected ligands: bellidifolin, naringenin, apigenin, and coptisine were predicted, all of them showed potential SOMs where these 9 isoforms of P450 enzymes can act on. For this reason, all these ligands or compounds might be metabolized well by the body.
Frontier's orbitals study or DFT calculation is an important method to predict the pharmacological properties of small molecules. HOMO and LUMO help to study the chemical reactivity and kinetic stability of the query molecules. The term "HOMO" represents the regions on a small molecule which might donate electrons during the formation of a complex and the term "LUMO" corresponds to the regions on a molecule that might receive electrons from the electron donor. The difference in HOMO and LUMO energy is regarded as gap energy that corresponds to the electronic excitation energy. The molecule that has greater orbital gap energy tends to be energetically more unfavorable to undergo a chemical reaction and vice versa [76,[113][114][115][116][117]. All of the 4 ligands were predicted to have significant gap energy representing their possibility to undergo a chemical reaction. However, among the four best selected ligands, coptisine, with its lowest gap energy, was predicted to have higher probability to undergo a chemical reaction [25,117] (Table 9 and Fig. 3).
Considering all the aspects of the study and comparing the four best selected compounds, it can be concluded that bellidifolin has the highest potentiality to be used as anti-AD drug. From Table 10, it is clear that bellidifolin generated the most sound and satisfactory performance among the four compounds, with no poor performance in any of the analysis. Even compared to the positive controls, donepezil and galantamine, bellidifolin generated very good results in almost all the analyses. Therefore, bellidifolin can be considered as the most suitable compound as a drug for AD treatment among the compounds considered in this study. However, further investigation could also be conducted on other compounds as they had also shown convincing results in the docking analysis (Table 3). And the authors recommend further in vivo and in vitro experiments to strengthen the findings of this study.

Conclusion
As Alzheimer's disease (AD) is one of the most rapidly growing dementia type diseases in the world, more and more research studies are being conducted to find out the best possible treatment for this disease. Although several treatments are already available, none of these treatments possesses significant satisfactory results. Plants have been known to contain many beneficial agents that can be used to treat a variety of diseases. In this experiment, 50 anti-AChE phytochemicals were selected to analyze their activities against the acetylcholinesterase (AChE) enzyme, which is a key enzyme responsible for the development of AD, by exploiting numerous approaches used in computer-aided drug design. Upon continuous computational experimentation, bellidifolin, naringenin, apigenin, and coptisine were predicted to be the best inhibitors of AChE. Then, their drug potentiality was checked in different post-screening studies where they were also predicted to show quite similar and sound performance, although in some aspects, their performances were not up to the mark. However, more wet-lab based studies must be performed on these 4 agents as well as the other remaining agents to finally confirm their potentiality, safety, and efficacy. Table 10 Overall comparison of the main features of bellidifolin, naringenin, apigenin, and coptisine as potential drugs to treat AD according to the study along with the 2 positive controls from comparison. All the compounds are ranked (excellent, good, and poor) based on their performance in every aspect