The medicinal use of Cannabis sativa L. dates back to ancient Chinese and Egyptian civilizations [1], with many compounds (more than 500) divided into flavonoids, terpenoids, and phytocannabinoids [2]. Δ-9-tetrahydrocannabinol (THC) and cannabidiol (CBD) are the most abundant natural molecules produced by cannabis [3]. Additionally, the plant produces more than 120 derivatives, generally called minor cannabinoids, in small quantities [4,5]. Recent research has shown potential antinociceptive [6], anti-inflammatory [7], and dermatological [8] activities of minor cannabinoids and other cannabinoid analogs, called cannabinoid-like molecules.

Cannabinoid-like molecules are described as compounds with cannabimimetic properties. These molecules are found in diverse plants and may share structural similarities with cannabinoids and minor cannabinoids. Examples include desmodianones and machaeriols [4].

An important field in medicinal cannabis research is the drug design for neurodegenerative diseases such as multiple sclerosis, Parkinson’s, and Alzheimer’s [9,10,11]. The multitarget effect of cannabinoid molecules is increasingly relevant for the design of novel medicines for neurodegenerative diseases. For example, in Alzheimer’s disease, some of these compounds can competitively inhibit the human acetylcholinesterase enzyme (hAChE) and bind endocannabinoid receptors, producing a neuroprotective effect [12,13,14,15]. Consequently, theoretical studies of the electronic structure of these compounds could provide insights into their reactivity, ligand–receptor interactions, and therapeutic potential [16]. By using machine learning and data mining tools that combine molecular and electronic structure properties, it is possible to develop new medicines affordably.

Machine learning (ML) methods aim to minimize the cost function, meaning they reduce the difference between the variable predicted by the model and the observed variable [17]. Deep learning (DL) methods are machine learning algorithms that involve constructing so-called deep artificial neural networks, which have several applications, such as classifying molecules, generating new molecules, and predicting physicochemical properties [18,19]. For instance, in computational biology, the successful development of the AlphaFold algorithm has enabled accurate prediction of protein structures [20]. Additionally, in organic synthesis, ML algorithms have been developed to predict more favorable routes in retrosynthetic reactions [21].

Machine learning methods are widely used in drug design [22,23,24], encompassing the discovery of new potential drugs, the prediction of biological activity, and the determination of molecular descriptors with significant biological activity. These ML models may take inputs such as in-line molecular simplifications, e.g., simplified molecular-input line-entry system (SMILES) notation, pharmacokinetic property descriptors (e.g., absorption, distribution, metabolism, excretion, and toxicity—ADMET), and electronic structure descriptors (e.g., HOMO and LUMO energies, gap values, and atomic charges). It is important to normalize or scale the data to ensure the accuracy and integrity of ML models [25], especially when the input data have different measurement units. In most cases, the input data present considerably different units of measurement, as mentioned in the examples above.

We used classical AChEIs with well-established clinical profiles, including tacrine (THA), the first FDA-approved AChEI [26], despite its hepatotoxicity concerns; huperzine A (HUP), a natural alkaloid from Huperzia serrata with neuroprotective properties [27]; huprine X (HUX) and huprine W (HUW), synthetic hybrids combining structural features of tacrine and huperzine A [28,29]; galantamine (GNT), a natural phenanthrene alkaloid that also modulates nicotinic receptors [30]; physostigmine (PHY), a reversible carbamate hAChE inhibitor derived from Physostigma venenosum [31]; and ladostigil (LADO), a dual hAChE-monoamine oxidase inhibitor with neuroprotective properties [32]. Additionally, we included the natural hAChE substrate acetylcholine (ACh) and dopamine (DOP) as neurochemical references.

Donepezil (DNP) was chosen as a key reference compound because of its favorable pharmacokinetic profile and current status as a first-line treatment for Alzheimer’s disease. We employed two machine learning algorithms to analyze molecular property relationships among cannabinoids, cannabinoid-like compounds, and these classical AChEIs: similarity matrix and K-means clustering. We implemented these algorithms specifically for the data structure to identify patterns within molecules of the same class and between different classes. Additionally, we used principal component analysis (PCA) to identify the molecular descriptors that most significantly contributed to the pharmacophoric profiles of the studied molecules. All algorithms were integrated into a Python package version 3.9 called MolSimEx (short for Molecular Similarity Explorer).

