Galantamine Based Novel Acetylcholinesterase Enzyme Inhibitors: A Molecular Modeling Design Approach

Citation: Silva, L.B.; Ferreira, E.F.B.; Maryam; Espejo-Román, J.M.; Costa, G.V.; Cruz, J.V.; Kimani, N.M.; Costa, J.S.; Bittencourt, J.A.H.M.; Cruz, J.N.; et al. Galantamine Based Novel Acetylcholinesterase Enzyme Inhibitors: A Molecular Modeling Design Approach. Molecules 2023,28, 1035. https://doi.org/10.3390/ molecules28031035 Academic Editor: Dariusz Matosiuk Received: 1 October 2022 Revised: 31 December 2022 Accepted: 9 January 2023 Published: 19 January 2023 Copyright: © 2023 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/). molecules Article Galantamine Based Novel Acetylcholinesterase Enzyme Inhibitors: A Molecular Modeling Design Approach Luciane B. Silva 1,2, Elenilze F. B. Ferreira 3, Maryam 4, JoséM. Espejo-Román5, Glauber V. Costa 2,6, Josiane V. Cruz 2, Njogu M. Kimani 7, Josivan S. Costa 2, JoséA. H. M. Bittencourt 2, Jorddy N. Cruz 2, Joaquín M. Campos 5and Cleydson B. R. Santos 1,2,* 1Graduate Program in Medicinal Chemistry and Molecular Modeling, Health Science Institute, Federal University of Pará, Belém 66075-110, Brazil 2Laboratory of Modeling and Computational Chemistry, Department of Biological and Health Sciences, Federal University of Amapá, Macapá68902-280, Brazil 3Laboratory of Organic Chemistry and Biochemistry, University of the State of Amapá, Macapá68900-070, Brazil 4COMSATS University Islamabad, Park Road, Tarlai Kalan, Islamabad 45550, Pakistan 5Department of Pharmaceutical Organic Chemistry, Faculty of Pharmacy, Campus of Cartuja, University of Granada, 18071 Granada, Spain 6Laboratory of Biotechnology in Natural Products, Department of Biological and Health Sciences, Federal University of Amapá, Macapá68902-280, Brazil 7Department of Physical Sciences, University of Embu, Embu P.O. Box 6-60100, Kenya *Correspondence: br[email protected]; Tel.: +55-(96)-4009-2699 Abstract: Acetylcholinesterase (AChE) enzymes play an essential role in the development of Alzheimer’s disease (AD). Its excessive activity causes several neuronal problems, particularly psychopathies and neuronal cell death. A bioactive pose on the hAChE B site of the human acetylcholinesterase (hAChE) enzyme employed in this investigation, which was obtained from the Protein Data Bank (PDB ID 4EY6), allowed for the prediction of the binding affinity and free binding energy between the protein and the ligand. Virtual screening was performed to obtain structures similar to Galantamine (GNT) with potential hAChE activity. The top 200 hit compounds were prioritized through the use of filters in ZincPharmer, with special features related to the pharmacophore. Critical analyses were carried out, such as hierarchical clustering analysis (HCA), ADME/Tox predictions, molecular docking, molecular simulation studies, synthetic accessibility (SA), lipophilicity, water solubility, and hot spots to confirm the stable binding of the two promising molecules (ZINC16951574-LMQC2, and ZINC08342556-LMQC5). The metabolism prediction, with metabolites M3-2, which is formed by Glutathionation reaction (Phase II), M1-2, and M2-2 formed from the reaction of S-oxidation and Aliphatic hydroxylation (Phase I), were both reactive but with no side effects. Theoretical synthetic routes and prediction of synthetic accessibility for the most promising compounds are also proposed. In conclusion, this study shows that in silico modeling can be used to create new drug candidate inhibitors for hAChE. The compounds ZINC16951574-LMQC2, and ZINC08342556-LMQC5 are particularly promising for oral administration because they have a favorable drug-likeness profile, excellent lipid solubility, high bioavailability, and adequate pharmacokinetics. Keywords: Alzheimer’s disease; ADME; molecular docking and molecular dynamics 1. Introduction Alzheimer’s disease (AD), discovered and named by the German scientist Alois Alzheimer (1907), is a progressive neurodegenerative disease. In the world, approximately 47 million people are affected by this disease, which is projected to increase by 62% before 2050 [ 1 ]. As a neurodegenerative, progressive, irreversible, AD develops due to the loss of neurons in the central nervous system and the malfunctioning of nerve cells. Earlyonset AD (EOAD) is generally hereditary autosomal dominant, constituting only 1–2% Molecules 2023,28, 1035. https://doi.org/10.3390/molecules28031035 https://www.mdpi.com/journal/molecules Molecules 2023,28, 1035 2 of 31 of AD, with genes including amyloid precursor protein (APP), presenilin 1 (PSEN1), and presenilin 2 (PSEN2) being considered the main factors. There is precedence in the literature for identifying two different types of mutations in this gene in patients with familial AD [ 2 ]. Several risk factors have been associated with Alzheimer’s disease, such as neuritic plaques or amyloid plaques and trisomy 21, a risk factor for early-onset dementia, in addition to rapid loss of synapses and degeneration of baseline cholinergic neurons [3]. Clinically, the most common diagnosis is indicated by progressive cognitive deficiencies or defined by loss of memory and learning ability, making it impossible to perform routine activities, and an adverse set of neuropsychiatric symptoms, such as apathy, verbal and physical agitation, irritability, anxiety, depression, delusions and hallucinations [4]. Besides, among the mechanisms related to the appearance and evolution of AD, the “cholinergic hypothesis” has emerged as a widely accepted therapeutic means to improve cognitive functions in AD, and different studies suggest that cholinergic inputs are present in this process [ 5 ]. It has been observed that cholinergic activity can influence amyloid processing. In the absence of muscarinic receptor activity, the amyloidogenic pathway is privileged, while the opposite occurs with regular receptor activity [ 6 ]. Another hypothesis is the “amyloid cascade,” which was first proposed in 1992. Other conjectures, such as the oligomeric and metallic hypotheses, can be extensions of the amyloid hypothesis and began taking on more significant proportions in the 1990s [7]. In recent years, with a comprehensive discussion about the amyloid cascade hypothesis, a growing body of evidence has come to suggest that endogenous metal ions, particularly those that have redox activity, such as copper (II) and iron (III), in addition to certain non-redox-active ions, such as zinc (II), can contribute to the evolution of neurodegenerative diseases, favoring the aggregation of Aβand increasing its toxicity [8,9]. AD is also associated with the deficit of the neurotransmitter Acetylcholine (ACh) and oxidative stress caused by an exacerbation of glutamatergic transmission [10]. Acetylcholinesterase (AChE) is the essential enzyme in the serine hydrolases family in cholinergic synapses that plays a crucial role in memory and cognition [ 11 ]. AChE inhibition has also been recognized as a therapeutic strategy for other disorders, such as dementia, myasthenia gravis, glaucoma, and Parkinson’s disease, in addition to AD [12]. Currently, AChE inhibitors, including Donepezil (Aricept ® ), Galantamine (GNT) (Reminyl ® ) (Figure 1), and Rivastigmine (Exelon ® ), are widely used in symptomatic treatments for AD. However, the effectiveness of these drugs is hindered by their side effects, such as gastrointestinal disorders, hepatotoxicity, dizziness, diarrhea, vomiting, nausea, pharmacokinetic disadvantages, and hypotension [ 4 ]. Thus, it is necessary to discover new, more effective compounds to reduce the risk. Molecules 2023, 28, x FOR PEER REVIEW 2 of 31 1. Introduction Alzheimer’s disease (AD), discovered and named by the German scientist Alois Alzheimer (1907), is a progressive neurodegenerative disease. In the world, approximately 47 million people are affected by this disease, which is projected to increase by 62% before 2050 [1]. As a neurodegenerative, progressive, irreversible, AD develops due to the loss of neurons in the central nervous system and the malfunctioning of nerve cells. Early-onset AD (EOAD) is generally hereditary autosomal dominant, constituting only 1–2% of AD, with genes including amyloid precursor protein (APP), presenilin 1 (PSEN1), and presenilin 2 (PSEN2) being considered the main factors. There is precedence in the literature for identifying two different types of mutations in this gene in patients with familial AD [2]. Several risk factors have been associated with Alzheimer’s disease, such as neuritic plaques or amyloid plaques and trisomy 21, a risk factor for early-onset dementia, in addition to rapid loss of synapses and degeneration of baseline cholinergic neurons [3]. Clinically, the most common diagnosis is indicated by progressive cognitive deficiencies or defined by loss of memory and learning ability, making it impossible to perform routine activities, and an adverse set of neuropsychiatric symptoms, such as apathy, verbal and physical agitation, irritability, anxiety, depression, delusions and hallucinations [4]. Besides, among the mechanisms related to the appearance and evolution of AD, the “cholinergic hypothesis” has emerged as a widely accepted therapeutic means to improve cognitive functions in AD, and different studies suggest that cholinergic inputs are present in this process [5]. It has been observed that cholinergic activity can influence amyloid processing. In the absence of muscarinic receptor activity, the amyloidogenic pathway is privileged, while the opposite occurs with regular receptor activity [6]. Another hypothesis is the “amyloid cascade,” which was first proposed in 1992. Other conjectures, such as the oligomeric and metallic hypotheses, can be extensions of the amyloid hypothesis and began taking on more significant proportions in the 1990s [7]. In recent years, with a comprehensive discussion about the amyloid cascade hypothesis, a growing body of evidence has come to suggest that endogenous metal ions, particularly those that have redox activity, such as copper (II) and iron (III), in addition to certain non-redox-active ions, such as zinc (II), can contribute to the evolution of neurodegenerative diseases, favoring the aggregation of Aβ and increasing its toxicity [8,9]. AD is also associated with the deficit of the neurotransmitter Acetylcholine (ACh) and oxidative stress caused by an exacerbation of glutamatergic transmission [10]. Acetylcholinesterase (AChE) is the essential enzyme in the serine hydrolases family in cholinergic synapses that plays a crucial role in memory and cognition [11]. AChE inhibition has also been recognized as a therapeutic strategy for other disorders, such as dementia, myasthenia gravis, glaucoma, and Parkinson’s disease, in addition to AD [12]. Currently, AChE inhibitors, including Donepezil (Aricept®), Galantamine (GNT) (Reminyl®) (Figure 1), and Rivastigmine (Exelon®), are widely used in symptomatic treatments for AD. However, the effectiveness of these drugs is hindered by their side effects, such as gastrointestinal disorders, hepatotoxicity, dizziness, diarrhea, vomiting, nausea, pharmacokinetic disadvantages, and hypotension [4]. Thus, it is necessary to discover new, more effective compounds to reduce the risk. Figure 1. (A) 2D and (B) 3D structure of the GNT compound. Green are carbon atoms, red oxygen atoms and blue nitrogen atoms. Figure 1. ( A ) 2D and ( B ) 3D structure of the GNT compound. Green are carbon atoms, red oxygen atoms and blue nitrogen atoms. Thus, this paper uses molecular modeling approaches to select new compounds with potential hAChE inhibition capacity. For this, we use the GNT as a control to filter new compounds. We sought to obtain new molecules with potential inhibitory activity hAChE via a pharmacophore model generation. The generated pharmacophore model was used to screen chemical databases virtually. The successful compounds were filtered by evaluating their drug-like properties, which were statistically assessed using Pearson’s correlations and hierarchical clustering analysis (HCA). Subsequently, they were submitted Molecules 2023,28, 1035 3 of 31 to pharmacokinetic, toxicological, and biological activity predictions. The successful compounds’ binding mode analyses were performed using docking and molecular dynamics simulations. Once identified, the compounds were subsequently subjected to metabolism prediction, ADME/Tox to metabolites, hot spots, and theoretical synthetic routes proposed for the most promising compounds. The methodological scheme of the steps used in this study is shown in Figure 2. Molecules 2023, 28, x FOR PEER REVIEW 3 of 31 Thus, this paper uses molecular modeling approaches to select new compounds with potential hAChE inhibition capacity. For this, we use the GNT as a control to filter new compounds. We sought to obtain new molecules with potential inhibitory activity hAChE via a pharmacophore model generation. The generated pharmacophore model was used to screen chemical databases virtually. The successful compounds were filtered by evaluating their drug-like properties, which were statistically assessed using Pearson’s correlations and hierarchical clustering analysis (HCA). Subsequently, they were submitted to pharmacokinetic, toxicological, and biological activity predictions. The successful compounds’ binding mode analyses were performed using docking and molecular dynamics simulations. Once identified, the compounds were subsequently subjected to metabolism prediction, ADME/Tox to metabolites, hot spots, and theoretical synthetic routes proposed for the most promising compounds. The methodological scheme of the steps used in this study is shown in Figure 2. Figure 2. Main methodological steps. (1) Database to search for the compounds, (2) (3) and (4) generation and validation of the pharmacophoric model, (5) virtual screening, (6) evaluation of the ADME/Tox in silico properties of the selected compounds, (7) evaluation of the mode of interaction by docking and molecular dynamics and finally (8) proposed compounds. 2. Results and Discussion 2.1. Pharmacophore Model Generation In PharmaGist the 15 selected molecules (Figure S1 in Supplementary Material) were used as the input, and GNT was added as a template, with a score of 32.692. Subsequently, Figure 2. Main methodological steps. (1) Database to search for the compounds, (2) (3) and (4) generation and validation of the pharmacophoric model, (5) virtual screening, (6) evaluation of the ADME/Tox in silico properties of the selected compounds, (7) evaluation of the mode of interaction by docking and molecular dynamics and finally (8) proposed compounds. 2. Results and Discussion 2.1. Pharmacophore Model Generation In PharmaGist the 15 selected molecules (Figure S1 in Supplementary Materials) were used as the input, and GNT was added as a template, with a score of 32.692. Subsequently, a matrix with the following pharmacophoric descriptors: atoms (ATM), spatial characteristics (SF), aromatic (ARO), hydrophobic (HYD), and Hydrogen-bond acceptor (ACC) was constructed, and their values evaluated using the statistical software STAT. The statistical method of Pearson’s correlation aimed to show the correlation between pharmacophoric descriptors and the inhibitory activity of the molecules (pIC50 =−logIC50) so that it was possible to evaluate between the correlations which values should be considered in the analysis as described by Ferreira et al. 2019 [13]. The pharmacophore descriptors and the biological activity described have a proportional correlation. For instance, a Hydrogen-bond acceptor (ACC) with a p-value of − 0.804 Molecules 2023,28, 1035 4 of 31 shows the number of hydrogen acceptors that correlate with the pIC 50 response of the selected molecules. Significantly, as with the descriptor Atoms (ATM) with a value of p=−0.634, each pharmacophoric descriptor contributes to an inhibitory activity (pIC50). Also, a note from Table 1, the correlation between pharmacological drug parameters is less than 0.928. In contrast, the correlation between inhibitory activity (pIC 50 ) is more significant than − 0.607, that is, an inversely proportional correlation, especially to the Hydrogen-bond acceptor (ACC). The selected pharmacophoric characteristics represent the necessary characteristics for generating pharmacophoric models in the search to identify potential compounds with AChE enzyme inhibitory activity. Table 1. Result of the matrix with 6 descriptors and the Pearson correlation of the variables resulting from the virtual screening ligand-based molecules. Structures ATM SF ARO HYD ACC pIC50 1 42 11 1 6 4 8.250 2 85 16 2 9 5 8.060 3 84 17 2 10 5 7.909 4 85 16 2 9 5 7.829 5 81 15 2 8 5 7.750 6 83 16 2 8 6 7.619 7 80 16 2 8 6 7.460 8 98 21 2 13 6 7.260 9 82 16 2 8 6 7.239 10 109 24 2 17 5 7.219 11 103 22 2 15 5 7.050 12 101 22 2 14 6 7.030 13 87 17 2 8 7 7.010 14 84 16 2 7 7 6.949 15 90 18 2 9 7 6.920 16 93 19 2 10 7 6.659 SF 0.928 - - - - - ARO 0.806 0.541 - - - - HYD 0.802 0.944 0.335 - - - ACC 0.338 0.187 0.501 −0.143 - - pIC50 −0.634 −0.607 −0.494 −0.355 −0.804 - 2.2. Evaluation of the Pharmacophore Model From the result of the Pearson correlation, the model can be evaluated by a chemometric study. A dendrogram (Figure 3) was obtained from a hierarchical clustering analysis (HCA) using the Pirouette 4.0 software. Confirmation of the data obtained by Pearson’s correlation was conducted by the generation of HCA pharmacophoric hypotheses in which a correlation of the biological activity (pIC 50 ) as an independent variable and the structural similarity cursor in categories, namely, more active (a) and less active (b), originating from six molecular descriptors, including Atoms (ATM), Spatial characteristics (SF), aromatic (ARO), hydrophobic (HYD), and Hydrogen-bond acceptor (ACC), was performed (Figure 3). The descriptor Hydrogen-bond acceptor (ACC) in the pink cluster shows a greater approximation with the pIC 50 . This confirms the value obtained in the Pearson correlation. The descriptors Hydrophobic (HYD), Spatial Resources (SF), Atom (ATM), and aromatic (ARO) are grouped into the same cluster (in blue). Molecules 2023,28, 1035 5 of 31 Molecules 2023, 28, x FOR PEER REVIEW 5 of 31 (SF), aromatic (ARO), hydrophobic (HYD), and Hydrogen-bond acceptor (ACC), was performed (Figure 3). Figure 3. HCA dendrogram, a correlation between pharmacophoric characteristics and pIC50. The descriptor Hydrogen-bond acceptor (ACC) in the pink cluster shows a greater approximation with the pIC50. This confirms the value obtained in the Pearson correlation. The descriptors Hydrophobic (HYD), Spatial Resources (SF), Atom (ATM), and aromatic (ARO) are grouped into the same cluster (in blue). The HCA technique showed a similar dendrogram in which the molecules were classified into two classes (more active and less active) according to their similarities. In the largest cluster, called the most active (in blue), are molecules 1–12 with the best inhibitory activity values. In the smallest cluster (in pink), molecules 13–16 are the least active (Figure 4). Figure 4. HCA dendrogram for 15 molecules derived from Gal with potential inhibitory activity of hAChE and GNT identified as number 1, the clusters are divided into (a) Most active represented in blue, and (b) less active, represented in pink. Thus, by the dendrogram in Figure 4, it is quickly perceived that molecules 1 and 9 are the most distinct compared to the other most active molecules (compounds 1, 2, and 3). However, they are very close, indicating a high chemical similarity. Figure 3. HCA dendrogram, a correlation between pharmacophoric characteristics and pIC50. The HCA technique showed a similar dendrogram in which the molecules were classified into two classes (more active and less active) according to their similarities. In the largest cluster, called the most active (in blue), are molecules 1–12 with the best inhibitory activity values. In the smallest cluster (in pink), molecules 13–16 are the least active (Figure 4). Molecules 2023, 28, x FOR PEER REVIEW 5 of 31 (SF), aromatic (ARO), hydrophobic (HYD), and Hydrogen-bond acceptor (ACC), was performed (Figure 3). Figure 3. HCA dendrogram, a correlation between pharmacophoric characteristics and pIC50. The descriptor Hydrogen-bond acceptor (ACC) in the pink cluster shows a greater approximation with the pIC50. This confirms the value obtained in the Pearson correlation. The descriptors Hydrophobic (HYD), Spatial Resources (SF), Atom (ATM), and aromatic (ARO) are grouped into the same cluster (in blue). The HCA technique showed a similar dendrogram in which the molecules were classified into two classes (more active and less active) according to their similarities. In the largest cluster, called the most active (in blue), are molecules 1–12 with the best inhibitory activity values. In the smallest cluster (in pink), molecules 13–16 are the least active (Figure 4). Figure 4. HCA dendrogram for 15 molecules derived from Gal with potential inhibitory activity of hAChE and GNT identified as number 1, the clusters are divided into (a) Most active represented in blue, and (b) less active, represented in pink. Thus, by the dendrogram in Figure 4, it is quickly perceived that molecules 1 and 9 are the most distinct compared to the other most active molecules (compounds 1, 2, and 3). However, they are very close, indicating a high chemical similarity. Figure 4. HCA dendrogram for 15 molecules derived from Gal with potential inhibitory activity of hAChE and GNT identified as number 1, the clusters are divided into (a) Most active represented in blue, and (b) less active, represented in pink. Thus, by the dendrogram in Figure 4, it is quickly perceived that molecules 1 and 9 are the most distinct compared to the other most active molecules (compounds 1, 2, and 3). However, they are very close, indicating a high chemical similarity. In addition, the Hydrogen-bond acceptor (ACC), Hydrophobic (HYD), and Atoms (ATM) acceptance number is a better prediction observed by molecule 10 (5 ACC, 17 HYD, and 109 ATM), which shows the highest values, as well as 11, 8, 12. Molecules +2 and +4 that are highly similar and have the same ACC, HYD, and ATM values also have better predictions. It is noticed that for the most active molecules, the descriptors ATM and Molecules 2023,28, 1035 6 of 31 HYDRO resources are relevant to the value of the inhibitory activity, shifting them to the most active side, confirming the importance of the number of Hydrogen-bond acceptors (ACC). However, in molecules +7 and +9 having the same number of ACC and HYD, their inhibitory activity is accentuated by the number of atoms in their structure (80 and 82). It is possible to observe the presence of the indole group in the most active molecules by the similarity that ends the side chain linked to the N atom. The indole fraction of one of the molecules binds in the same region, where A β join hAChE, where they share the same binding site as GNT [ 14 ]. The indole fraction was selected as a suitable parameter for binding to aromatic residues in PAS due to its ability to participate in hydrophobic interactions [15]. 2.3. Pharmacophoric-Based Hierarchical Virtual Screening The pharmacological model obtained through PharmaGist was submitted to the ZINCpharmer server to obtain the spatial coordinates of the pharmacophore. The aligned molecules shared 3 pharmacophoric characteristics: Aromatic (ARO), Hydrophobic (HYD), and Hydrogen-bond acceptor (ACC) (Table 2and Figure 5), generating a model with the following coordinates: Table 2. Pharmacophore characteristics. Pharmacophoric Characteristics Coordinates X Y Z Radius Molecules 2023, 28, x FOR PEER REVIEW 6 of 31 In addition, the Hydrogen-bond acceptor (ACC), Hydrophobic (HYD), and Atoms (ATM) acceptance number is a better prediction observed by molecule 10 (5 ACC, 17 HYD, and 109 ATM), which shows the highest values, as well as 11, 8, 12. Molecules +2 and +4 that are highly similar and have the same ACC, HYD, and ATM values also have better predictions. It is noticed that for the most active molecules, the descriptors ATM and HYDRO resources are relevant to the value of the inhibitory activity, shifting them to the most active side, confirming the importance of the number of Hydrogen-bond acceptors (ACC). However, in molecules +7 and +9 having the same number of ACC and HYD, their inhibitory activity is accentuated by the number of atoms in their structure (80 and 82). It is possible to observe the presence of the indole group in the most active molecules by the similarity that ends the side chain linked to the N atom. The indole fraction of one of the molecules binds in the same region, where Aᵦ join hAChE, where they share the same binding site as GNT [14]. The indole fraction was selected as a suitable parameter for binding to aromatic residues in PAS due to its ability to participate in hydrophobic interactions [15]. 2.3. Pharmacophoric-Based Hierarchical Virtual Screening The pharmacological model obtained through PharmaGist was submitted to the ZINCpharmer server to obtain the spatial coordinates of the pharmacophore. The aligned molecules shared 3 pharmacophoric characteristics: Aromatic (ARO), Hydrophobic (HYD), and Hydrogen-bond acceptor (ACC) (Table 2 and Figure 5), generating a model with the following coordinates: Table 2. Pharmacophore characteristics. Pharmacophoric Characteristics Coordinates X Y Z Radius Aromatic (ARO) 8.20 −59.95 −3.93 1.1 Hydrophobic (HYD) 8.49 −62.48 −25.88 1.0 Hydrogen bond acceptor (ACC) 7.59 −61.43 −22.89 0.5 Figure 5. Drug model selected for virtual screening. Orange sphere (Aromatic Ring), Pink (Hydrophobic), and White (Hydrogen-bond acceptor). (A) The pharmacophoric hypothesis was tested on a set of selected compounds. (B) Hypothesis and pharmacophoric analysis of the active site: the three characteristics are in regions of interactions with seven essential amino acids at the hAChE receptor. As a result of the virtual screening used in medicines, Top Hit 200 compounds were prioritized through filters in ZincPharmer, with special features related to the pharmacophore. Success for the compound shows that hAChE inhibitory activity depends on the Aromatic (ARO) 8.20 − 59.95 −3.93 1.1 Hydrophobic (HYD) 8.49 − 62.48 − 25.88 1.0 Hydrogen bond acceptor (ACC) 7.59 − 61.43 − 22.89 0.5 Molecules 2023, 28, x FOR PEER REVIEW 6 of 31 In addition, the Hydrogen-bond acceptor (ACC), Hydrophobic (HYD), and Atoms (ATM) acceptance number is a better prediction observed by molecule 10 (5 ACC, 17 HYD, and 109 ATM), which shows the highest values, as well as 11, 8, 12. Molecules +2 and +4 that are highly similar and have the same ACC, HYD, and ATM values also have better predictions. It is noticed that for the most active molecules, the descriptors ATM and HYDRO resources are relevant to the value of the inhibitory activity, shifting them to the most active side, confirming the importance of the number of Hydrogen-bond acceptors (ACC). However, in molecules +7 and +9 having the same number of ACC and HYD, their inhibitory activity is accentuated by the number of atoms in their structure (80 and 82). It is possible to observe the presence of the indole group in the most active molecules by the similarity that ends the side chain linked to the N atom. The indole fraction of one of the molecules binds in the same region, where Aᵦ join hAChE, where they share the same binding site as GNT [14]. The indole fraction was selected as a suitable parameter for binding to aromatic residues in PAS due to its ability to participate in hydrophobic interactions [15]. 2.3. Pharmacophoric-Based Hierarchical Virtual Screening The pharmacological model obtained through PharmaGist was submitted to the ZINCpharmer server to obtain the spatial coordinates of the pharmacophore. The aligned molecules shared 3 pharmacophoric characteristics: Aromatic (ARO), Hydrophobic (HYD), and Hydrogen-bond acceptor (ACC) (Table 2 and Figure 5), generating a model with the following coordinates: Table 2. Pharmacophore characteristics. Pharmacophoric Characteristics Coordinates X Y Z Radius Aromatic (ARO) 8.20 −59.95 −3.93 1.1 Hydrophobic (HYD) 8.49 −62.48 −25.88 1.0 Hydrogen bond acceptor (ACC) 7.59 −61.43 −22.89 0.5 Figure 5. Drug model selected for virtual screening. Orange sphere (Aromatic Ring), Pink (Hydrophobic), and White (Hydrogen-bond acceptor). (A) The pharmacophoric hypothesis was tested on a set of selected compounds. (B) Hypothesis and pharmacophoric analysis of the active site: the three characteristics are in regions of interactions with seven essential amino acids at the hAChE receptor. As a result of the virtual screening used in medicines, Top Hit 200 compounds were prioritized through filters in ZincPharmer, with special features related to the pharmacophore. Success for the compound shows that hAChE inhibitory activity depends on the Figure 5. Drug model selected for virtual screening. Orange sphere (Aromatic Ring), Pink (Hydrophobic), and White (Hydrogen-bond acceptor). ( A ) The pharmacophoric hypothesis was tested on a set of selected compounds. ( B ) Hypothesis and pharmacophoric analysis of the active site: the three characteristics are in regions of interactions with seven essential amino acids at the hAChE receptor. As a result of the virtual screening used in medicines, Top Hit 200 compounds were prioritized through filters in ZincPharmer, with special features related to the pharmacophore. Success for the compound shows that hAChE inhibitory activity depends on the precision and specificity of the activated pharmacophore [ 13 ]. ZINC includes 1.4 billion compounds, which continue to grow in size with ZINC20 and 1.3 billion of which are purchasable, sourced from 310 catalogs from 150 companies. Over 90% of catalogs are refreshed every 90 days, and over 90% of compounds have been verified as purchasable within the last three months [ 16 ]. A virtual screening through the library composed of commercial ZINC molecules and internal molecules (real/virtual files extended from the real scaffold) as an internal 3D database prepared for the virtual tracking of any model [ 16 ]. Thus, a molecular adjustment was applied to filtered molecules from the virtual screening based on the pharmacophore model. Molecules 2023,28, 1035 7 of 31 2.4. Pharmacokinetic Predictions Our pharmacokinetic selection was designed to determine a molecule with physicochemical properties equal to or better than the reference molecule, GNT. In Table 3, it can be seen that the predictions of the pharmacokinetic properties of the selected molecules from the virtual screening reduced the number of compounds in a top hit 200 to obtain, in this phase, 8 molecules through the following physical-chemical descriptors. Table 3. Pharmacokinetic predictions for GAL and molecules obtained by virtual screening. Compounds #Star a RO5 b%HOA cQPP Caco dQPP MDCK eQPlog Po/w fCNS gQPlogBB hMW i 0–4 <25 Poor >80 Great <25 Poor >500 Great <25 Poor >500 Great <5 −2 (Inactive) +2 (Active) <1 150–500 GNT (control) 1 0 90.0 716.5000 381.7100 2.042 1 0.368 287.350 ZINC86196920 0 0 100 4677.623 2621.609 2.786 1 0.105 208.257 ZINC16951574 (LMQC2) 0 0 100 1838.444 1056.986 3.174 2 0.611 313.397 ZINC91960073 0 0 100 5925.536 7466.586 2.815 1 0.410 253.685 ZINC08342556 (LMQC5) 0 0 100 2928.309 2764.505 3.828 1 0.144 340.452 ZINC86199797 0 0 100 4673.581 4673.581 3.166 1 0.029 222.284 ZINC13108311 0 0 100 2387.068 5990.536 3.254 1 0.081 350.234 ZINC13362890 0 0 100 3705.456 2037.983 2.415 1 0.111 230.267 ZINC21657754 0 0 100 3792.205 2089.602 2.136 1 0.065 246.266 a Number of computed properties that fall outside the required range for 95% of the known drug; b Number of violations of Lipinski’s ‘Rule of Five.’ c Percentage of human oral absorption (%HOA) (acceptable range: <25% is poor and >80% is high); d Predicted Caco-2 cell permeability in nm/s (acceptable range: <25 is poor and >500 is great); e Predicted apparent MDCK cell permeability in nm/s (acceptable range: <25 is poor and >500 is great); f the apparent permeability of compound between octanol/water (QPlogPo/w) < 5; g CNS Predicted activity of the central nervous system on a scale of − 2 (inactive) to +2 (active). h the apparent permeability of compound in blood–the brain barrier (QPlogBB). iMolar weight. The #Star parameter evaluated the pharmacokinetic profile of the compounds with a general pharmacokinetic acceptance rating for the drug similarity parameter, which indicates property descriptors for value ranges optimized for 95% of the drugs [ 17 ]. In this analysis, all selected compounds showed values equal to zero, i.e., no violation, indicating an important similarity with the commercially available drugs. However, GNT exhibited 1 violation within these parameters (Star). There were no violations of Lipinski’s rule of five, according to the parameter for which all the molecules have values equal to zero [ 18 ]. Additionally, all the molecules have favorable properties concerning the percentage of oral absorption, a desired characteristic for developing drug candidates that can be efficient for oral administration. This data was corroborated by the cell models QPP Caco (nm/s) and QPP MDCK (nm/s), which also refers to the prediction of oral absorption of the drug, used as a model for the intestinal blood barrier. The QPP Caco results show values above 1838.444, more than double that of the reference molecule (GNT), values that estimate excellent intestinal absorption. Since the permeability of the blood-brain barrier cells is stimulated by the MDCK cells (QPP MDCK ), all compounds exhibited good results within parameters considered to be above 500 nm/s. For the results obtained by the parameter QplogBB (blood and brain partition), values less than 1 must be used, and for QPlogPo/w (octanol/water), values must be less than 5; Positive values of the log BB mean that these compounds prefer brain tissue to blood, while a value close to zero means that a compound has a Brain/Blood concentration ratio, respectively. The results of QPlogBB for selected compounds showed values less than 0.611, while QPlogPo/w values were less than 3.828; both are good predictions for the compounds [ 17 ]. MW, corresponding to the molecular weight of the molecule in (g · mol −1 ), must obey the molar mass (MW) ≤ 500 g/mol for good oral bioavailability. The analysis of MW distributions showed a peak value between 208.257 and 350.234 g · mol −1 . All the compounds do not violate this MW parameter, presenting values for good oral availability. Molecules 2023,28, 1035 8 of 31 2.5. Toxicological Predictions Derek software [ 19 , 20 ] classifies potentially toxic compounds into six categories: right, probable, plausible, implausible, improbable, or impossible (Table 4). The compounds were classified in the plausible category. A few selected compounds were similar to the reference molecule (GNT) and showed no toxicity alerts. The molecules that showed the best results in this step are ZINC86196920, ZINC16951574 (LMQC2), ZINC08342556 (LMQC5), ZINC86199797, ZINC13362890, and ZINC21657754. ZINC15910273 molecule generated a toxicophoric alert of skin sensitization due to the presence of hydrazine or precursor in its structures. A structural skin sensitization alert within a molecule indicates which molecule has the potential to cause skin sensitization. Whether a molecule is a skin sensitizer will also depend on its percutaneous absorption, having lipophilic characteristics that are absorbed more quickly by the skin and, therefore, are more likely to cause sensitization [ 21 ]. Therefore, these molecules were excluded due to their undesired potential toxic activities. Table 4. Toxicity prediction by toxicophoric identification of acute oral and hepatotoxicity. Compounds Toxicity Prediction Alert (Lhasa Prediction) Toxicophoric Group LD50 Predicted in Rodents (mg/kg) Toxicity Class * Hepatotoxicity GNT (control) - - 85 III Inactive ZINC86196920 - - 362 IV Inactive ZINC16951574 (LMQC2) - - 520 IV Inactive ZINC08342556 (LMQC5) - - 1000 IV Inactive ZINC13362890 - - 360 IV Inactive ZINC86199797 - - 362 IV Inactive ZINC21657754 - - 800 IV Active ZINC13108311 Skin sensitization Thiol or thiol exchange agent Molecules 2023, 28, x FOR PEER REVIEW 9 of 31 ZINC21657754 - - 800 IV Active ZINC13108311 Skin sensitization Thiol or thiol exchange agent 1000 IV Active ZINC15910273 Hydrazine or precursor 1000 IV Active * If swallowed = Class I: fatal (LD50 ≤ 5). Class II: fatal (5 < LD50 ≤ 50). Class III: toxic (50 < LD50 ≤ 300). Class IV: harmful (300 < LD50 ≤ 2000). Class V: may be harmful (2000 < LD50 ≤ 5000), Class VI: nontoxic (LD50 > 5000). 2.6. Prediction of Biological Activity (PASS) The prediction of potential biological activity was performed using PASS [23] resulting in six selected compounds. Table 5 shows the biological activities that were considered: a CYP2D6 substrate, a CYP3A4 substrate, treatments for Alzheimer’s disease, and a CYP3A4 inducer. The values of Pa and Pi can vary from 0.000 to 1000, recorded as the compound’s probability of being inactive or active, respectively [24]. The biological activities obtained from the PASS were related to aspects of similarity with other bioactive substances. Table 5. Prediction of biological activity of compounds used for virtual screening using the PASS online server. Compound Pa Pi Biological Activity GNT (control) 0.957 0.003 Substrate CYP2D6 0.553 0.052 Substrate CYP3A4 0.428 0.055 Treatment of Alzheimer’s disease ZINC86196920 0.408 0.056 Substrate CYP2D6 0.413 0.085 Substrate CYP3A4 0.223 0.124 Treatment of Alzheimer’s disease ZINC16951574 (LMQC2) 0.553 0.024 Substrate CYP2D6 0.264 0.187 Substrate CYP3A4 ZINC08342556 (LMQC5) - - - ZINC86199797 0.544 0.025 Substrate CYP2D6 0.532 0.055 Substrate CYP3A4 ZINC13362890 0.191 0.165 Treatment of Alzheimer’s disease ZINC21657754 0.192 0.164 Treatment of Alzheimer’s disease Pa = probability of being active (Pa > 0.000 or Pa = 1.000); Pi = Probability of being inactive (Pi = 0.000). Our results revealed that GNT has activity against the substrate CYP2D6 (Pa = 0.957) and CYP3A4 (Pa = 0.553) and treatment of Alzheimer’s disease (Pa = 0.428), confirming the results in silico. AChE has low pharmacokinetics in its interactions. However, since GNT is metabolized in the liver through CYP2D6 and CYP3A4, its hepatic metabolism can be affected by specific substrates, inhibitors, or enhancers of the same enzymes [25]. According to Tables 5 and 6, the biological activity of these compounds revealed that ZINC86196920, ZINC16951574 (LMQC2), and ZINC86199797 have CYP2D6 and CYP3A4 substrate activities, all with Pa ˃ 0.5. This implies they have the greatest chance of being similar to other bioactive compounds. However, compounds with Pa < 0.5 did not show any reports about their activity in the database, that is, they do not have reported activities in the literature; they were also selected for the next stages of the process when their Pa ˃ Pi, as reported by Rodrigues et al. 2013 [26]. 1000 IV Active ZINC15910273 Hydrazine or precursor Molecules 2023, 28, x FOR PEER REVIEW 9 of 31 ZINC21657754 - - 800 IV Active ZINC13108311 Skin sensitization Thiol or thiol exchange agent 1000 IV Active ZINC15910273 Hydrazine or precursor 1000 IV Active * If swallowed = Class I: fatal (LD50 ≤ 5). Class II: fatal (5 < LD50 ≤ 50). Class III: toxic (50 < LD50 ≤ 300). Class IV: harmful (300 < LD50 ≤ 2000). Class V: may be harmful (2000 < LD50 ≤ 5000), Class VI: nontoxic (LD50 > 5000). 2.6. Prediction of Biological Activity (PASS) The prediction of potential biological activity was performed using PASS [23] resulting in six selected compounds. Table 5 shows the biological activities that were considered: a CYP2D6 substrate, a CYP3A4 substrate, treatments for Alzheimer’s disease, and a CYP3A4 inducer. The values of Pa and Pi can vary from 0.000 to 1000, recorded as the compound’s probability of being inactive or active, respectively [24]. The biological activities obtained from the PASS were related to aspects of similarity with other bioactive substances. Table 5. Prediction of biological activity of compounds used for virtual screening using the PASS online server. Compound Pa Pi Biological Activity GNT (control) 0.957 0.003 Substrate CYP2D6 0.553 0.052 Substrate CYP3A4 0.428 0.055 Treatment of Alzheimer’s disease ZINC86196920 0.408 0.056 Substrate CYP2D6 0.413 0.085 Substrate CYP3A4 0.223 0.124 Treatment of Alzheimer’s disease ZINC16951574 (LMQC2) 0.553 0.024 Substrate CYP2D6 0.264 0.187 Substrate CYP3A4 ZINC08342556 (LMQC5) - - - ZINC86199797 0.544 0.025 Substrate CYP2D6 0.532 0.055 Substrate CYP3A4 ZINC13362890 0.191 0.165 Treatment of Alzheimer’s disease ZINC21657754 0.192 0.164 Treatment of Alzheimer’s disease Pa = probability of being active (Pa > 0.000 or Pa = 1.000); Pi = Probability of being inactive (Pi = 0.000). Our results revealed that GNT has activity against the substrate CYP2D6 (Pa = 0.957) and CYP3A4 (Pa = 0.553) and treatment of Alzheimer’s disease (Pa = 0.428), confirming the results in silico. AChE has low pharmacokinetics in its interactions. However, since GNT is metabolized in the liver through CYP2D6 and CYP3A4, its hepatic metabolism can be affected by specific substrates, inhibitors, or enhancers of the same enzymes [25]. According to Tables 5 and 6, the biological activity of these compounds revealed that ZINC86196920, ZINC16951574 (LMQC2), and ZINC86199797 have CYP2D6 and CYP3A4 substrate activities, all with Pa ˃ 0.5. This implies they have the greatest chance of being similar to other bioactive compounds. However, compounds with Pa < 0.5 did not show any reports about their activity in the database, that is, they do not have reported activities in the literature; they were also selected for the next stages of the process when their Pa ˃ Pi, as reported by Rodrigues et al. 2013 [26]. 1000 IV Active * If swallowed = Class I: fatal (LD 50 ≤ 5). Class II: fatal (5 < LD 50 ≤ 50). Class III: toxic (50 < LD 50 ≤ 300). Class IV: harmful (300 < LD50 ≤2000). Class V: may be harmful (2000 < LD50 ≤5000), Class VI: non-toxic (LD50 > 5000). The ProTox-II [ 22 ] was used to determine the predicted oral toxicity of compounds based on 2D structure similarities with 33,000 compounds and their associated LD 50 values, organ toxicity (hepatotoxicity), and toxicological endpoints (such as cytotoxicity) (Table 4). The predicted rodent LD 50 was lowest for ZINC08342556 (LMQC5) at 1000 mg/kg with inactivity for hepatotoxicity and class IV, with most other candidates in the 362–800 mg/kg range. By comparison, the predicted LD 50 of Pro-Tox-II for GNT was 85 mg/kg and class III (Table 4). However, only the molecules ZINC86196920, ZINC16951574 (LMQC2), ZINC08342556 (LMQC5), ZINC13362890, ZINC86199797, and ZINC21657754 presented class IV—if ingested, harmful (Table 4). Compounds ZINC13108311 and ZINC15910273 presented LD 50 of 1000 mg/kg and class IV, but showed hepatotoxicity, corroborating the results obtained in the Derek of the toxicophoric alert. Additionally, the compound ZINC21657754 also showed active hepatotoxicity. 2.6. Prediction of Biological Activity (PASS) The prediction of potential biological activity was performed using PASS [ 23 ] resulting in six selected compounds. Table 5shows the biological activities that were considered: a CYP2D6 substrate, a CYP3A4 substrate, treatments for Alzheimer’s disease, and a CYP3A4 inducer. The values of Pa and Pi can vary from 0.000 to 1000, recorded as the compound’s probability of being inactive or active, respectively [ 24 ]. The biological activities obtained from the PASS were related to aspects of similarity with other bioactive substances. Molecules 2023,28, 1035 9 of 31 Table 5. Prediction of biological activity of compounds used for virtual screening using the PASS online server. Compound Pa Pi Biological Activity GNT (control) 0.957 0.003 Substrate CYP2D6 0.553 0.052 Substrate CYP3A4 0.428 0.055 Treatment of Alzheimer’s disease ZINC86196920 0.408 0.056 Substrate CYP2D6 0.413 0.085 Substrate CYP3A4 0.223 0.124 Treatment of Alzheimer’s disease ZINC16951574 (LMQC2) 0.553 0.024 Substrate CYP2D6 0.264 0.187 Substrate CYP3A4 ZINC08342556 (LMQC5) - - - ZINC86199797 0.544 0.025 Substrate CYP2D6 0.532 0.055 Substrate CYP3A4 ZINC13362890 0.191 0.165 Treatment of Alzheimer’s disease ZINC21657754 0.192 0.164 Treatment of Alzheimer’s disease Pa = probability of being active (Pa > 0.000 or Pa = 1.000); Pi = Probability of being inactive (Pi = 0.000). Our results revealed that GNT has activity against the substrate CYP2D6 (Pa = 0.957) and CYP3A4 (Pa = 0.553) and treatment of Alzheimer’s disease (Pa = 0.428), confirming the results in silico. AChE has low pharmacokinetics in its interactions. However, since GNT is metabolized in the liver through CYP2D6 and CYP3A4, its hepatic metabolism can be affected by specific substrates, inhibitors, or enhancers of the same enzymes [25]. According to Tables 5and 6, the biological activity of these compounds revealed that ZINC86196920, ZINC16951574 (LMQC2), and ZINC86199797 have CYP2D6 and CYP3A4 substrate activities, all with Pa > 0.5. This implies they have the greatest chance of being similar to other bioactive compounds. However, compounds with Pa < 0.5 did not show any reports about their activity in the database, that is, they do not have reported activities in the literature; they were also selected for the next stages of the process when their Pa > Pi , as reported by Rodrigues et al. 2013 [26]. It is important to remember that the ZINC16951574 (LMQC2) proposal has a relatively good Pa value. At the same time, Pi tends to be 0 with high relevance and exhibits activities for the substrate CYP2D6 (Pa = 0.408) and CYP3A4 (Pa = 0.413). In addition to the results obtained from the in vitro study [ 27 ], this suggests that CYP2D6 and CYP3A4 metabolize about 75% of the drug. This proposal also contains the fact that it is modeled to respect the physicochemical, steric, and electronic properties important for the purpose for which these molecules were designed. 2.7. Evaluation of Generated Metabolites and Investigation of Their Properties ADMETox 2.7.1. Metabolism Prediction of the Most Promising Compounds We used SMARTCyp [ 28 ] to evaluate the probable metabolites generated in phase I and II metabolisms undergone by the most promising compounds (ZINC16951574 (LMQC2) and ZINC8342556 (LMQC5)). The predictions for the compounds show that the following metabolic reactions occur: C-Oxidation, S-Oxidation, Aliphatic Hydroxylation, Epoxidation, and Glutathione, by the CYP 450 enzyme [29] (Figure 6and Table 7). Molecules 2023,28, 1035 16 of 31 TYR-337B, and HIS-447B participate in the most significant number of interactions with the probe’s molecules. Several other probes also form hydrogen bonds. These are the same residues that ZINC16951574 (LMQC2) and ZINC08342556 (LMQC5) interacted with in our molecular docking results (Figures 8and 9). This suggests that the main compounds selected by this criterion may have great potential to act as inhibitors of the enzyme hAChE. Molecules 2023, 28, x FOR PEER REVIEW 16 of 31 frequency of contacts between the probe molecules and each residue at the site of binding to the hAChE crystalline structure. As directed, GLY-121B, TYR-124B, GLU-202B, TYR337B, and HIS-447B participate in the most significant number of interactions with the probe’s molecules. Several other probes also form hydrogen bonds. These are the same residues that ZINC16951574 (LMQC2) and ZINC08342556 (LMQC5) interacted with in our molecular docking results (Figures 8 and 9). This suggests that the main compounds selected by this criterion may have great potential to act as inhibitors of the enzyme hAChE. Figure 11. Connection at the site of binding: Distribution of unbound interactions between residues and probes in the consensus (CS) shown from the mapping of the crystalline structure hAChE (PDB ID: 4EY6) (A) Intermolecular interactions not linked between probes and waste; (B) Hydrogen bonds, exhibiting five important interactions GLY-121B, TYR-124, GLU-202, TYR-337B, and HIS-447B. Figure 12. The FTMap identifies a cluster number per protein at the active site; the crystallographic structure of the hAChE protein is shown in the surface representation. (A) The FTMAP identified the hot spot region on the surface. (B) Two hot spots were identified and displayed on the ribbon line. 2.9. MD Simulations and Affinity Energy Calculations Figure 11. Connection at the site of binding: Distribution of unbound interactions between residues and probes in the consensus (CS) shown from the mapping of the crystalline structure hAChE (PDB ID: 4EY6) ( A ) Intermolecular interactions not linked between probes and waste; ( B ) Hydrogen bonds, exhibiting five important interactions GLY-121B, TYR-124, GLU-202, TYR-337B, and HIS-447B. Molecules 2023, 28, x FOR PEER REVIEW 16 of 31 frequency of contacts between the probe molecules and each residue at the site of binding to the hAChE crystalline structure. As directed, GLY-121B, TYR-124B, GLU-202B, TYR337B, and HIS-447B participate in the most significant number of interactions with the probe’s molecules. Several other probes also form hydrogen bonds. These are the same residues that ZINC16951574 (LMQC2) and ZINC08342556 (LMQC5) interacted with in our molecular docking results (Figures 8 and 9). This suggests that the main compounds selected by this criterion may have great potential to act as inhibitors of the enzyme hAChE. Figure 11. Connection at the site of binding: Distribution of unbound interactions between residues and probes in the consensus (CS) shown from the mapping of the crystalline structure hAChE (PDB ID: 4EY6) (A) Intermolecular interactions not linked between probes and waste; (B) Hydrogen bonds, exhibiting five important interactions GLY-121B, TYR-124, GLU-202, TYR-337B, and HIS-447B. Figure 12. The FTMap identifies a cluster number per protein at the active site; the crystallographic structure of the hAChE protein is shown in the surface representation. (A) The FTMAP identified the hot spot region on the surface. (B) Two hot spots were identified and displayed on the ribbon line. 2.9. MD Simulations and Affinity Energy Calculations Figure 12. The FTMap identifies a cluster number per protein at the active site; the crystallographic structure of the hAChE protein is shown in the surface representation. ( A ) The FTMAP identified the hot spot region on the surface. (B) Two hot spots were identified and displayed on the ribbon line. Molecules 2023,28, 1035 17 of 31 2.9. MD Simulations and Affinity Energy Calculations 2.9.1. MD Simulations The MD simulations were performed to study the binding strength and mode of action of the GNT, ZINC16951574 (LMQC2), and ZINC08342556 (LMQC5) interacting with hAChE, molecular target for DA. The trajectories analysis of MD assisted in evaluating the conformation changes with time in the protein-ligand complexes and control compounds. We have evaluated the dynamic change by means of RMSD, RMSF, and the radius of gyration (Rg) during the 100 ns simulation period. Root-mean-square deviations (RMSD) were plotted by the mean of the carbon alpha (C α ) atom of the protein backbone. The RMSD value for the protein in the complex with ligands was approximately 3 Å as shown (Figure 13). Molecules 2023, 28, x FOR PEER REVIEW 17 of 31 2.9.1. MD Simulations The MD simulations were performed to study the binding strength and mode of action of the GNT, ZINC16951574 (LMQC2), and ZINC08342556 (LMQC5) interacting with hAChE, molecular target for DA. The trajectories analysis of MD assisted in evaluating the conformation changes with time in the protein-ligand complexes and control compounds. We have evaluated the dynamic change by means of RMSD, RMSF, and the radius of gyration (Rg) during the 100 ns simulation period. Root-mean-square deviations (RMSD) were plotted by the mean of the carbon alpha (Cα) atom of the protein backbone. The RMSD value for the protein in the complex with ligands was approximately 3 Å as shown (Figure 13). Figure 13A shows that the RMSD value starts steadily increase from 1 Å to 2 Å for the simulation time scale (0 to 10 ns) and then keeps fluctuating between 1.5 Å to 2 Å during the whole simulation time. The point fluctuation RMSD of both complexes is 1.49 ± 0.17 nm. No greater fluctuation was observed during whole trajectories, which showed that ligands remained bound to the protein’s active site and are stable to compounds ZINC16951574 (LMQC2), and ZINC08342556 (LMQC5) at 1.49 Å ± 0.17. Figure 13. (A) Root-mean-square deviation RMSD plots of hAChE/GNT (control), hAChE/ZINC08342556 (LMQC5), and hAChE/ZINC16951574 (LMQC2) for 100 ns MD-simulation. (B) Root-mean-square fluctuation (RMSF) plots of hAChE/GNT (control), hAChE/ZINC08342556 (LMQC5), and hAChE/ZINC16951574 (LMQC2). (C) the radius of gyration (Rg) plots of hAChE/GNT (control), hAChE hAChE/ZINC16951574 (LMQC2), and hAChE/ZINC08342556 (LMQC5). For all figures, the red color represents ZINC16951574 (LMQC2), and the green color represents ZINC08342556 (LMQC5). Figure 13. ( A ) Root-mean-square deviation RMSD plots of hAChE/GNT (control), hAChE/ZINC 08342556 (LMQC5), and hAChE/ZINC16951574 (LMQC2) for 100 ns MD-simulation. ( B ) Rootmean-square fluctuation (RMSF) plots of hAChE/GNT (control), hAChE/ZINC08342556 (LMQC5), and hAChE/ZINC16951574 (LMQC2). ( C ) the radius of gyration (Rg) plots of hAChE/GNT (control), hAChE hAChE/ZINC16951574 (LMQC2), and hAChE/ZINC08342556 (LMQC5). For all figures, the red color represents ZINC16951574 (LMQC2), and the green color represents ZINC08342556 (LMQC5). Figure 13A shows that the RMSD value starts steadily increase from 1 Å to 2 Å for the simulation time scale (0 to 10 ns) and then keeps fluctuating between 1.5 Å to 2 Å during the whole simulation time. The point fluctuation RMSD of both complexes is 1.49 ± 0.17 nm. No greater fluctuation was observed during whole trajectories, which showed that ligands remained bound to the protein’s active site and are stable to compounds ZINC16951574 (LMQC2), and ZINC08342556 (LMQC5) at 1.49 Å ±0.17. Molecules 2023,28, 1035 18 of 31 RMSF of all the residues was calculated to investigate the local protein mobility. As seen in Figure 13B, the protein’s ligand-binding site region is more rigid than the end terminal of the domain. Ligands binding residues (GLU202, TRP86, TYR337, Gly121, Gly122, and His447) are more rigid than other protein portions. A greater fluctuation is observed in Residue GLY256, CYS257 of the LMCQ5 complex since residue 256, 257 appears in the target protein’s loop region and does not impair the interaction of the ligands’ inactive site. The radius of gyration (R g ) values over the simulation time scale are calculated for complexes and are presented in Figure 13C. It has been observed that the R g value ranges from 22.8 Å to 23.3 Å in a complex (hAChE/ZINC16951574-LMQC2) and 22.7 Å to 23.1 Å in a complex (hAChE-ZINC08342556/LMCQ5) compared with hAChE/GNT (control) and remain constant with little fluctuating between these ranges in the whole trajectory. Ligands generally showed a stable interaction profile throughout the MD simulations period. Moreover, the RMSD of only ligands was plotted to evaluate the interaction of ligands in the binding pocket (Figure 14). A peak of the RMSD of selected compounds and control is obtained between 0 Å and 1 Å, revealing that ligands (ZINC08342556-LMQC5, and ZINC16951574-LMQC2), showed a similar interaction pattern to observe for GNT (control compound) in the target protein. Molecules 2023, 28, x FOR PEER REVIEW 18 of 31 RMSF of all the residues was calculated to investigate the local protein mobility. As seen in Figure 13B, the protein’s ligand-binding site region is more rigid than the end terminal of the domain. Ligands binding residues (GLU202, TRP86, TYR337, Gly121, Gly122, and His447) are more rigid than other protein portions. A greater fluctuation is observed in Residue GLY256, CYS257 of the LMCQ5 complex since residue 256, 257 appears in the target protein’s loop region and does not impair the interaction of the ligands’ inactive site. The radius of gyration (Rg) values over the simulation time scale are calculated for complexes and are presented in Figure 13C. It has been observed that the Rg value ranges from 22.8 Å to 23.3 Å in a complex (hAChE/ZINC16951574-LMQC2) and 22.7 Å to 23.1 Å in a complex (hAChE-ZINC08342556/LMCQ5) compared with hAChE/GNT (control) and remain constant with little fluctuating between these ranges in the whole trajectory. Ligands generally showed a stable interaction profile throughout the MD simulations period. Moreover, the RMSD of only ligands was plotted to evaluate the interaction of ligands in the binding pocket (Figure 14). A peak of the RMSD of selected compounds and control is obtained between 0 Å and 1 Å, revealing that ligands (ZINC08342556-LMQC5, and ZINC16951574-LMQC2), showed a similar interaction pattern to observe for GNT (control compound) in the target protein. Figure 14. Root-mean-square deviation RMSDs plot of ligands GNT (control), ZINC08342556 (LMQC5), and ZINC16951574 (LMQC2) for 100 ns MD simulations. 2.9.2. MMPBSA Binding Free Energy Calculation MMPBSA results rank the compounds ZINC08342556 (LMQC5) < ZINC16951574 (LMQC2) < GNT (control). In control, compound electrostatics is a major energy contribution, whereas, in other compounds, ΔEvdW is dominant over other energies term. Moreover, for ranking, binding energy criteria are used as a component of the energies ΔEvdW and ΔEelectrostatic; the polar and non-polar contributions are ΔEEPB and ΔENPOLAR (Table 12 and Figure S4). ZINC16951574 (LMQC2) and ZINC08342556 (LMQC5) show stronger and more significant free energy values compared with GNT (Template). In comparison, ZINC08342556 (LMQC5) has the highest ΔEvdW, ΔEelectrostatic, and ΔEEPB values. While ZINC16951574 (LMQC2) has the lowest ∆Gtotal due to the lowest electrostatic contribution, all compounds show favorable binding interaction and have significant binding energy values with the target protein. The ΔEEPB also remained significantly unfavorable, showing that ΔEvdW has a major contribution to binding with the target protein. The total energy of the ZINC16951574 (LMQC2) and ZINC08342556 (LMQC5) is significantly different from that of control compounds GNT (control). In particular, the ZINC08342556 (LMQC5) complex shows a higher energy value of −54.2359 kcal/mol, indicating stronger interactions with receptor proteins in the active site. Other energy components include polar and non-polar solvation and van der Waals interaction forces. All ligands had stability during simulation time and showed negative binding energy values. Figure 14. Root-mean-square deviation RMSDs plot of ligands GNT (control), ZINC08342556 (LMQC5), and ZINC16951574 (LMQC2) for 100 ns MD simulations. 2.9.2. MMPBSA Binding Free Energy Calculation MMPBSAresultsrankthecompoundsZINC08342556(LMQC5)<ZINC16951574 (LMQC2) < GNT (control). In control, compound electrostatics is a major energy contribution, whereas, in other compounds, ∆ E vdW is dominant over other energies term. Moreover, for ranking, binding energy criteria are used as a component of the energies ∆ E vdW and ∆ Eelectrostatic; the polar and non-polar contributions are ∆EEPB and ∆ENPOLAR (Table 12 and Figure S4). Table 12. Binding free energy (kcal/mol) of the compounds studied. Compound ∆Gtotal (a) ∆EVDW (b) ∆Eelectrostatic (c) ∆EEPB (d) ∆EENPOLAR (e) GNT (control) −44.4665 −48.1589 −7.5486 15.7695 −4.5285 ZINC16951574 (LMQC2) −52.3996 −46.3026 −135.8650 135.2977 −5.5297 ZINC08342556 (LMQC5) −54.2359 −58.9034 −1.4877 11.1959 −5.0406 (a) ∆ G total (total free energy); (b) ∆ E VDW (van der Waals total free energy); (c) ∆ E electrostatic (electrostatics interactions); (d) ∆ E EPB (polar solvation free energy); (e) ∆ E NPOLAR (non-polar solvation free energy). GNT. Data are expressed as kcal/mol. ZINC16951574 (LMQC2) and ZINC08342556 (LMQC5) show stronger and more significant free energy values compared with GNT (Template). In comparison, ZINC08342556 (LMQC5) has the highest ∆ E vdW , ∆ Eelectrostatic, and ∆ E EPB values. While ZINC16951574 (LMQC2) has the lowest ∆ G total due to the lowest electrostatic contribution, all compounds Molecules 2023,28, 1035 19 of 31 show favorable binding interaction and have significant binding energy values with the target protein. The ∆ E EPB also remained significantly unfavorable, showing that ∆ E vdW has a major contribution to binding with the target protein. The total energy of the ZINC16951574 (LMQC2) and ZINC08342556 (LMQC5) is significantly different fromthatof control compoundsGNT(control). Inparticular, the ZINC08342556 (LMQC5) complex shows a higher energy value of − 54.2359 kcal/mol, indicating stronger interactions with receptor proteins in the active site. Other energy components include polar and non-polar solvation and van der Waals interaction forces. All ligands had stability during simulation time and showed negative binding energy values. The energy components ∆ E vdW and ∆ E electrostatic were more involved in all complexes. In contrast, ∆ E EPB polar energy value 135.865 was observed in ZINC16951574 (LMQC2), which is significantly higher than the control compound (GNT), which showed that this compound has more polar interactions within the active site of the receptor. Additionally, ∆ E NPOLAR was found to be very similar in control and screening compounds. This indicates that ZINC16951574 (LMQC2) and ZINC08342556 (LMQC5) can be strong inhibitors of the hAChE enzyme receptor. 2.10. Prediction of Lipophilicity and Water Solubility and Structure-Activity Relationship of the Promising Molecule Before proposing molecules’ theoretical synthetic routes for the compounds, we believe this methodological proposal will help medicinal chemists elaborate solutions based on their solubility, considering future in vivo tests [40]. Lipophilicity is a key parameter in drug discovery, based on the complements of the successful and informative physicochemical properties in medicinal chemistry [ 41 ]. Additionally, a drug meant for parenteral usage must be highly soluble in water to deliver a sufficient quantity of active ingredients in the small volume of such a pharmaceutical dosage [42]. Two topological methods for predicting water solubility are included in SwissADME. The first is an implementation of the ESOL36 model, and the second is an adaptation by Ali et al. 2012 [43]. The third SwissADME predictor for solubility was developed by SILICOS-IT [42]. The LogP parameter affects the ability of a molecule to decompose in non-polar versus aqueous environments [ 44 ]. GNT showed the highest consensual LogPo/w with a value of 3.33 (see Table 13), low solubility in water, and requiring organic solvents for solubilization. GNT has been reported to be soluble in organic solvents such as ethanol, DMSO, and dimethyl formamide (DMF) and is sparingly soluble in aqueous buffers [45–47]. Table 13. Lipophilicity and Water Solubility for Promising Compounds *. Compound iLOG XLOGP WLOGP MLOGP SILICOS-IT Consensus LogP GNT (control) 3.37 3.36 2.66 3.38 3.87 3.33 ZINC16951574 (LMQC2) 2.64 1.84 1.32 1.74 2.03 1.91 ZINC08342556 (LMQC5) 3.27 2.85 2.57 2.13 3.25 2.82 * iLOGP: physics-based method relying on free energies of solvation in n-octanol and water calculated by the Generalized-Born and solvent-accessible surface area model; XLOGP: an atomistic method including corrective factors and knowledge-based library; WLOGP, implementation of a purely atomistic method based on the fragmental system of Wildman and Crippen; MLOGP: an archetype of a topological method relying on a linear relationship with 13 molecular descriptors; SILICOS-IT: a hybrid method is relying on 27 fragments and 7 topological descriptors. The promising molecules investigated in the present study, the ZINC16951574 (LMQC2), presented (1.91) the lowest consensus value of LogPo/w. Then, ZINC08342556 (LMQC5) showed a value of 2.82 of LogPo/w. A range with the template molecule of ± 1.42 and 0.51, respectively. Positive values indicate that all highly lipophilic molecules meet an essential criterion for a drug candidate. To estimate water solubility, we can classify on a qualitative scale based on the log S scale: highly soluble > 0 > very soluble > − 2 > soluble > − 4 > moderately soluble > − 6 > slightly soluble > − 10 > insoluble [ 48 ]. Compounds ZINC16951574 (LMQC2) and Molecules 2023,28, 1035 20 of 31 ZINC08342556 (LMQC5) are moderately soluble (Ms) or soluble (S) in at least two solubility methods compared to GNT (control). In this study, only negative logS values in the range − 6.00 to − 6.75 were found to be promising compounds that showed consensus regarding logS values in the range, as shown in Table 14. Table 14. Prediction of solubility through the free web tool SwissADME *. Compound ESOL Ali SILICOS-IT Consensus LogS GNT (control) −4.27 (ms) −4.67 (ms) −4.82 (ms) −6.00 (ps) ZINC16951574 (LMQC2) −3.51 (s) −3.34 (s) −4.52 (ms) −6.75 (ps) ZINC08342556 (LMQC5) −2.93 (s) −2.34 (s) −2.96 (s) −6.19 (ps) * Solubility (s = Soluble, ms = Moderately soluble, ps = Poorly soluble), and Consensus LogS. At the end of the virtual screening, the two compounds with the most promising results were subjected to a search on SciFinder ® (https://scifinder.cas.org/ (accessed on 30 September 2022)). No additional information on the compounds selected in the research was found; only information on some physical and chemical properties was already reported in the ZINCPharmer database. The ZINC16951574 (LMQC2) molecule has a characteristic dimethoxy group, very similar to a new class of medicines developed as a treatment for senile dementia of Alzheimer’s type. Its action increases central cholinergic activity, inhibiting AChE (AChE) in the brain [49]. The ZINC08342556 (LMQC5) molecule has in its structure a triazolothiadiazine group, which can serve as MTDLs (New multipurpose targeted ligands), as is evident in in vitro and in vivo studies from Sagar et al. 2018 [ 50 ] Additionally, it revealed important antiamyloid, neuroprotective and anti-amnesic properties [50]. In general, the results of the present study suggest that the selected compounds can be tested for biological activities with good evidence of reproducing the results in silico. Therefore, future studies are needed to confirm the inhibitory activity of the AChE enzyme by these molecules. 2.11. Prediction of Synthetic Accessibility (SA) and Theoretical Synthetic Routes 2.11.1. Prediction of Synthetic Accessibility AMBIT and SwissADME were used to evaluate the SA of ZINC16951574 (LMQC2), and ZINC08342556 (LMQC5) (see Table 15). Table 15. Synthetic Accessibility (SA) Prediction for Selected Compounds *. Compound SA (%) (a) SA Score (%) (b) GNT (control) 64.33 40.76 ZINC16951574 (LMQC2) 36.05 40.57 ZINC08342556 (LMQC5) 56.32 40.55 * Synthetic Accessibility (SA). (a) AMBIT web server ranges from easy accessibility (score ≥ 50), median accessibility (10 < score ≤ 49), and difficult accessibility (score ≤ 10). (b) SwissADME—SA scores range from 10 (very easy) to 100 (very difficult). Only the ZINC08342556 (LMQC5) molecule was the only one that presented the synthetic accessibility predicted as easy, getting a score of 56.318. ZINC16951574 (LMQC2) showed SA, achieving a score above 36.045, indicating median accessibility for synthesis. By comparison, the accessibility prediction for GNT (control) was 64.3331. The result obtained by SwissADME for the ZINC16951574 (LMQC2) presented an SA score of 40.57%, and the ZINC08342556 (LMQC5) presented a score of 40.55% (Table 15). Compared to the template composite (GNT), which presented an SA score of 40.76%, the SA values were close, ranging between ±0.19 and ±0.21. ZINC16951574 (LMQC2), and ZINC08342556 (LMQC5) can be considered difficult to synthesize, considering the results obtained and the data found in the literature [42]. Molecules 2023,28, 1035 21 of 31 With relation to GNT for the ZINC08342556 (LMQC5). Showed a variation of ± 8.01 to AMBIT web server and ± 0.21 to SwissADME, in which we propose in this study two theoretical synthetic routes; the first proposal is the formation of 4-amino-5-(tetrahydrofuran2-yl)-4H-1,2,4-triazole-3-thiol (III) a key intermediate (Figures 15 and 16). However, this step can provide the final compound in three stages with an excellent theoretical yield; we designed a second synthetic theoretical route composed of only a step to improve the yield of ZINC08342556 (LMQC5). Molecules 2023, 28, x FOR PEER REVIEW 21 of 31 Only the ZINC08342556 (LMQC5) molecule was the only one that presented the synthetic accessibility predicted as easy, getting a score of 56.318. ZINC16951574 (LMQC2) showed SA, achieving a score above 36.045, indicating median accessibility for synthesis. By comparison, the accessibility prediction for GNT (control) was 64.3331. The result obtained by SwissADME for the ZINC16951574 (LMQC2) presented an SA score of 40.57%, and the ZINC08342556 (LMQC5) presented a score of 40.55% (Table 15). Compared to the template composite (GNT), which presented an SA score of 40.76%, the SA values were close, ranging between ±0.19 and ±0.21. ZINC16951574 (LMQC2), and ZINC08342556 (LMQC5) can be considered difficult to synthesize, considering the results obtained and the data found in the literature [42]. With relation to GNT for the ZINC08342556 (LMQC5). Showed a variation of ±8.01 to AMBIT web server and ±0.21 to SwissADME, in which we propose in this study two theoretical synthetic routes; the first proposal is the formation of 4-amino-5-(tetrahydrofuran-2-yl)-4H-1,2,4-triazole-3-thiol (III) a key intermediate (Figures 15 and 16). However, this step can provide the final compound in three stages with an excellent theoretical yield; we designed a second synthetic theoretical route composed of only a step to improve the yield of ZINC08342556 (LMQC5). 2.11.2. Theoretical Synthetic Routes Proposed for Compounds ZINC16951574 (LMQC2), and ZINC08342556 (LMQC5) The theoretical proposed synthesis of ZINC16951574 (LMQC2) (see Figure 15) starts with the formation of benzofuran acetaldehyde intermediate XI according to studies by G. Stork et al. 2009 [51]. Subsequently, protected dienyl XIII can be formed by a hydrozirconation of XII with Schwartz reagent, in situ reaction with XI, and finally trapping with chlorotrimethylsilane (TESCl). Figure 15. Theoretical synthetic route for the preparation of compound ZINC16951574 (LMQC2). Cycloaddition of XIII with decalin and triethylamine (Et3N) will produce phenanthrofuran XIV [51]. Conversion of the hydroxyl group in XIV to the methoxy group can Figure 15. Theoretical synthetic route for the preparation of compound ZINC16951574 (LMQC2). Molecules 2023, 28, x FOR PEER REVIEW 22 of 31 be carried out by standard methylation to form XV [52]. Styrene XVII will be synthesized by reacting with Wittig salt under strongly basic conditions [53] from aldehyde XVI previously furnished from ester XV. The deprotection of XVII and following Dess-Martin oxidation will allow obtaining the precursor XVIII. To conclude, reductive amination of carbonyl [54] XVIII will furnish the desired final product ZINC16951574 (LMQC2). We propose two theoretical synthetic routes for the compounds ZINC08342556 (LMQC5). First (Figure 16), the formation of 4-amino-5-(tetrahydrofuran-2-yl)-4H-1,2,4triazole-3-thiol (III) as a key intermediate is proposed. The readily available hydrazide I first reacted with carbon disulfide (CS2) in alcoholic potassium hydroxide (KOH) to give dithiocarbazate derivative II. Triazole III will then be formed from the reaction of II with hydrazide hydrate (N2H4 . H2O) under reflux conditions [55,56]. Materials I, IV, V, VI, and XII are commercially available. Finally, treating III with substituted phenacyl bromide IV under reflux and using ethanol (EtOH) as solvent [56] will produce the final triazolothiadiazine ZINC08342556 (LMQC5). Figure 16. First theoretical synthetic route for the preparation of compound LQMC 5. Although this proposed synthetic route can give the final compound in three steps with a good theoretical yield, we designed a second synthetic theoretical route consisting of only one step to improve the yield of ZINC08342556 (LMQC5) (see Figure 17). Cyclocondensation of V with substituted phenacyl bromide IV under reflux conditions in an ethanolic medium [57] will form the target compound ZINC08342556 (LMQC5). Materials I, IV, V, VI, and XII are commercially available. II Figure 17. Second synthetic theoretical route (alternative) for the preparation of the compound LQMC5. Figure 16. First theoretical synthetic route for the preparation of compound LQMC 5. 2.11.2. Theoretical Synthetic Routes Proposed for Compounds ZINC16951574 (LMQC2), and ZINC08342556 (LMQC5) The theoretical proposed synthesis of ZINC16951574 (LMQC2) (see Figure 15) starts with the formation of benzofuran acetaldehyde intermediate XI according to studies by G. Stork et al. 2009 [ 51 ]. Subsequently, protected dienyl XIII can be formed by a hydrozir- Molecules 2023,28, 1035 22 of 31 conation of XII with Schwartz reagent, in situ reaction with XI, and finally trapping with chlorotrimethylsilane (TESCl). Cycloaddition of XIII with decalin and triethylamine (Et 3 N) will produce phenanthrofuran XIV [ 51 ]. Conversion of the hydroxyl group in XIV to the methoxy group can be carried out by standard methylation to form XV [ 52 ]. Styrene XVII will be synthesized by reacting with Wittig salt under strongly basic conditions [ 53 ] from aldehyde XVI previously furnished from ester XV. The deprotection of XVII and following Dess-Martin oxidation will allow obtaining the precursor XVIII. To conclude, reductive amination of carbonyl [ 54 ] XVIII will furnish the desired final product ZINC16951574 (LMQC2). We propose two theoretical synthetic routes for the compounds ZINC08342556 (LMQC5). First (Figure 16), the formation of 4-amino-5-(tetrahydrofuran-2-yl)-4H-1,2,4-triazole-3-thiol (III) as a key intermediate is proposed. The readily available hydrazide I first reacted with carbon disulfide (CS2) in alcoholic potassium hydroxide (KOH) to give dithiocarbazate derivative II. Triazole III will then be formed from the reaction of II with hydrazide hydrate (N 2 H 4 . H 2 O) under reflux conditions [ 55 , 56 ]. Materials I, IV, V, VI, and XII are commercially available. Finally, treating III with substituted phenacyl bromide IV under reflux and using ethanol (EtOH) as solvent [ 56 ] will produce the final triazolothiadiazine ZINC08342556 (LMQC5). Although this proposed synthetic route can give the final compound in three steps with a good theoretical yield, we designed a second synthetic theoretical route consisting of only one step to improve the yield of ZINC08342556 (LMQC5) (see Figure 17). Cyclocondensation of V with substituted phenacyl bromide IV under reflux conditions in an ethanolic medium [ 57 ] will form the target compound ZINC08342556 (LMQC5). Materials I, IV, V, VI, and XII are commercially available. II Molecules 2023, 28, x FOR PEER REVIEW 22 of 31 be carried out by standard methylation to form XV [52]. Styrene XVII will be synthesized by reacting with Wittig salt under strongly basic conditions [53] from aldehyde XVI previously furnished from ester XV. The deprotection of XVII and following Dess-Martin oxidation will allow obtaining the precursor XVIII. To conclude, reductive amination of carbonyl [54] XVIII will furnish the desired final product ZINC16951574 (LMQC2). We propose two theoretical synthetic routes for the compounds ZINC08342556 (LMQC5). First (Figure 16), the formation of 4-amino-5-(tetrahydrofuran-2-yl)-4H-1,2,4triazole-3-thiol (III) as a key intermediate is proposed. The readily available hydrazide I first reacted with carbon disulfide (CS2) in alcoholic potassium hydroxide (KOH) to give dithiocarbazate derivative II. Triazole III will then be formed from the reaction of II with hydrazide hydrate (N2H4 . H2O) under reflux conditions [55,56]. Materials I, IV, V, VI, and XII are commercially available. Finally, treating III with substituted phenacyl bromide IV under reflux and using ethanol (EtOH) as solvent [56] will produce the final triazolothiadiazine ZINC08342556 (LMQC5). Figure 16. First theoretical synthetic route for the preparation of compound LQMC 5. Although this proposed synthetic route can give the final compound in three steps with a good theoretical yield, we designed a second synthetic theoretical route consisting of only one step to improve the yield of ZINC08342556 (LMQC5) (see Figure 17). Cyclocondensation of V with substituted phenacyl bromide IV under reflux conditions in an ethanolic medium [57] will form the target compound ZINC08342556 (LMQC5). Materials I, IV, V, VI, and XII are commercially available. II Figure 17. Second synthetic theoretical route (alternative) for the preparation of the compound LQMC5. Figure 17. Second synthetic theoretical route (alternative) for the preparation of the compound LQMC5. 3. Materials and Methods 3.1. Compound Selection Initially, GNT was selected as a template compound and 15 structures based on studies carried out by Atanasova et al., 2015 [ 14 ] (Table 16 and Figure S1) which has inhibitory activities (IC 50 ) against the hAChE enzyme determined under the same biological test conditions, were divided into two classes according to their activities: ( − ) less active (those with IC 50 ≥ 99.1 nM) and (+) most active (those with IC 50 < 94.1 nM). The chemical structures of these selected inhibitors were designed using ChemSketch 12 and subsequently optimized using Discovery Studio 4.0. Molecules 2023,28, 1035 23 of 31 Table 16. Structures and IC50 values for the test set of GNT derivatives. Molecules 2023, 28, x FOR PEER REVIEW 23 of 31 3. Materials and Methods 3.1. Compound Selection Initially, GNT was selected as a template compound and 15 structures based on studies carried out by Atanasova et al., 2015 [14] (Table 16 and Figure S1) which has inhibitory activities (IC50) against the hAChE enzyme determined under the same biological test conditions, were divided into two classes according to their activities: (−) less active (those with IC50 ≥ 99.1 nM) and (+) most active (those with IC50 < 94.1 nM). The chemical structures of these selected inhibitors were designed using ChemSketch 12 and subsequently optimized using Discovery Studio 4.0. Table 16. Structures and IC50 values for the test set of GNT derivatives. O H N O OH (CH2)n-O - R No R n IC50, nM No R n IC50, nM 1 N 6 5.62 9 N 10 61.20 2 N 6 6.52 10 N 9 88.10 3 N 4 8.86 11 NO 4 94.10 4 N 4 12.40 12 NO 4 99.10 5 N 2 24.50 13 NO 6 113.00 6 N 8 34.70 14 N O 6 122.00 7 N O 4 53.80 15 N 6 222.00 8 N 12 58.10 - - - - 3.2. Pharmacophore Model Generation PharmaGist [58] generated the pharmacophoric pattern of GNT and 15 structures. The method essentially aligns and superimposes the molecules with the template (GNT), identifying the molecules with the greatest number of common characteristics and highest punctuation values, being considered the best candidates for pharmacophoric models [59] 3.3. Validation of the Pharmacophore Model (a) No R (b) n (c) IC50 (d), nM No R (b) N (c) IC50 (d), nM 1 Molecules 2023, 28, x FOR PEER REVIEW 23 of 31 3. Materials and Methods 3.1. Compound Selection Initially, GNT was selected as a template compound and 15 structures based on studies carried out by Atanasova et al., 2015 [14] (Table 16 and Figure S1) which has inhibitory activities (IC50) against the hAChE enzyme determined under the same biological test conditions, were divided into two classes according to their activities: (−) less active (those with IC50 ≥ 99.1 nM) and (+) most active (those with IC50 < 94.1 nM). The chemical structures of these selected inhibitors were designed using ChemSketch 12 and subsequently optimized using Discovery Studio 4.0. Table 16. Structures and IC50 values for the test set of GNT derivatives. O H N O OH (CH2)n-O - R No R n IC50, nM No R n IC50, nM 1 N 6 5.62 9 N 10 61.20 2 N 6 6.52 10 N 9 88.10 3 N 4 8.86 11 NO 4 94.10 4 N 4 12.40 12 NO 4 99.10 5 N 2 24.50 13 NO 6 113.00 6 N 8 34.70 14 N O 6 122.00 7 N O 4 53.80 15 N 6 222.00 8 N 12 58.10 - - - - 3.2. Pharmacophore Model Generation PharmaGist [58] generated the pharmacophoric pattern of GNT and 15 structures. The method essentially aligns and superimposes the molecules with the template (GNT), identifying the molecules with the greatest number of common characteristics and highest punctuation values, being considered the best candidates for pharmacophoric models [59] 3.3. Validation of the Pharmacophore Model 6 5.62 9 Molecules 2023, 28, x FOR PEER REVIEW 23 of 31 3. Materials and Methods 3.1. Compound Selection Initially, GNT was selected as a template compound and 15 structures based on studies carried out by Atanasova et al., 2015 [14] (Table 16 and Figure S1) which has inhibitory activities (IC50) against the hAChE enzyme determined under the same biological test conditions, were divided into two classes according to their activities: (−) less active (those with IC50 ≥ 99.1 nM) and (+) most active (those with IC50 < 94.1 nM). The chemical structures of these selected inhibitors were designed using ChemSketch 12 and subsequently optimized using Discovery Studio 4.0. Table 16. Structures and IC50 values for the test set of GNT derivatives. O H N O OH (CH2)n-O - R No R n IC50, nM No R n IC50, nM 1 N 6 5.62 9 N 10 61.20 2 N 6 6.52 10 N 9 88.10 3 N 4 8.86 11 NO 4 94.10 4 N 4 12.40 12 NO 4 99.10 5 N 2 24.50 13 NO 6 113.00 6 N 8 34.70 14 N O 6 122.00 7 N O 4 53.80 15 N 6 222.00 8 N 12 58.10 - - - - 3.2. Pharmacophore Model Generation PharmaGist [58] generated the pharmacophoric pattern of GNT and 15 structures. The method essentially aligns and superimposes the molecules with the template (GNT), identifying the molecules with the greatest number of common characteristics and highest punctuation values, being considered the best candidates for pharmacophoric models [59] 3.3. Validation of the Pharmacophore Model 10 61.20 2 Molecules 2023, 28, x FOR PEER REVIEW 23 of 31 3. Materials and Methods 3.1. Compound Selection Initially, GNT was selected as a template compound and 15 structures based on studies carried out by Atanasova et al., 2015 [14] (Table 16 and Figure S1) which has inhibitory activities (IC50) against the hAChE enzyme determined under the same biological test conditions, were divided into two classes according to their activities: (−) less active (those with IC50 ≥ 99.1 nM) and (+) most active (those with IC50 < 94.1 nM). The chemical structures of these selected inhibitors were designed using ChemSketch 12 and subsequently optimized using Discovery Studio 4.0. Table 16. Structures and IC50 values for the test set of GNT derivatives. O H N O OH (CH2)n-O - R No R n IC50, nM No R n IC50, nM 1 N 6 5.62 9 N 10 61.20 2 N 6 6.52 10 N 9 88.10 3 N 4 8.86 11 NO 4 94.10 4 N 4 12.40 12 NO 4 99.10 5 N 2 24.50 13 NO 6 113.00 6 N 8 34.70 14 N O 6 122.00 7 N O 4 53.80 15 N 6 222.00 8 N 12 58.10 - - - - 3.2. Pharmacophore Model Generation PharmaGist [58] generated the pharmacophoric pattern of GNT and 15 structures. The method essentially aligns and superimposes the molecules with the template (GNT), identifying the molecules with the greatest number of common characteristics and highest punctuation values, being considered the best candidates for pharmacophoric models [59] 3.3. Validation of the Pharmacophore Model 6 6.52 10 Molecules 2023, 28, x FOR PEER REVIEW 23 of 31 3. Materials and Methods 3.1. Compound Selection Initially, GNT was selected as a template compound and 15 structures based on studies carried out by Atanasova et al., 2015 [14] (Table 16 and Figure S1) which has inhibitory activities (IC50) against the hAChE enzyme determined under the same biological test conditions, were divided into two classes according to their activities: (−) less active (those with IC50 ≥ 99.1 nM) and (+) most active (those with IC50 < 94.1 nM). The chemical structures of these selected inhibitors were designed using ChemSketch 12 and subsequently optimized using Discovery Studio 4.0. Table 16. Structures and IC50 values for the test set of GNT derivatives. O H N O OH (CH2)n-O - R No R n IC50, nM No R n IC50, nM 1 N 6 5.62 9 N 10 61.20 2 N 6 6.52 10 N 9 88.10 3 N 4 8.86 11 NO 4 94.10 4 N 4 12.40 12 NO 4 99.10 5 N 2 24.50 13 NO 6 113.00 6 N 8 34.70 14 N O 6 122.00 7 N O 4 53.80 15 N 6 222.00 8 N 12 58.10 - - - - 3.2. Pharmacophore Model Generation PharmaGist [58] generated the pharmacophoric pattern of GNT and 15 structures. The method essentially aligns and superimposes the molecules with the template (GNT), identifying the molecules with the greatest number of common characteristics and highest punctuation values, being considered the best candidates for pharmacophoric models [59] 3.3. Validation of the Pharmacophore Model 9 88.10 3 Molecules 2023, 28, x FOR PEER REVIEW 23 of 31 3. Materials and Methods 3.1. Compound Selection Initially, GNT was selected as a template compound and 15 structures based on studies carried out by Atanasova et al., 2015 [14] (Table 16 and Figure S1) which has inhibitory activities (IC50) against the hAChE enzyme determined under the same biological test conditions, were divided into two classes according to their activities: (−) less active (those with IC50 ≥ 99.1 nM) and (+) most active (those with IC50 < 94.1 nM). The chemical structures of these selected inhibitors were designed using ChemSketch 12 and subsequently optimized using Discovery Studio 4.0. Table 16. Structures and IC50 values for the test set of GNT derivatives. O H N O OH (CH2)n-O - R No R n IC50, nM No R n IC50, nM 1 N 6 5.62 9 N 10 61.20 2 N 6 6.52 10 N 9 88.10 3 N 4 8.86 11 NO 4 94.10 4 N 4 12.40 12 NO 4 99.10 5 N 2 24.50 13 NO 6 113.00 6 N 8 34.70 14 N O 6 122.00 7 N O 4 53.80 15 N 6 222.00 8 N 12 58.10 - - - - 3.2. Pharmacophore Model Generation PharmaGist [58] generated the pharmacophoric pattern of GNT and 15 structures. The method essentially aligns and superimposes the molecules with the template (GNT), identifying the molecules with the greatest number of common characteristics and highest punctuation values, being considered the best candidates for pharmacophoric models [59] 3.3. Validation of the Pharmacophore Model 4 8.86 11 Molecules 2023, 28, x FOR PEER REVIEW 23 of 31 3. Materials and Methods 3.1. Compound Selection Initially, GNT was selected as a template compound and 15 structures based on studies carried out by Atanasova et al., 2015 [14] (Table 16 and Figure S1) which has inhibitory activities (IC50) against the hAChE enzyme determined under the same biological test conditions, were divided into two classes according to their activities: (−) less active (those with IC50 ≥ 99.1 nM) and (+) most active (those with IC50 < 94.1 nM). The chemical structures of these selected inhibitors were designed using ChemSketch 12 and subsequently optimized using Discovery Studio 4.0. Table 16. Structures and IC50 values for the test set of GNT derivatives. O H N O OH (CH2)n-O - R No R n IC50, nM No R n IC50, nM 1 N 6 5.62 9 N 10 61.20 2 N 6 6.52 10 N 9 88.10 3 N 4 8.86 11 NO 4 94.10 4 N 4 12.40 12 NO 4 99.10 5 N 2 24.50 13 NO 6 113.00 6 N 8 34.70 14 N O 6 122.00 7 N O 4 53.80 15 N 6 222.00 8 N 12 58.10 - - - - 3.2. Pharmacophore Model Generation PharmaGist [58] generated the pharmacophoric pattern of GNT and 15 structures. The method essentially aligns and superimposes the molecules with the template (GNT), identifying the molecules with the greatest number of common characteristics and highest punctuation values, being considered the best candidates for pharmacophoric models [59] 3.3. Validation of the Pharmacophore Model 4 94.10 4 Molecules 2023, 28, x FOR PEER REVIEW 23 of 31 3. Materials and Methods 3.1. Compound Selection Initially, GNT was selected as a template compound and 15 structures based on studies carried out by Atanasova et al., 2015 [14] (Table 16 and Figure S1) which has inhibitory activities (IC50) against the hAChE enzyme determined under the same biological test conditions, were divided into two classes according to their activities: (−) less active (those with IC50 ≥ 99.1 nM) and (+) most active (those with IC50 < 94.1 nM). The chemical structures of these selected inhibitors were designed using ChemSketch 12 and subsequently optimized using Discovery Studio 4.0. Table 16. Structures and IC50 values for the test set of GNT derivatives. O H N O OH (CH2)n-O - R No R n IC50, nM No R n IC50, nM 1 N 6 5.62 9 N 10 61.20 2 N 6 6.52 10 N 9 88.10 3 N 4 8.86 11 NO 4 94.10 4 N 4 12.40 12 NO 4 99.10 5 N 2 24.50 13 NO 6 113.00 6 N 8 34.70 14 N O 6 122.00 7 N O 4 53.80 15 N 6 222.00 8 N 12 58.10 - - - - 3.2. Pharmacophore Model Generation PharmaGist [58] generated the pharmacophoric pattern of GNT and 15 structures. The method essentially aligns and superimposes the molecules with the template (GNT), identifying the molecules with the greatest number of common characteristics and highest punctuation values, being considered the best candidates for pharmacophoric models [59] 3.3. Validation of the Pharmacophore Model 4 12.40 12 Molecules 2023, 28, x FOR PEER REVIEW 23 of 31 3. Materials and Methods 3.1. Compound Selection Initially, GNT was selected as a template compound and 15 structures based on studies carried out by Atanasova et al., 2015 [14] (Table 16 and Figure S1) which has inhibitory activities (IC50) against the hAChE enzyme determined under the same biological test conditions, were divided into two classes according to their activities: (−) less active (those with IC50 ≥ 99.1 nM) and (+) most active (those with IC50 < 94.1 nM). The chemical structures of these selected inhibitors were designed using ChemSketch 12 and subsequently optimized using Discovery Studio 4.0. Table 16. Structures and IC50 values for the test set of GNT derivatives. O H N O OH (CH2)n-O - R No R n IC50, nM No R n IC50, nM 1 N 6 5.62 9 N 10 61.20 2 N 6 6.52 10 N 9 88.10 3 N 4 8.86 11 NO 4 94.10 4 N 4 12.40 12 NO 4 99.10 5 N 2 24.50 13 NO 6 113.00 6 N 8 34.70 14 N O 6 122.00 7 N O 4 53.80 15 N 6 222.00 8 N 12 58.10 - - - - 3.2. Pharmacophore Model Generation PharmaGist [58] generated the pharmacophoric pattern of GNT and 15 structures. The method essentially aligns and superimposes the molecules with the template (GNT), identifying the molecules with the greatest number of common characteristics and highest punctuation values, being considered the best candidates for pharmacophoric models [59] 3.3. Validation of the Pharmacophore Model 4 99.10 5 Molecules 2023, 28, x FOR PEER REVIEW 23 of 31 3. Materials and Methods 3.1. Compound Selection Initially, GNT was selected as a template compound and 15 structures based on studies carried out by Atanasova et al., 2015 [14] (Table 16 and Figure S1) which has inhibitory activities (IC50) against the hAChE enzyme determined under the same biological test conditions, were divided into two classes according to their activities: (−) less active (those with IC50 ≥ 99.1 nM) and (+) most active (those with IC50 < 94.1 nM). The chemical structures of these selected inhibitors were designed using ChemSketch 12 and subsequently optimized using Discovery Studio 4.0. Table 16. Structures and IC50 values for the test set of GNT derivatives. O H N O OH (CH2)n-O - R No R n IC50, nM No R n IC50, nM 1 N 6 5.62 9 N 10 61.20 2 N 6 6.52 10 N 9 88.10 3 N 4 8.86 11 NO 4 94.10 4 N 4 12.40 12 NO 4 99.10 5 N 2 24.50 13 NO 6 113.00 6 N 8 34.70 14 N O 6 122.00 7 N O 4 53.80 15 N 6 222.00 8 N 12 58.10 - - - - 3.2. Pharmacophore Model Generation PharmaGist [58] generated the pharmacophoric pattern of GNT and 15 structures. The method essentially aligns and superimposes the molecules with the template (GNT), identifying the molecules with the greatest number of common characteristics and highest punctuation values, being considered the best candidates for pharmacophoric models [59] 3.3. Validation of the Pharmacophore Model 2 24.50 13 Molecules 2023, 28, x FOR PEER REVIEW 23 of 31 3. Materials and Methods 3.1. Compound Selection Initially, GNT was selected as a template compound and 15 structures based on studies carried out by Atanasova et al., 2015 [14] (Table 16 and Figure S1) which has inhibitory activities (IC50) against the hAChE enzyme determined under the same biological test conditions, were divided into two classes according to their activities: (−) less active (those with IC50 ≥ 99.1 nM) and (+) most active (those with IC50 < 94.1 nM). The chemical structures of these selected inhibitors were designed using ChemSketch 12 and subsequently optimized using Discovery Studio 4.0. Table 16. Structures and IC50 values for the test set of GNT derivatives. O H N O OH (CH2)n-O - R No R n IC50, nM No R n IC50, nM 1 N 6 5.62 9 N 10 61.20 2 N 6 6.52 10 N 9 88.10 3 N 4 8.86 11 NO 4 94.10 4 N 4 12.40 12 NO 4 99.10 5 N 2 24.50 13 NO 6 113.00 6 N 8 34.70 14 N O 6 122.00 7 N O 4 53.80 15 N 6 222.00 8 N 12 58.10 - - - - 3.2. Pharmacophore Model Generation PharmaGist [58] generated the pharmacophoric pattern of GNT and 15 structures. The method essentially aligns and superimposes the molecules with the template (GNT), identifying the molecules with the greatest number of common characteristics and highest punctuation values, being considered the best candidates for pharmacophoric models [59] 3.3. Validation of the Pharmacophore Model 6 113.00 6 Molecules 2023, 28, x FOR PEER REVIEW 23 of 31 3. Materials and Methods 3.1. Compound Selection Initially, GNT was selected as a template compound and 15 structures based on studies carried out by Atanasova et al., 2015 [14] (Table 16 and Figure S1) which has inhibitory activities (IC50) against the hAChE enzyme determined under the same biological test conditions, were divided into two classes according to their activities: (−) less active (those with IC50 ≥ 99.1 nM) and (+) most active (those with IC50 < 94.1 nM). The chemical structures of these selected inhibitors were designed using ChemSketch 12 and subsequently optimized using Discovery Studio 4.0. Table 16. Structures and IC50 values for the test set of GNT derivatives. O H N O OH (CH2)n-O - R No R n IC50, nM No R n IC50, nM 1 N 6 5.62 9 N 10 61.20 2 N 6 6.52 10 N 9 88.10 3 N 4 8.86 11 NO 4 94.10 4 N 4 12.40 12 NO 4 99.10 5 N 2 24.50 13 NO 6 113.00 6 N 8 34.70 14 N O 6 122.00 7 N O 4 53.80 15 N 6 222.00 8 N 12 58.10 - - - - 3.2. Pharmacophore Model Generation PharmaGist [58] generated the pharmacophoric pattern of GNT and 15 structures. The method essentially aligns and superimposes the molecules with the template (GNT), identifying the molecules with the greatest number of common characteristics and highest punctuation values, being considered the best candidates for pharmacophoric models [59] 3.3. Validation of the Pharmacophore Model 8 34.70 14 Molecules 2023, 28, x FOR PEER REVIEW 23 of 31 3. Materials and Methods 3.1. Compound Selection Initially, GNT was selected as a template compound and 15 structures based on studies carried out by Atanasova et al., 2015 [14] (Table 16 and Figure S1) which has inhibitory activities (IC50) against the hAChE enzyme determined under the same biological test conditions, were divided into two classes according to their activities: (−) less active (those with IC50 ≥ 99.1 nM) and (+) most active (those with IC50 < 94.1 nM). The chemical structures of these selected inhibitors were designed using ChemSketch 12 and subsequently optimized using Discovery Studio 4.0. Table 16. Structures and IC50 values for the test set of GNT derivatives. O H N O OH (CH2)n-O - R No R n IC50, nM No R n IC50, nM 1 N 6 5.62 9 N 10 61.20 2 N 6 6.52 10 N 9 88.10 3 N 4 8.86 11 NO 4 94.10 4 N 4 12.40 12 NO 4 99.10 5 N 2 24.50 13 NO 6 113.00 6 N 8 34.70 14 N O 6 122.00 7 N O 4 53.80 15 N 6 222.00 8 N 12 58.10 - - - - 3.2. Pharmacophore Model Generation PharmaGist [58] generated the pharmacophoric pattern of GNT and 15 structures. The method essentially aligns and superimposes the molecules with the template (GNT), identifying the molecules with the greatest number of common characteristics and highest punctuation values, being considered the best candidates for pharmacophoric models [59] 3.3. Validation of the Pharmacophore Model 6 122.00 7 Molecules 2023, 28, x FOR PEER REVIEW 23 of 31 3. Materials and Methods 3.1. Compound Selection Initially, GNT was selected as a template compound and 15 structures based on studies carried out by Atanasova et al., 2015 [14] (Table 16 and Figure S1) which has inhibitory activities (IC50) against the hAChE enzyme determined under the same biological test conditions, were divided into two classes according to their activities: (−) less active (those with IC50 ≥ 99.1 nM) and (+) most active (those with IC50 < 94.1 nM). The chemical structures of these selected inhibitors were designed using ChemSketch 12 and subsequently optimized using Discovery Studio 4.0. Table 16. Structures and IC50 values for the test set of GNT derivatives. O H N O OH (CH2)n-O - R No R n IC50, nM No R n IC50, nM 1 N 6 5.62 9 N 10 61.20 2 N 6 6.52 10 N 9 88.10 3 N 4 8.86 11 NO 4 94.10 4 N 4 12.40 12 NO 4 99.10 5 N 2 24.50 13 NO 6 113.00 6 N 8 34.70 14 N O 6 122.00 7 N O 4 53.80 15 N 6 222.00 8 N 12 58.10 - - - - 3.2. Pharmacophore Model Generation PharmaGist [58] generated the pharmacophoric pattern of GNT and 15 structures. The method essentially aligns and superimposes the molecules with the template (GNT), identifying the molecules with the greatest number of common characteristics and highest punctuation values, being considered the best candidates for pharmacophoric models [59] 3.3. Validation of the Pharmacophore Model 4 53.80 15 Molecules 2023, 28, x FOR PEER REVIEW 23 of 31 3. Materials and Methods 3.1. Compound Selection Initially, GNT was selected as a template compound and 15 structures based on studies carried out by Atanasova et al., 2015 [14] (Table 16 and Figure S1) which has inhibitory activities (IC50) against the hAChE enzyme determined under the same biological test conditions, were divided into two classes according to their activities: (−) less active (those with IC50 ≥ 99.1 nM) and (+) most active (those with IC50 < 94.1 nM). The chemical structures of these selected inhibitors were designed using ChemSketch 12 and subsequently optimized using Discovery Studio 4.0. Table 16. Structures and IC50 values for the test set of GNT derivatives. O H N O OH (CH2)n-O - R No R n IC50, nM No R n IC50, nM 1 N 6 5.62 9 N 10 61.20 2 N 6 6.52 10 N 9 88.10 3 N 4 8.86 11 NO 4 94.10 4 N 4 12.40 12 NO 4 99.10 5 N 2 24.50 13 NO 6 113.00 6 N 8 34.70 14 N O 6 122.00 7 N O 4 53.80 15 N 6 222.00 8 N 12 58.10 - - - - 3.2. Pharmacophore Model Generation PharmaGist [58] generated the pharmacophoric pattern of GNT and 15 structures. The method essentially aligns and superimposes the molecules with the template (GNT), identifying the molecules with the greatest number of common characteristics and highest punctuation values, being considered the best candidates for pharmacophoric models [59] 3.3. Validation of the Pharmacophore Model 6 222.00 8 Molecules 2023, 28, x FOR PEER REVIEW 23 of 31 3. Materials and Methods 3.1. Compound Selection Initially, GNT was selected as a template compound and 15 structures based on studies carried out by Atanasova et al., 2015 [14] (Table 16 and Figure S1) which has inhibitory activities (IC50) against the hAChE enzyme determined under the same biological test conditions, were divided into two classes according to their activities: (−) less active (those with IC50 ≥ 99.1 nM) and (+) most active (those with IC50 < 94.1 nM). The chemical structures of these selected inhibitors were designed using ChemSketch 12 and subsequently optimized using Discovery Studio 4.0. Table 16. Structures and IC50 values for the test set of GNT derivatives. O H N O OH (CH2)n-O - R No R n IC50, nM No R n IC50, nM 1 N 6 5.62 9 N 10 61.20 2 N 6 6.52 10 N 9 88.10 3 N 4 8.86 11 NO 4 94.10 4 N 4 12.40 12 NO 4 99.10 5 N 2 24.50 13 NO 6 113.00 6 N 8 34.70 14 N O 6 122.00 7 N O 4 53.80 15 N 6 222.00 8 N 12 58.10 - - - - 3.2. Pharmacophore Model Generation PharmaGist [58] generated the pharmacophoric pattern of GNT and 15 structures. The method essentially aligns and superimposes the molecules with the template (GNT), identifying the molecules with the greatest number of common characteristics and highest punctuation values, being considered the best candidates for pharmacophoric models [59] 3.3. Validation of the Pharmacophore Model 12 58.10 - - - - (a) The structures of the compounds, (b) Derivatives with indole moiety ending the side chain attached to N-atom, (c) Substituent fixed as an indole moiety ending group in combination with various C substituents, (d) Inhibitory concentration-50%. 3.2. Pharmacophore Model Generation PharmaGist [ 58 ] generated the pharmacophoric pattern of GNT and 15 structures. The method essentially aligns and superimposes the molecules with the template (GNT), identifying the molecules with the greatest number of common characteristics and highest punctuation values, being considered the best candidates for pharmacophoric models [ 59 ] 3.3. Validation of the Pharmacophore Model The pharmacophoric models indicated by PharmaGist were characterized according to their physical-chemical and structural properties, such as Atoms (ATM), Spatial characteristics (SF), Aromatic (ARO), Hydrophobic (HYD), and Hydrogen-bond acceptor (ACC). Pearson’s correlation was used to identify the relationship between the pharmacophoric properties of the 15 structures associated with the values of pIC 50 . The correlation cutoff point was 0.3, according to previous studies by Ferreira et al., 2019 [ 13 ] and Silva Costa et al., 2018 [ 60 ]. The experimental values of IC 50 were converted into pIC 50 ( − log IC 50 ) to reduce the inconsistencies caused by the statistical steps. Molecules 2023,28, 1035 24 of 31 Hierarchical Cluster Analysis (HCA) was also used to verify and group the pharmacophoric models based on their chemical similarities. Euclidean distance was used to measure similarity. We used the most active structures (5.62 nM to 94.1 nM) and least active (99.1 nM to 222 nM) in the function of pIC50 (nM). 3.4. Pharmacophoric-Based Hierarchical Virtual Screening The pharmacophoric model was submitted to ZINCpharmer [ 16 ] to discover new potential ligands, and the following databases were used: Zinc_Purchasable, Zinc_Drug Database, Zinc_In_Man, Zinc_Natural_Derivatives, Zinc_Natural_Products and Zinc_FDA_BindingD of the ZincPharmer database. The screening process was carried out using the filters: (a) a number of different orientations of the same conformation (max hits per conf): 1; (b) a number of different orientations of different conformations of the same molecule (max hits per mol): 1; (c) molecular weight: 200 to 500 DA; (d) top hit 200; and (d) number of rotatable bonds in the molecule: 1 to 10. 3.5. In Silico Evaluation of Pharmacokinetic and Toxicological Properties of Promising Compounds ADMETox and hepatotoxicity properties were evaluated using software such as QikProp, PreADMET, Derek, and ProTox-II. The pharmacokinetic properties were analyzed with QikProp: #star, Linpiski rule, human intestinal oral absorption; QPPCaco; QPPMDCK; QPlogPo/w; CNS; QPlogBB and MW. PreADMET calculates pharmacokinetic properties such as Human Intestinal Absorption (HIA), the permeability of CaCO 2 cells in vitro (P CaCO2 ), the permeability of Maden Darby Canine Kidney cell (P MDCK ), skin permeability (P Skin ), plasma protein binding (PPB) and penetration through the blood/brain barrier (C Brain /C Blood ), and toxicological properties such as mutagenic and carcinogenic effects. The toxicological properties were predicted using the Derek software 10.0. This software is based on the search for the 2D similarity of the molecule in question or its fragments with fragments from a database with already recognized toxicity (toxicophore). It correlates the structures with various toxicological characteristics, such as mutagenicity, carcinogenicity, skin sensitization, irritation, reproductive and developmental toxicity, neurotoxicity, etc. The compounds had their hepatotoxicity evaluated using ProTox-II [ 22 ]. It allows the enabling the prediction of the largest number of toxicity endpoints consisting of 33 models and their associated LD50 values. 3.6. Prediction of Biological Activity The prediction of the biological activity of the compounds was performed using the PASS software [ 23 ], which predicts with a high accuracy of up to 2000 biological activities that are possible for chemical compounds. The approach used in the PASS is based on the suggestion that the chemical structure of compounds is closely related to their biological activity. The software was used to predict the potential biological activity of hAChE ligands against Alzheimer’s and cytochromes 3A4 and 2D6. 3.7. Molecular Docking Simulations and FTMap 3.7.1. Docking Simulations The crystallographic structure of hAChE interacting with GNT was retrieved from the Protein Data Bank (PDB), where it can be located with the PDB ID: 4EY6 [ 61 ] (see Table 17). Table 17. Protocol data used in the validation of molecular docking. Receptor Ligand Coordinates of the Grid Center (Angstrom) Grid Dimensions (Angstrom) hAChE (PDB ID: 4EY6) GNT X = 8.817 Y = −60.624 Z = −23.964 X = 18.522 Y = 23.934 Z = 11.629 Molecules 2023,28, 1035 25 of 31 Molecular docking calculations were performed using AutoDock 4.2/PyRx [ 62 , 63 ] with parameters of the genetic algorithm (with a population size of 150), a maximum number of evaluations of 250,000, a maximum number of generations of 27,000, and a crossing rate of 0.8. Interactions between inhibitors and hAChE were visualized using Discovery Studio 4.0. The calculation of the binding affinity ( ∆ G) was also carried out to compare the actual data collected and the values obtained in silico, which was the same protocol adopted by Santos et al., 2020 [14], according to Equation (1). ∆G = −RT ln Ki(1) where R (gas constant) is 1.987 × 10 −3 kcal · mol −1· K −1 , T (temperature) is 303 K, and Ki (inhibition constant) is 61.96 ×10−9M for GNT (reported in the PDB) [61]. The best free energy of binding values was obtained in the PyRx tool GUI and log files. The calculation of free energy of binding that we will try to estimate is calculated according to Morris et al., 2009 [64] and Lai et al., 2012 [65] (Equation (2)). ∆Gbind =∆H−T∆S (2) ∆ H represents the enthalpic contribution, T ∆ S represents the entropic contribution, and T is the temperature in Kelvin. ∆S are connected via the Gibbs equation. 3.7.2. Active Site Mapping—Fragment Binding Hot Spots The FTMap [ 66 ] identifies protein druggable hot spots using Fourier domain correlation techniques; the FTMAP method consists of accelerated MD simulations to calculate and equilibrate the structure on which the surface is interacting with a series of probes. For the FTMap analysis, the crystalline structure hAChE in the complex with GNT was submitted and analyzed in the protein interaction mode, according to the published protocol [66,67]. 3.8. Metabolites and Their Properties ADMETox and Carcinogenicity The compounds ZINC16951574 (LMQC2) and ZINC08342556 (LMQC5) were investigated for the generation of metabolites; ADMETox and carcinogenic properties of these metabolites using SMARTCyp [28], PreADMETox [68], and Metatox [69]. 3.9. MD Simulations The MD simulations were run using MDWeb, based on simulation software such as NAMD, Amber, and Gromacs. The complexes that were submitted to MD simulations were obtained via the molecular docking of the ligands GNT, LMCQ5, and LMCQ2 with hAChE as molecular targets. The parameters of the ligands were defined by GAFF [ 70 ], proteins were then treated with FF14SB [ 71 ], and the ionization state of their amino acids was studied with PROPKA [ 72 ] at a neutral pH. The complexes were solvated in an octahedral periodic box dimension of 12 Å defined by the TIP3P [ 73 ] water model. The partial charges of the complexes were neutralized by adding counter ions. The energy of the systems was mimicked in five stages, where each one was executed 3000 cycles using the steepest descent method and 5000 cycles using the conjugate gradient algorithm. In the first stage, the hydrogen of the water molecules was minimized, and in the second stage, the water molecules and ions were minimized soon after the hydrogen atoms of the protein were minimized, and finally, all the solvents and solutes had their energy minimized. The systems were heated gradually to 300 K for 800 ps and were performed using an NVT ensemble. The collision frequency used was 3 ps −1, and the Langevin thermostat was used for temperature control [ 74 ]. The particle mesh Ewald method [ 75 ] was used to calculate electrostatic interactions, and the bonds involving hydrogen atoms were restricted by means of the SHAKE algorithm [76].