Acetylcholinesterase Inhibitory Activities of Essential Oils from Vietnamese Traditional Medicinal Plants

Essential oils are promising as environmentally friendly and safe sources of pesticides for human use. Furthermore, they are also of interest as aromatherapeutic agents in the treatment of Alzheimer’s disease, and inhibition of the enzyme acetylcholinesterase (AChE) has been evaluated as an important mechanism. The essential oils of some species in the genera Callicarpa, Premna, Vitex and Karomia of the family Lamiaceae were evaluated for inhibition of electric eel AChE using the Ellman method. The essential oils of Callicarpa candicans showed promising activity, with IC50 values between 45.67 and 58.38 μg/mL. The essential oils of Callicarpa sinuata, Callicarpa petelotii, Callicarpa nudiflora, Callicarpa erioclona and Vitex ajugifolia showed good activity with IC50 values between 28.71 and 54.69 μg/mL. The essential oils Vitex trifolia subsp. trifolia and Callicarpa rubella showed modest activity, with IC50 values of 81.34 and 89.38, respectively. trans-Carveol showed an IC50 value of 102.88 µg/mL. Molecular docking and molecular dynamics simulation were performed on the major components of the studied essential oils to investigate the possible mechanisms of action of potential inhibitors. The results obtained suggest that these essential oils may be used to control mosquito vectors that transmit pathogenic viruses or to support the treatment of Alzheimer’s disease.


Introduction
Acetylcholinesterase (AChE) is a cholinergic enzyme primarily found at postsynaptic neuromuscular junctions, especially in muscles and nerves. It immediately breaks down or hydrolyzes acetylcholine (ACh), a naturally occurring neurotransmitter, into acetic acid and choline. The primary role of AChE is to terminate neuronal transmission and signaling between synapses to prevent ACh dispersal and activation of nearby receptors [1]. Hydrolysis of acetylcholine is required to allow a cholinergic neuron to return to its resting state after activation [2].
When AChE is inactivated, e.g., by an organophosphorus or carbamate ester, the enzyme is no longer able to hydrolyze ACh; the concentration of ACh in the junction Table 1. Major components of the leaf essential oils of Callicarpa, Premna, Vitex and Karomia species from Vietnam.