The frontier molecular orbital values and electronic descriptors are provided in Tables S3 and S4 of the Supplementary Material. The molecular orbital energy values of cannabinoids and cannabinoid-like molecules did not significantly differ between the maximum and minimum values of each orbital of the studied molecules. The most significant difference was found in the LUMO orbital (1.70 eV) between D9CT (−0.17 eV) and DMDD (−1.87 eV), while the smallest difference was found in the LUMO+4 orbital (0.49 eV), between D9CT (0.47 eV) and RDLH (−0.02 eV), excluding the energy values of classical drugs. The slight difference between the highest and lowest values for the HOMO-4 orbital (1.45 eV) between DMDE (−6.68 eV) and CNBA (−8.13 eV) is similar to the difference observed for the LUMO orbital, reflecting the similar electronic distribution of the molecules. Owing to the electronic similarities among cannabinoids, distinguishing classes without machine learning tools is challenging. This topic is particularly important, as ES descriptors are relevant in modeling the pharmacophoric profile and defining drug classes [33].

The general descriptors and ADMET data (Tables S5 and S6) are highly dispersed, as shown in a similarity matrix depicting the ES and ADMET property vectors for each molecule, expressed as percentiles. The similarity matrix is shown in the Supplementary Material (Figures S1–S3). Figure 1 shows the correlation between the ADMET and ES data spaces across the set of molecules.

When the data are normalized (Figure 1), outliers can be interpreted as having significant positive or negative correlations (ρ) for molecules with the highest and lowest property vectors. Therefore, unique outliers above the upper limits, such as those observed in the DNP for ADMET data and in the THA for ES data, correspond to the autocorrelation values of these molecules in the similarity matrix. It is worth noting that HUW exhibits a substantial increase in correlation within the ES data space compared with the ADMET data space. Conversely, molecules such as GNT and MCRA show lower correlation in the ES data space than in the ADMET data space. In contrast, others, such as RDLJ, maintain a nearly constant correlation with both ADMET and ES data.

Figure S6 displays the results of applying the K-means algorithm to the ADMET and ES data separately, using a correlation matrix, Mcorr, as the data space, with k = 5. The feature space illustrated in the image is based on the first row (x-axis) and the last row (y-axis) of the dataset (Tables S3–S6 of the Supplementary Material).

A comparative analysis between ADMET and ES data, using k = 5, revealed significant differences (Figures S6 and S7 of the Supplementary Material) in the positioning of specific molecules within the data space, suggesting that similarities between molecules may change depending on the data type considered. For instance, OMCNB is positioned close to DNP in the ADMET data space, whereas it is close to THA in the ES data space. The ADMET data showed an even distribution of cluster members, whereas the ES data showed a concentration of molecules primarily in clusters 1 and 4 (Figure S6). Cluster 3 in the ES data has increased considerably, driven by the inclusion of cannabinoids such as THC, THCP, HD11, and D9CT, along with ACh. Moreover, the THA in the ES data was migrated to cluster 2. Although MCDC and CNBA remained grouped in cluster 1, other molecular groups varied between the datasets. For example, HUP clustered with cannabinoids such as THCVM2, TCNT, and CMG, as well as the classical AChEI THA, which is in cluster 2 of the ES data. At the same time, classical inhibitors HUX and HUW co-clustered with cannabinoid-like molecules DMDD and DMDE, along with the prenyl bibenzyls RDLH, RDLI, and RDHJ. These associations suggest potential targets for drug repositioning studies.

A set of molecules showed unexpected clustering patterns when the ADMET and ES data spaces were separated, with k = 5. Therefore, the next step was to use ADMET and ES data simultaneously for the K-means algorithm. For instance, separating RDLA from other prenyl bibenzyls was only possible when ADMET and ES data were combined, with k≤5, as shown in Table S7 and Figure 2. Figure S8 in the Supplementary Material displays the results of clustering using both ADMET and EE data.