Acetylcholinesterase Inhibitory Activity
Acetylcholinesterase inhibitory activity of essential oils of Callicarpa, Premna and Karomia species from Vietnam and their major components have been presented in Tables  2 and 3. The essential oils C. candicans (Da Nang), C. candicans (Quang Nam), C. erioclona, C. sinuata, C. petelotii, C. nudiflora, and V. ajugifolia showed good activity, with IC50 values between 28.71 ± 3.85 and 54.69 ± 3.05 μg/mL. The essential oils of C. candicans (Nghe An)   To the best of our knowledge, this is the first time trans-carveol has been evaluated for AChE inhibitory activity, and it exhibited an IC 50 value of 676 ± 52 µM. Studying the AChE inhibition of monoterpenoids with a p-menthane skeleton observed that monoterpene ketones showed stronger inhibition than the alcohols, and the presence of an isopropenyl group improved the inhibitory strength [30]. Limonene, in the study by Miyazawa and co-authors, showed inhibition at a concentration of 1.2 mM in the range of 22.0-25.0%; however, our study recorded an IC 50 value of 390.2 ± 30.0 µM. This difference may have been due to the experimental method and origin of the AChE used.
Several studies have shown that essential oils with high concentrations of sesquiterpene derivatives exhibit moderate and weak AChE inhibitory activity. A study by Karakaya and co-authors clearly showed that an increase in the concentration of sesquiterpenoids reduced the AChE inhibition of the essential oil, at a concentration of 200 µg/mL, the essential oil from the aerial parts of Salvia verticillata subsp. amasiaca with 60.1% sesquiterpenoids inhibited 51.65 ± 2.05% while its floral essential oil with 78.9% sesquiterpenoids showed inhibition of 42.19 ± 1.55% [33]. Other studies support this trend, such as Salinas et al. 2020 [37], Ali et al. 2012 [38], and Siebert et al. 2015 [39]. Some of the results in this study that were consistent with this trend were those for C. longifolia, C. macrophylla, P. cambodiana, V. pinnata, K. fragrans, which had IC 50 values between 144.3 and 221.85 µg/mL.
Some of the essential oils in this study, although characterized by absolute predominance of sesquiterpenoids, exhibited good and potent AChE inhibitory activity (see above). Similarly, the leaf essential oil of Annona cherimola, characterized by 73.87% of sesquiterpenoids including germacrene D (28.77%), bicyclogermacrene (11.12%), E-β-caryophyllene (10.52%), sabinene (9.05%), and β-pinene (7.93%), showed strong AChE inhibition with an IC 50 value of 41.51 ± 1.02 µg/mL [40]. The essential oil of Eugenia riedeliana was characterized by an absolute predominance of sesquiterpenoids (94.2%), but also exhibited strong inhibition with an IC 50 value of 67.3 µg/mL [41]. When increasing the content of sesquiterpenoids compared with monoterpenoids concentration, an increase in the ability of the essential oil to inhibit AChE was observed [42]. Studies by other groups such as da Silva Barbosa et al. [43] and Miyazawa et al. [44] have supported this trend. Research has shown that the main compound did not contribute to the activity of the essential oil containing it, and in such cases minor components were responsible [45]. Essential oils are a complex mixture of many compounds, most of which have not been studied for their AChE inhibitory activity or studies for their synergistic or antagonistic effects.
Furanosesquiterpenoids are a specific class of compound that has been reported to have AChE inhibitory activity in an in vivo model [46]. Commiterpenes A-C have been reported to have neuroprotective effects [47]. The leaf essential oil of Eugenia uniflora contains atractylone and 3-furanoeudesmene. Both the essential oil and a mixture of the two components showed antinociceptive activity [48]. The furan ring has been suggested to be involved in the inhibitory effect of AChE in previous studies [49]. The essential oil of C. candicans, rich in atractylone, showed excellent larvicidal activity against Aedes aegypti (LC 50 = 2.7-5.34 µg/mL at 24 h) and Culex quinquefasciatus (LC 50 = 1.20-2.04 µg/mL at 24 h) [27].
Callicarpa nudiflora essential oil is characterized by 62.1% monoterpenoid compounds, and this may be one of the reasons responsible for its potent AChE inhibitory activity.