As illustrated in Figure S9 of Supplementary Material, the Pearson similarity matrix has correctly grouped the classical acetylcholinesterase inhibitors (AChEIs) based on their ADMET and ES data. For example, there is a positive correlation of 0.51 between the alkaloid cholinomimetics physostigmine and galantamine. Furthermore, the synthetic compounds HUX and HUW, derived from a combination of tacrine and huperzine, also show a positive correlation of 0.57.

Table 1 shows the similarity index higher than or equal to 0.75 for the ADMET, ES, and ADMET + ES data. As expected, acetylcholine (ACh) shows no significant positive correlations with minor cannabinoids and cannabinoid analogs when analyzed with ADMET data (as illustrated in Figure S2 of Supplementary Material).

The largest vectors in loadings (Figure S6 of Supplementary Material) for the first five PCs were HOMO-4 and HOMO-3 (0.28) for PC1, LUMO+3 (0.36) for PC2, TPSA (0.47) for PC3, charge on heteroatoms and volume (0.36) for PC4, and drugscore (0.67) for PC5.

The distribution of molecules in the ADMET and ES data space can be visualized using the first three principal components (PCs) of the molecular space, as illustrated in Figure 3. This visualization uses bi-dimensional scatter plots generated from PCA, following dimensionality reduction through loading analysis. In Figure 3, each class is represented by a distinct color and symbol. Figure 3 illustrates the molecular distribution in the reduced-dimensional space of the remaining descriptors after PCA dimensionality reduction. The diagonal panels display the density distributions of each cluster projected onto individual principal components. Higher-variance principal components (particularly the first component) provide enhanced cluster separation, as evidenced by the approximately Gaussian distributions of individual clusters along these dimensions and the reduced overlap between different clusters.

The properties retained after PCA—such as molecular weight (MW(g/mol)), number of aromatic rings, blood–brain barrier (BBB) penetration (binary descriptor), and binding score (kcal/mol) (Table 2)—demonstrate a strong ability to cluster most of the molecules examined effectively.

The Elbow and PCA methods were applied to evaluate the k value in the K-means algorithm. The elbow method was applied to the K-means algorithm, using separate ADMET and ES properties, while PCA was applied to all properties. In the PCA approach, the value of k was assumed to be equal to the number of principal components (PCs) that collectively contributed to almost all of the data variance (k = 5), as shown in Figure S4.

The elbow method is typically used to identify the ideal number of clusters in the K-means algorithm [34]. The value of k obtained by the elbow method was k=5 for ADMET and ES data, with the maximum value of the WSS being greater (>2000) for the ADMET data than for the ES data, as shown in Figure S5. Therefore, using the variance from the dimensionality reduction of the principal components (PCs) to determine k (Figure S4) is a valid approach, which is supported by the established mathematical relationship between PCA and k-means clustering, where both methods optimize similar objective functions based on data variance minimization [35]. Specifically, PCA looks for directions (principal components) that maximize variance (minimize reconstruction error). Otherwise, k-means minimizes within-cluster variance. The eigenvalue spectrum from PCA shows a natural clustering structure, as principal components with higher eigenvalues correspond to directions of greater data separability, directly informing optimal cluster number selection.

Figure 2 suggests that the clusters up to group 4 share similarities, regardless of whether k = 5 or k = 6 for the ADMET + ES data in Table S7. In contrast to the formation of these groups, ACh is isolated from the other molecules, which is significant because it serves as a functional negative control in this study. It means that ACh does not inhibit hAChE since it is the natural substrate of the enzyme. This outcome demonstrates that using combined ADMET and ES data resulted in a clustering profile that aligns more closely with the expected outcomes. Moreover, the known molecules used in the treatment of Alzheimer’s are grouped in cluster 0.

However, DNP, HUX, and HUW are grouped with RDLJ despite k = 5 or k = 6, indicating a promising possibility of RDLJ, a cyclopropanochroman derivative from Radula javanica with less psychoactive effects [36]. The RDLJ has the same methoxymethyl ether group also found in DNP. Another outcome is that ARAC1 and ARAC2 are grouped together with HUW, HUX, and DNP for k = 5. ARAC1 and ARAC2 are prenylated stilbenoids with the same skeleton as stilbene and other stilbenoids, which are linked to the treatment of Alzheimer’s disease using natural products [37,38].

The highest correlation value observed for Ach, cannabinoid, and cannabinoid-like molecules was 0.33 with TCNT (Figure S4 in the Supplementary Material). However, when examining the similarity matrix using ES data (Figure S3 of Supplementary Material), ACh exhibited positive correlations with THC (0.45), THCP (0.45), THCVM2 (0.40), D9CT (0.55), HD11 (0.60), and R6E (0.57), as well as an increased correlation with TCNT (0.54). Interestingly, the correlation between hAChE inhibitors (AChEIs) decreases significantly when only the ADMET data is taken into account. These differences result in a different cluster classification in the K-means algorithm, as shown in Figure S8 for ADMET + ES properties.

The similarity index shown by the ADMET + ES (Table 1) groups most similar cannabinoid and cannabinoid-like molecules, while separating GNT, PHYSO, and RIVA. However, all cannabinoid and cannabinoid-like molecules have a positive correlation, albeit small, with the classical AChEIs in the ADMET + ES data. For this set, the largest values are found for THCMV2 (0.35 with GNT), RDLH (0.37 with HUW), RDLJ (0.41 with GNT), and TCNT (0.41 with DOP). The observed clustering of cannabinoid-like compounds with classical AChEIs can be rationalized through shared pharmacodynamic and structural features that may influence acetylcholinesterase inhibition.

Structurally, both classical AChEIs and several cannabinoid derivatives possess aromatic systems and hydrogen bond acceptor/donor groups, where the presence of such groups facilitates binding to the acetylcholinesterase active site [39,40]. Furthermore, cannabinoids have demonstrated cholinergic modulation, e.g., in the direct AChE inhibition, as reported for THC [41], and the indirect cholinergic enhancement via CB1 receptor-mediated modulation of acetylcholine release [42,43]. The prenylated compounds, especially those derived from Radula species, share structural motifs with established AChEIs, including phenolic hydroxyl groups and extended conjugated systems; such groups may contribute to enzyme–substrate interactions [44]. These similarities support the computational clustering results, suggesting that the identified correlations may reflect pharmacological relationships rather than statistical artifacts.

As illustrated in Figure 3, this is particularly suitable in the apparent grouping of ACh in class 3 (depicted by the triangle symbol with a salmon color palette). Moreover, a distinct separation is observed between class 0, which contains classical inhibitors like DNP, GNT, and HUP, and both class 2, which consists of cannabinoids such as CMG, CNBRPL, and D9CT (represented by the grey diamond symbol), and class 1, which includes cannabinoid-like molecules such as DCBA, DDAA, GLDC, and OMCNB (marked by the light blue square symbol). From Figure 3, we can also see that the molecules in cluster 4 (marked by a red inverted triangle), which consist of cannabinoid-like compounds such as ARAC1, ARAC2, and DMDD, as well as the molecules in cluster 2, are closely grouped with the classic AChEIs in cluster 0. These findings emphasize the need to treat these properties as essential descriptors of molecular characteristics, highlighting their pivotal role in molecular analysis.

Although the loading analysis retained only ADMET properties in our final model (Table 2), it is worth noting that all electronic structure descriptors exhibited high loadings (Figure S10, Supplementary Material). This feature indicates that electronic descriptors also contributed significantly to the principal components, even though they were not retained in the final dimensionally reduced model.

Furthermore, the boxplot analysis (Figure 1) showed that electronic structure data produced more positive and negative outliers than ADMET descriptors, as indicated by the Pearson correlation coefficient. These outliers represent molecules that exhibited high positive or negative correlations with the reference molecule in each boxplot analysis. This suggests that electronic structure descriptors can identify molecular similarity and dissimilarity within our dataset, thereby separating cannabinoids from classical AChE inhibitors.

The complementary nature of these descriptor types is a challenge. While ADMET properties provided the most consistent dimensionality reduction (as shown by the loadings analysis), electronic descriptors exhibit a broader distribution of correlations, suggesting that both are valuable for molecular characterization.