Homology Modeling Study of AChE1 Enzyme
The physicochemical properties of AChE1 analyzed using the Protparam webserver are presented in Table 4. On the other hand, various data on the secondary structures of AChE1, including alpha helix, extended strands and random coil, were predicted using SOPMA (Table 5). According to the data analyzed in Table 4, the AChE1 enzyme contains 91 amino acids, and its molecular weight is 10,716.71 Da. The theoretical pI value was predicted to be 4.76, which suggests that this enzyme has a negative charge. According to Gupta et al. [54], the aliphatic index shows the relative volume occupied by aliphatic residues such as lysine, valine, leucine, and isoleucine. The low value predicted for this parameter (44.95) means that the enzyme cannot be stable in a wide range of temperature. The obtained results presented in Table 5 indicate that the random coil and alpha helix are the main components of the AChE1 enzyme (45.05% and 36.26%, respectively). The constructed model did not possess a beta sheet or turns in the secondary structure.
The modeled structure of AChE1 with a potential binding pocket was predicted using CASTp tool (Figure 2A). A model can be considered to be good when the QMEAN4 score varies between 0 and 1 [55]. In this study, the QMEAN4 score of AChE1 was 0.45, thus, it might be argued that the predicted model is reliable for the performance of further docking studies. In addition, a more negative value of QMEAN4 Z-score indicates the low quality of the constructed model compared with the template structure. The QMEAN4 Z-score of AChE1 model was recorded to have a value of −1.06, suggesting that its quality is significantly high compared to the template structure 6ARX_A ( Figure 2B). The predicted structure of AChE1 was also superimposed on the template model using the TM-align tool [56], and the obtained RMSD value was 0.11 Å. The homology modeling structure was then validated using the SAVES webserver (https://saves.mbi.ucla.edu (accessed on: 20 July 2022). The obtained Ramachandran plot shows that 94.8% of the residues were located in the most favorable region, and only 5.2% of the residues were located in the allowed region ( Figure 2C). No residues were found in the outlier region. The ERRAT program gives an overall quality factor of 100, suggesting that this model is highly reliable ( Figure 2D). The QMEAN4 Z-score of AChE1 model was recorded to have a value of −1.06, suggesting that its quality is significantly high compared to the template structure 6ARX_A ( Figure  2B). The predicted structure of AChE1 was also superimposed on the template model using the TM-align tool [56], and the obtained RMSD value was 0.11 Å. The homology modeling structure was then validated using the SAVES webserver (https://saves.mbi.ucla.edu (accessed on: 20 July 2022). The obtained Ramachandran plot shows that 94.8% of the residues were located in the most favorable region, and only 5.2% of the residues were located in the allowed region ( Figure 2C). No residues were found in the outlier region. The ERRAT program gives an overall quality factor of 100, suggesting that this model is highly reliable ( Figure 2D). For further analysis, EeAChE and AChE1 enzyme models were superimposed to validate the accuracy of the homology modeling process. The structure validation between the EeAChE and AChE1 enzyme models was executed using Chimera 1.13.1. Obtained results indicate these models matched with high identity ( Figure 3); therefore, the For further analysis, EeAChE and AChE1 enzyme models were superimposed to validate the accuracy of the homology modeling process. The structure validation between the EeAChE and AChE1 enzyme models was executed using Chimera 1.13.1. Obtained results indicate these models matched with high identity ( Figure 3); therefore, the EeAChE model will be chosen as a representative structure for docking studies in the latter stage.
EeAChE model will be chosen as a representative structure for docking s latter stage.

Molecular Docking Studies
The best docking conformation of galantamine inside human AChE (P was superimposed with the native ligand. The Root Mean Square Devia value obtained for the superimposition was 0.714 Å (Figure 4).

Molecular Docking Studies
The best docking conformation of galantamine inside human AChE (PDB ID: 4EY6) was superimposed with the native ligand. The Root Mean Square Deviation (RMSD) value obtained for the superimposition was 0.714 Å (Figure 4).
x FOR PEER REVIEW 9 of 22 EeAChE model will be chosen as a representative structure for docking studies in the latter stage.

Molecular Docking Studies
The best docking conformation of galantamine inside human AChE (PDB ID: 4EY6) was superimposed with the native ligand. The Root Mean Square Deviation (RMSD) value obtained for the superimposition was 0.714 Å (Figure 4). The binding free energies, ligand efficacy, and residues participating in the interaction of the studied compounds are tabulated in Table 6.  The binding free energies, ligand efficacy, and residues participating in the interaction of the studied compounds are tabulated in Table 6. Lowest-energy docked poses of the studied compounds suggested by docking simulation are presented in Figure 5. Lowest-energy docked poses of the studied compounds suggested by docking simulation are presented in Figure 5.   Autodock4 is a common docking tool with approximately 6000 citations since 2009 [57]. This program is an open-source package that can predict the binding affinity and dock pose of ligands toward a specific protein target [58]. It has been reported by Gohlke et al. that ligand partial charge calculated using the PM6 method has been shown to greatly increase the docking accuracy and cluster population of the most accurate docking [59]. In this study, the main components of the studied essential oils were studied using Autodock4 to allow a deeper insight into their mechanism of inhibition against the EeAChE and AChE1 enzymes. Initially, galantamine was redocked towards the human AChE enzyme to validate the docking procedure. According to Gowthaman et al., when the RMSD of the dock pose of the co-crystallized ligand is less than 2.0 Å in relation to the native crystallographic pose, docking validation can be considered to be satisfactory [60]. Retrieving the dock pose of co-crystallized ligands, it was possible to validate the docking protocol.
According to the ranking criteria of Autodock4, the more negative the value of the docking score, the higher the binding affinity of the compound towards the targeted receptor [61,62]. The obtained dock score of galantamine was −12.76 kcal/mol; thus, any ligands whose docking energy was close to this threshold would be assumed to exhibit good binding affinity toward the targeted enzyme. As indicated in Table 6, compounds with low docking scores might be assumed to be potential EeAChE inhibitors; these results show good agreement with the inhibition activities obtained from the enzymatic studies.
Binding orientation analysis of galantamine indicated that Glu202 and Ser203 are the key residues participating in forming H-bonds with the reference ligand. This interaction was further stabilized by the hydrophobic bond with Trp86, Gly121, Gly122, Tyr124, Phe295, Phe297, Tyr337, Phe338, His447 ( Figure 5). It should be noted that Trp86, Phe297, Tyr337, Tyr341 and His447 have been reported to participate in constituting the active site of the EeAChE enzyme [63,64].
Of all the docked compounds, limonene exhibited the most negative value of binding free energy (−9.78 kcal/mol) towards EeAChE, suggesting that it binds to this enzyme with the highest binding affinity. Dock pose analysis with the targeted enzyme revealed that Phe297, Phe338 and Tyr341 were the key residues contributing significantly to achieving good docking scores.
The remaining components were ranked in the following order: β-pinene, (E)-β-caryophyllene, α-humulene, trans-carveol and caryophyllene oxide, with Auto-dock4 docking scores of −9.24; −8.97; −6.95; −6.35 and −5.00 kcal/mol, respectively. An array of hydrophobic interactions was observed, which were contributed by Trp86, Tyr337, His447 with β-pinene. Binding mode analysis of (E)-β-caryophyllene showed this compound shared a common non-polar interaction with essential residue Phe297 in Autodock4 is a common docking tool with approximately 6000 citations since 2009 [57]. This program is an open-source package that can predict the binding affinity and dock pose of ligands toward a specific protein target [58]. It has been reported by Gohlke et al. that ligand partial charge calculated using the PM6 method has been shown to greatly increase the docking accuracy and cluster population of the most accurate docking [59]. In this study, the main components of the studied essential oils were studied using Autodock4 to allow a deeper insight into their mechanism of inhibition against the EeAChE and AChE1 enzymes. Initially, galantamine was redocked towards the human AChE enzyme to validate the docking procedure. According to Gowthaman et al., when the RMSD of the dock pose of the co-crystallized ligand is less than 2.0 Å in relation to the native crystallographic pose, docking validation can be considered to be satisfactory [60]. Retrieving the dock pose of co-crystallized ligands, it was possible to validate the docking protocol.
According to the ranking criteria of Autodock4, the more negative the value of the docking score, the higher the binding affinity of the compound towards the targeted receptor [61,62]. The obtained dock score of galantamine was −12.76 kcal/mol; thus, any ligands whose docking energy was close to this threshold would be assumed to exhibit good binding affinity toward the targeted enzyme. As indicated in Table 6, compounds with low docking scores might be assumed to be potential EeAChE inhibitors; these results show good agreement with the inhibition activities obtained from the enzymatic studies.
Binding orientation analysis of galantamine indicated that Glu202 and Ser203 are the key residues participating in forming H-bonds with the reference ligand. This interaction was further stabilized by the hydrophobic bond with Trp86, Gly121, Gly122, Tyr124, Phe295, Phe297, Tyr337, Phe338, His447 ( Figure 5). It should be noted that Trp86, Phe297, Tyr337, Tyr341 and His447 have been reported to participate in constituting the active site of the EeAChE enzyme [63,64].
Of all the docked compounds, limonene exhibited the most negative value of binding free energy (−9.78 kcal/mol) towards EeAChE, suggesting that it binds to this enzyme with the highest binding affinity. Dock pose analysis with the targeted enzyme revealed that Phe297, Phe338 and Tyr341 were the key residues contributing significantly to achieving good docking scores.
The remaining components were ranked in the following order: β-pinene, (E)-βcaryophyllene, α-humulene, trans-carveol and caryophyllene oxide, with Autodock4 docking scores of −9.24; −8.97; −6.95; −6.35 and −5.00 kcal/mol, respectively. An array of hydrophobic interactions was observed, which were contributed by Trp86, Tyr337, His447 with β-pinene. Binding mode analysis of (E)-β-caryophyllene showed this compound shared a common non-polar interaction with essential residue Phe297 in comparison to galantamine. The docked pose was strengthened by additional binding toward Trp286. Regarding the remaining compounds, including α-humulene, trans-carveol, and caryophyllene oxide, despite of their docking conformation with important residues within the binding site of targeted enzyme, the high value of docking scores suggests they have lower potential to be considered as inhibitors for EeAChE.

FPL Simulation
Although the docking protocol often produces suitable results when compared with the experimental results, this method does not consider the receptor dynamics and limiting the number trial position of ligands, which might ultimately result in inaccurate predictions. Thus, a more accurate and precise method could be employed to refine the docking observation. Among various available techniques, the fast pulling of ligand (FPL) technique is a very efficient method with low required CPU time consumption that is able to provide results with high accuracy and precision. In particular, a ligand is forced to travel from bound to unbound states via a harmonic external force ( Figure 6). The physical details during the unbinding process reveal the binding affinity and mechanism of a ligand to the AChE enzyme.

FPL Simulation
Although the docking protocol often produces suitable result the experimental results, this method does not consider the recep iting the number trial position of ligands, which might ultimate predictions. Thus, a more accurate and precise method could be e docking observation. Among various available techniques, the (FPL) technique is a very efficient method with low required CPU is able to provide results with high accuracy and precision. In forced to travel from bound to unbound states via a harmonic ex The physical details during the unbinding process reveal the mechanism of a ligand to the AChE enzyme. The relative binding affinity of a ligand to the AChE enzym eight independent FPL calculations. The pulling force was recor other metrics were monitored every 10 ps. All of the computed over eight independent trajectories. The recorded pulling force a Figure 7. The relative binding affinity of a ligand to the AChE enzyme was estimated using eight independent FPL calculations. The pulling force was recorded every 0.1 ps, and other metrics were monitored every 10 ps. All of the computed values were averaged over eight independent trajectories. The recorded pulling force and work are shown in Figure 7.
The maximum pulling force (F max ), called the rupture force, and the recorded pulling work (W) were used as a criterion for evaluating the ligand affinity. As mentioned in the previous work [65], the pulling work is more appropriate than the rupture force, as it is directly associated with ligand-binding free energy via the isobaric-isothermal Jarzynski equality. The data obtained showed that the pulling forces continuously increased to maximum values before rapidly dropping to zero after the nonbonded contacts between the ligand and the protein were terminated after 400 ps to 700 ps.

Drug-Likeness Studies
As a follow-up to the obtained docking results, the ADMET properties of the "hit" compounds were then analyzed, including human intestinal absorption (HIA), blood-brain barrier (BBB), carcinogenicity, and tumorigenesis potential (Table 7). The maximum pulling force (Fmax), called the rupture force, and the recorded pulling work (W) were used as a criterion for evaluating the ligand affinity. As mentioned in the previous work [65], the pulling work is more appropriate than the rupture force, as it is directly associated with ligand-binding free energy via the isobaric-isothermal Jarzynski equality. The data obtained showed that the pulling forces continuously increased to maximum values before rapidly dropping to zero after the nonbonded contacts between the ligand and the protein were terminated after 400 ps to 700 ps.  It is estimated that a compound with logPS > −3 can be considered to penetrate the central nervous system [66]. The obtained data indicate that all studied compounds satisfy the basic criteria to be a potential pesticide or drug, with HIA values ranging from 94.09% to 96.46%, and the CNS values were determined to be within the range −1.847 to −2.637. It should be noted that caryophyllene oxide was also predicted to have tumorigenic effect, suggesting further caution be used when applied in treatment.

Conclusions
Essential oils from traditional Vietnamese medicinal plants were studied for their AChE inhibitory activity, with most of them showing good inhibitory activity. Species such as C. candicans, C. erioclona, C. rubella, C. nudiflora, P. corymbosa, V. trifolia subsp. litoralis and V. trifolia subsp. trifolia are widely distributed, and their essential oils may be considered to be renewable raw materials. Most of the essential oils in this study were characterized by an absolute predominance of sesquiterpenoids, and C. erioclona, C. sinuata and C. petelotii essential oils exhibited the strongest AChE inhibition. For all of the essential oil components studied, reasonably low docking scores were obtained, which is in good agreement with enzymatic studies showing them to be potential EeAChE and AChE1 enzyme inhibitors. Current knowledge is insufficient to explain these cases, however, and further studies are needed on AChE inhibition by purified monoterpenoids and sesquiterpenoids and their synergistic and antagonistic effects.

Source of Essential Oil
Details of the isolation methods of essential oils are available in our previous articles [27][28][29]. The extraction yield of the essential oils is presented in Table 8.

Acetylcholinesterase (AChE) Inhibition Assay
Acetylcholinesterase (AChE) inhibitory activity of essential oil was determined according to the method described by Ellman and our previous study [67,68]. The stock solution was obtained by dissolving the essential oil in DMSO (Merck), which was then diluted with H 2 O (deionized distilled water) to obtain different experimental concentrations. Each solution mixture consisted of 140 µL of phosphate buffer solution (pH: 8), 20 µL of essential oil at concentrations of 500, 100, 20, and 4 µg/mL, and 20 µL of the enzyme AChE (0.25 IU/mL, Sigma-Aldrich, St. Louis, MO, USA). The reaction mixtures were transferred to the test wells of a 96-well microtiter plate and incubated at 25 • C for 15 min. Then, solutions of 10 µL dithiobisnitrobenzoic acid (DTNB, 2.5 mM, Sigma-Aldrich) and 10 µL acetylthiocholine iodide (ACTI, 2.5 mM, Sigma-Aldrich) were added to each of the test wells and incubation continued for 10 min at 25 • C. At the end of the experiment, the absorbance of each solution was measured at 405 nm. Galantamine was used as a positive control. The negative control well did not contain the test sample. Each test was carried out in triplicate.

Molecular Docking Studies
The major components in the essential oil of studied species were selected for the docking study. The three-dimensional structures of the studied compounds were prepared using MarvinSketch 19.27.0 and PyMOL version 1.3r1 [69]. Energy minimization of studied ligands were conducted using MM2 force field and quantum chemical calculations were performed using the PM6 semiempirical method implemented in Gaussian 09 [70]. Galantamine, a tertiary alkaloid acetylcholinesterase inhibitor, was used as the reference ligand.
The X-ray crystal structure of EeAChE from Electrophorus electricus was retrieved from the Protein Data Bank (RCSB) with PDB ID: 1C2B [71]. The tertiary structure of acetylcholinesterase 1 (AChE1) from Aedes aegypti has not been well determined in previous studies; thus, the structure was constructed by comparative modeling using the SWISS-MODEL webserver (http://swissmodel.expasy.org (accessed on: 20 July 2022). The PDB entry 6ARX_A was selected as template structure for modeling AChE1 enzyme. The amino acid sequence of AChE1 enzyme was already determined and its information is publicity at UniProtKB, archieved under entry ID: UniProtKB-Q8MYC0. The Protparam webserver (https://web.expasy.org/protparam (accessed on: 20 July 2022) was used to calculate the physicochemical properties of the enzyme. To predict the secondary structure of the AChE1 enzyme, the SOPMA tool was used (https://npsa-prabi.ibcp.fr/cgi-bin/npsa_ automat.pl?page=/NPSA/npsa_sopma.html (accessed on: 20 July 2022). The predicted 3D structure was validated using PROCHECK to evaluate backbone conformation based on Psi/Phi Ramachandran plot analysis. AutoDockTools (The Scripps Research Institute, CA, USA) [72] was employed to set up and performed docking calculation. To turn the protein molecule into a free receptor, the heteroatoms including water molecules were deleted. Then, the Kollman charges and solvation parameters were assigned, and Gasteiger charges were added to each atom.
The molecular docking study uses AutoDock4 with Lamarckian genetic algorithm (LGA) to search for the optimum dock pose together with the scoring function to calculate the binding affinity. The binding sites of EeAChE and AChE1 were enclosed in a box with a number of grid points in the x × y × z directions of 50 × 50 × 50, and a grid spacing of 0.375 Å. Initially, AutoGrid was run to generate the grid map of various atoms of the ligands and receptor. After the completion of the grid map, AutoDock was run by using autodock parameters as follows: GA population size, 300; maximum number of energy evaluations, 25,000,000; and number of generations, 27,000. A maximum of 50 conformers were considered for each molecule, and the root-mean-square (RMS) cluster tolerance was set to 2.0 Å in each run. Ligand conformation with the lowest free energy of binding, chosen from the most favored cluster, was selected for the further analysis. The outputs of AutoDock modeling studies were analyzed using PyMOL version 1.3r1 [69] and Discovery Studio Visualizer (San Diego, CA, USA) [73].
As no experimental structure of EeAChE from Electrophorus electricus complexed with an inhibitor is available on the RCSB archive, galantamine was redocked over itself in the crystal structure of human AChE complexed with galantamine (PDB ID: 4EY6) [74] to validate the accuracy of the docking protocol.

Binding Affinity Calculation
To follow up the Autodock 4.2 docking, the Molecular Mechanics/Generalized Born and Surface Area (MM-GBSA) method was used to predict binding free energy using Schrödinger software 2020-2 [74]. This energy is determined by the difference between the complex and the specific energy of the protein and ligand (Equation (1)). This total energy has its heat of reaction calculated by the summary of the solvation free energy (∆Gsol), the gas-phase interaction energy (∆Egas) and the entropy terms (∆S) (Equation (2)). At constant pressure, the heat of the reaction is approximately equal to the change of internal energy, and thus ∆Eint is neglected (Equation (3)). The remaining form of energy includes all intermolecular interactions, for example electrostatic interactions, protein-ligand vdW interactions, ligand desolation, internal strain energies with OPLS2005 force field. The solvation free energy (∆Gsol) consists of non-polar (∆Gsurf) and polar (∆GGB) energy forms and their solvation energy is calculated by sum of the solvent-accessible surface area and GB model (Equation (4)). In general, the binding free energy (∆Gbind) is determined by the sum of solvation free energy and gas-phase interaction energy. At the same time, the entropy term is neglected in the calculation for relative free binding energies [75,76].
∆G bind = ∆H − T∆S ≈ ∆E gas + ∆G sol − T∆S (2) ∆G sol = ∆G GB + ∆G surf (4) In summary, the MM-GBSA calculations are not in agreement with experimental binding affinities, but they show the tendency to bind and the reasonable correlation with the experiment values, and the more negative value of MM-GBSA indicates more potent approximate free energies of binding.

Fast Pulling Ligand Simulation
The structure of the complex was obtained via the molecular docking method. Caver 2.1 [77] was employed to evaluate the disassociate direction of a ligand as suggested by previous works [78]. The complex was then aligned to be the disassociate pathway oriented to the Z-axis. The complex was inserted into a periodic boundary condition rectangular box (7.26 × 9.23 × 12.35 nm) consisting of an enzyme AChE, a ligand, 22,761 water molecules, and 7 counterbalanced ions (Na + ). In particular, the protein was parameterized via the AMBER99SB-ILDN force field [79], and the ligand was represented using general Amber force field [80]. The water molecule was topologized using TIP3P water model [81].
GROMACS version 2016 was employed to carry out the molecular dynamics (MD) simulations. The simulation was performed according to the following four steps: energic minimization, NVT, NPT, and steered MD (SMD) simulations. Specifically, the non-covalent pair was affected within a range of 0.9 nm, and the pair list was updated every 5 fs. The particle mesh Ewald method [82] was employed to mimic the electrostatic interaction with cut-off of 0.9 nm. The van der Waals interaction was affected in a range of 0.9 nm. Both of NVT and NPT simulations were carried out with a length of 100 ps at 300 K. During simulations, the AChE C α atoms were restrained by using a slight harmonic force. The last snapshot of NPT simulations was used as initial structure of SMD simulations. In this scheme, a harmonic external force with a cantilever spring constant of k = 600 kJ/mol/nm 2 and pulling speed v = 0.005 nm/ps was put on the ligand center of mass along the Zdirection. The ligand was disassociated from the binding site of AChE enzyme over the SMD simulations. The data were recorded every 0.1 ps during SMD simulations. The calculations were repeated independently eight times with the same initial conformation to guarantee sufficient sampling.

Data Analysis
Inhibitory data were analyzed by log-probit analysis [83] to acquire IC 50 value as well as 95% confidence limits using Minitab ® version 19.2020.1 (Minitab, LLC, State College, PA, USA).