Future studies could use these two approaches—dimensionality reduction via loading analysis and correlation coefficient variation analysis—by developing more robust search algorithms, such as simplex optimization or Monte Carlo simulations, to better capture the combined effects of both ADMET and electronic structure descriptors in molecular clustering.

We have identified a set of 30 molecules, out of 253 derived from phytocannabinoids, flavonoids, and terpenoids, as potential new inhibitors against the hAChE (6O4W) enzyme. We employed two machine learning algorithms to categorize the molecules: similarity matrix and K-means. These algorithms were implemented specifically for the data structure to identify patterns within molecules of the same class and between different classes. Principal component analysis was also used to find the most determining molecular descriptors.

From the K-means classification, one of the prenyl bibenzyl molecules, the cyclopropanochromane derivative radulanin J, presented the highest correlation with the AChEI profile of the drug donepezil (0.41, based on AMDET+ES data). The molecules tetrahydrocannabivarine-M2, radunalin H, radunalin J, and trans-cannabitriol showed positive correlations with the classical AChEIs in the ADMET similarity matrix and electronic structure data. The TCNT-DOP and RDLJ-DNP pairs showed a Pearson correlation coefficient of 0.41, indicating a significant correlation driven by the use of ADMET and electronic structure data.

It is important to note that, despite the promising results, particularly for radulanin J, the present computational framework includes virtual screening, requiring future experimental studies to validate the proposed insights, such as in vitro AChE inhibition assays for direct activity confirmation and IC50 determination, cytotoxicity studies using neuronal cell lines to assess safety profiles, and blood–brain barrier permeability assays to evaluate central nervous system penetration potential.

Additionally, while our rigid docking protocol provided consistent ranking criteria for the machine learning pipeline, future studies focusing on the most promising candidates would benefit from more sophisticated computational approaches, including induced-fit docking, ensemble receptor conformations, explicit water modeling, and enhanced scoring methods, such as the molecular mechanics energies combined with the generalized Born and surface area continuum solvation (MM/GBSA) calculation, to reduce potential docking bias and provide more accurate binding predictions. These validations would provide essential pharmacological data for advancing the developed model.

Moreover, in the current MolSimEx implementation, the ASED metric is effectively equivalent to the L2 norm for clustering, yielding assignments and centroids that are indistinguishable from those of standard K-Means. Accordingly, potential advantages over the classical L2 metric have yet to emerge. Future research that improves centroid estimation and aligns the objective function with ASED may advance MolSimEx, enabling the new norm to exceed L2 in specific cases.

Although these findings are promising, they are based on a relatively narrow dataset, and broader evaluations will be needed to validate the benefits of the metrics. The increasing availability of reliable electronic structure calculations is an open field for drug discovery. Future research should broaden these analyses to include larger and more diverse chemical libraries. This research would enable the identification of even more robust and generalizable molecular similarity patterns for the rational design of new drugs.

The authors also acknowledge the support of DPI/BCE/UnB. The authors have reviewed and edited the output and take full responsibility for the content of this publication.

The following supporting information can be downloaded at: https://www.mdpi.com/article/10.3390/ijms1010000/s1. References [59,100,101,102,103,104,105,106,107,108,109,110,111,112] are cited in the supplementary materials.

Conceptualization, R.C.V., É.C.M.N. and J.B.L.M.; Formal analysis, R.C.V., É.C.M.N. and J.B.L.M.; Investigation, R.C.V.; Methodology, É.C.M.N. and J.B.L.M.; Resources, J.B.L.M.; Supervision, J.B.L.M.; Visualization, R.C.V.; Writing—original draft, R.C.V.; Writing—review and editing, É.C.M.N. and J.B.L.M. All authors have read and agreed to the published version of the manuscript.

Not applicable.

Not applicable.

The source code (MolSimEx, Molecular Similarity Explorer) of this study is freely available at GitHub (https://github.com/Rafael-Campos-unb/MolSimEx) (accessed on 3 November 2025).

The authors declare no conflicts of interest.

This work was supported by the Brazilian funding agencies CNPq (Grant No. 306682/2021-4 and 308419/2025-1) and FAPDF (Grant No. 00193–00000869/2021–31).