A strategy to find novel candidate anti-Alzheimer’s disease drugs by constructing interaction networks between drug targets and natural compounds in medical plants

Data collection and preprocessing

More than 60,000 natural compounds from 8,529 different plants were from the TCM Database@Taiwan (http://tcm.cmu.edu.tw/). This web-based database is the most comprehensive non-commercial database of TCM(Chen, 2011). We obtained the 3D structures of the compounds from the database as mol2 files and converted them to the pdbqt format and SMILES string using Open Babel toolbox v2.3.1(O’Boyle et al., 2011).

The validated therapeutic AD target proteins and pathways were provided by Thomson Reuters Integrity database (https://integrity.thomson-pharma.com/integrity/). Studies on the target proteins were obtained using Pubmed and PMC (https://www.ncbi.nlm.nih.gov/pubmed). The target proteins structures were obtained from the Protein Data Bank (PDB) database (http://www.rcsb.org/pdb/home/). We selected target proteins of interest by determining whether studies relating them to AD and crystal protein structures (with ligand present) were available. Inhibitors or agonists were used to confirm the ligands in the references we identified. We ultimately selected 30 target proteins. Information regarding the 30 target proteins, such as the Integrity name, corresponding name in the Uniprot database (http://www.uniprot.org/), Uniprot ID, PDB ID, ligand ID and document IDs in Pubmed and PMC are available in Table S1. The 3D structures of the proteins are available as files in the pdb format. This format was converted to the pdbqt format using AutoDock tools v1.5.6 (Morris et al., 2009), and the 3D view was generated by Discovery Studio v3.1 (http://accelrys.com/products/collaborative-science/biovia-discovery-studio/).

Molecular docking between natural compounds and AD target proteins

Docking is tantamount to position the ligand in different orientations and conformations within the binding site to calculate optimal binding geometries and energies. Interactions between natural compounds and the AD target proteins were previously predicted using AutoDock Vina 1.1.2 (Trott & Olson, 2010). The docking binding site center for each target protein is the structural binding center of the ligand present in the crystal structure. The ligands were confirmed by the studies we identified. The coordinates of the docking center, ligand ID in PDB database and supporting documents are shown in Table S1. To allow free rotation of the compounds, the search space was set to 25 × 25 × 25 Å in each axis. The default settings were used for all of the other docking parameters. Each docking was performed by a command that contained the space size and three-dimensional coordinate of the docking center. For each compound, the binding pose with the lowest energy for each docking test was considered the best binding mode. A lower energy score indicated a stronger binding affinity between the ligand and receptor. The compounds with the top 0.5% docking score were chosen as the candidate ligands for each target protein.

Validation of the docking results

To validate the docking results, three methods were used. First, we manually checked the docking results and visualized the interaction between the compound and receptor to verify that the compound was in the binding pocket. Second, the original ligand in the crystal structure was set as the reference. The docking energy of candidate anti-AD compounds should be better than or comparable to that of original ligand. Third, we check whether the top 0.5% TCM compounds are similar to the existing drugs that have been reported in anti-AD research.

Interactions among the target proteins, compounds and plants

The interaction between each target protein and its best-binding TCM compound was displayed using the PyMOL (PyMOL Molecular Graphics System, version 1.7) program (https://pymol.org). The pharmacophore was displayed using Discovery Studio v3.1.

The target-compound and target-plant networks were constructed using Cytoscape v3.4.0 (Shannon et al., 2003). In these networks, the target protein and compound were considered to be connected if the compound successfully docked to the target protein, and the target protein and the plant were considered to be connected if the plant with the compound was able to interact with the target protein. The strength of the links is represented by the line’ thickness, which indicates the number of compounds shared between the target protein and plant.

Collection of anti-AD plants from Chinese medicine prescription

The term “senile dementia” was searched in the subject column of the CNKI database (http://www.cnki.net/) to retrieve Chinese medicine prescriptions for anti-AD from the relevant Chinese articles. Articles from clinical studies between 2007 and 2017 were selected. Chinese medicine prescriptions and the usage frequency were also obtained from these articles. The common anti-AD plants in traditional Chinese clinical medicines were identified from the prescriptions. For reference, Table S2 shows information regarding the Chinese version of the raw prescription data with corresponding English, the Latin name of the anti-AD plants in each prescription, the patient number (male and female if available), the article title of the study, published data (years) and the article format.

Similarity between candidate compounds and existing drugs

The Tanimoto coefficient (Tc) and Pybel (O’Boyle, Morley & Hutchison, 2008) Python package were used to measure the structural similarities between compounds. The fingerprint FP2 implemented in Pybel was generated for each structure and used to calculate Tc. Tc is defined as Tc = C(i, j)/U(i, j), where C(i, j) is the number of common features in the fingerprints of molecules i and j and U(i, j) is the number of all of the features in the union of the fingerprints of molecules i and j. If the fingerprints of two compounds are Tc = 1, even if they differ among themselves by isolated instances of C, N or O atoms, they will be considered structurally identical.

Cytoscape v3.4.0 was used to construct a network linking the candidate compounds, their target proteins and structurally identical drugs. A natural compound and an existing drug in the DrugBank database (Wishart et al., 2006) were considered to be connected if their Tc score was 1. The natural compounds and their target proteins were also connected in this network.

Clusters of potential candidate compounds for AD

One-thousand-four-hundred-seventy-six compounds located in the top 0.5% of all compounds that docked with 30 target proteins were regarded as potential candidate compounds for AD. The cluster ligands protocol in BOVIA Pipeline Pilot V8.5 (http://accelrys.com/products/collaborative-science/biovia-pipeline-pilot/) was used to cluster the compounds. During clustering, a set of compounds was assigned to different clusters based on the similarity of their properties. Clustering was based on the root-mean-square (RMS) difference of the descriptor properties or the Tanimoto distance for fingerprints. In our study, we clustered the compounds based on the Tanimoto distance using the fingerprint FP2. Cluster ligands were performed by the number of size or the number of molecules per cluster. The default parameter in Pipeline Pilot V8.5 was set to fixed number of 10 clusters. For simplicity, we used default parameters to cluster our compounds. The molecule with the lowest total distance to all other members of the cluster was considered the cluster center.

ADMET and logP properties for candidate compounds for AD

The ADMET and logP properties of candidate compounds for AD were estimated using Discovery Studio. ADMET refers to absorption, distribution, metabolism, excretion and toxicity and logP refers to the logarithm of the partition coefficient. The ADMET properties, including aqueous solubility, blood brain barrier penetration (BBB), human intestinal absorption (HIA), plasma protein binding (PPB) and hepatotoxicity, as well as logP were used to filter the compounds. The values of these properties were set as the controlled parameters, which were 3∼4 (3: good; 4: optimal) for aqueous solubility, 1∼2 (1: high; 2: medium) for BBB, 0 (0: good) for HIA, FALSE for both PPB and hepatotoxicity and logP <5.

Molecular docking of natural compounds and embedded ligands to the 30 AD target proteins

Of the 60,000 compounds in the TCM Database, 30,438 contain plant information, and these compounds were docked with the 30 selected AD target proteins. The original ligands of each target protein were also docked to their corresponding target proteins. The docking results are shown in Table 1. The docking scores of ligands embedded in the protein crystal structure ranged from −3.31 to −12.65 (kcal/mol). The lowest docking energy scores for the 30 target proteins ranged from −8.44 to −14.5 (kcal/mol). Some target proteins, such as Caspase-3, QC, IDO and GLP-1R, were able to bind over 20,000 natural compounds with docking scores superior to those of their embedded ligands. However, the docking scores of the target RAR with natural compounds were inferior to those of its embedded ligand.

Because many TCM compounds can bind to AD target proteins, for each target protein, only the top 0.5% of compounds in terms of docking scores (a total of 1,476 compounds) were selected as candidate compounds for AD. The original study of AutoDockVina showed that the success rate of Vina is 80% (RMSD < 2) for an independent validation dataset. Furthermore, Vina achieves a low standard error of 2.85 kcal/mol compared with the experimental free energies. There is a highly positive correlation between the predicted and experimental free energies of binding. In our study, almost all of the docking energies of the top 0.5% of compounds bound to target proteins were superior to those of their embedded ligands (Fig. 1). Thus, the 1,476 compounds are likely candidate compounds for AD.

Analysis of the interactions between target proteins and ligands including TCM compounds, original ligands and approved AD drugs

The docking pose interactions between the target proteins and ligands including their best-binding TCM compounds and original ligands are shown in Fig. 2A and Figs. S1 to S29. The figures show side-by-side-comparisons of best TCM ligands and original ligands. The best TCM ligands and original ligands are located in the same binding pocket for each target protein and they all have some common residues. Taking AchE as an example: the best-binding TCM compound for AchE is 24,829, the original ligand is Huperzine A, and their common residues are TYR-124, PHE-297, PHE-338, TYR-337, ASP-74. Furthermore, we compared the structures of the best-binding TCM compounds and previously known ligands (Table S3) and found the structures to be different. The Tc scores are between 0.06 and 0.48. Therefore, most of the best-binding TCM compounds are novel scaffolds.

We analyzed the pharmacophore of SIRT1 using the top 10 TCM compounds binding the target protein (Fig. 2B). The pharmacophore model consists of one hydrogen bond acceptor (HBA, green) and five hydrophobic centers (blue); therefore, we think that compounds that have this model may easily bind with SIRT1.

We also compared the docking energy scores between the approved AD drugs and the candidate compounds for the protein AchE (Table S4). The approved AD drugs that we analyzed were Donepezil, Galantamine and Rivastigmine. All the three drugs were AchE inhibitors. The average docking energy score of the candidate compounds was superior to that of the approved AD drugs for the protein AchE. The docking pose interaction of the three drugs and the best-binding TCM compound with the protein AchE is shown in Fig. S30. The three drugs and best-binding TCM compound were located in the same binding pocket, and their docking poses were different because their structures were different.

Analysis of single-target and multi-target compounds

Among the 1,476 candidate compounds for AD, there were 850 compounds with a single target and 626 compounds with multiple targets (see Figs. S31 and S32). The single-target compounds corresponding to each target varied widely in number. For example, in single-target networks, target SIRT1 corresponds to 30 compounds, whereas target lyn just corresponds to four compounds. The multi-target compounds were classified into 18 networks based on their corresponding target numbers, which ranged from 2 to 24. As the number of targets per compound increased, the number of compounds in that category decreased. For example, the two-target network contained 260 compounds, whereas the three-target and four-target networks contained 90 and 77 compounds, respectively. Finally, we observed that compound 24,508 could bind to 25 AD target proteins. The structure and network of compound 24,508 are shown in Fig. S32.

Candidate AD compounds and their enrichment plants

We mapped 1476 candidate AD compounds (corresponding to 30 AD target proteins) to the 334 plants. The plant numbers for each target protein ranged from 42 to 71, whereas the compound numbers for each target protein ranged from 62 to 132 (Fig. 3 and Fig. S33).

We selected 101 clinically related studies out of over 10,000 senile dementia-related articles and identified 141 anti-AD traditional Chinese plants from the clinical prescriptions in these articles. The 141 traditional Chinese anti-AD plants were classified based on their functional properties in the TCM database. Most of the 141 anti-AD plants were in the ‘Tonifying, Replenishing’ category, and the plants in this category accounted for 28.45% of all anti-AD plants (Fig. S34).

The best-associated plant for each target protein contained the greatest number of compounds capable of docking with the target protein (Table 2 and Fig. 4). The number of compounds in Table 2 was based on the top 0.5% of compounds successfully docking to each target protein. Thus, 30 target proteins corresponded to 16 best-associated plants, of which the top five plants were anti-AD traditional Chinese plants, including Panax and Morus, which corresponded to 7 and 5 target proteins, respectively, as well as Salvia, Rheum and Paeonia.

Similarities between candidate compounds and existing drugs

A structural comparison between all existing approved drugs recorded in the DrugBank and the top 0.5% of the natural compounds tested demonstrated that some compounds were identical to existing drugs (Tc = 1). The connection network among candidate compounds, existing drugs and AD target proteins was established, and the chemical structures of the compounds are also shown in Fig. 5. There were 17 candidate compounds, 14 existing drugs and 25 AD-associated target proteins in the network. The 14 drugs included Lutein (DB00137), Vitamin A (DB00162), Vitamin E (DB00163), Azelaic Acid (DB00548), Ergotamine (DB00696), Estradiol (DB00783), Menthol (DB00825), Drostanolone (DB00858), Glyburide (DB01016), Tubocurarine (DB01199), Metocurine (DB01336), Yohimbine (DB01392), Lactose (DB04465), Artemether (DB06697).

Most of the drugs listed above have been reported to have a positive effect in AD. Some studies have shown that Lutein is involved in preventing cognitive decline and decreasing the risk of AD; thus, Lutein may contribute to the treatment of AD (Kiko et al., 2012; Min & Min, 2014; Xu & Lin, 2015). Similarly, Vitamin A, Vitamin E, Estradiol, Menthol, Glyburide and Yohimbine are also considered useful in the prevention and treatment of AD (Bhadania et al., 2012; Dysken et al., 2014; Lamkanfi et al., 2009; Lan et al., 2016; Mohamd et al., 2011; Ono & Yamada, 2012; Peskind et al., 1995; Takasaki et al., 2011). Therefore, compounds with structures similar to the existing drugs may also have anti-AD function by interacting with similar target proteins.

Of the 17 candidate compounds, 11 can only bind with one target protein, whereas the rest, which are similar to the 14 drugs discussed above, can interact with more than one target protein. For example, compound 18491 which is similar to Menthol can only interact with the target Ftase and compounds 19,476 and 19,477 are similar to Tubocurarine and Metocurine, respectively, whereas compound 18,582 which is similar to Ergotamine can bind with 17 target proteins.

Structure clusters of candidate AD compounds

To compare the structural features of the candidate AD compounds, the 1,476 candidate compounds were assigned to 10 clusters. The structures of the cluster center compounds and maximal common substructure of each cluster are shown in Fig. 6. The Generate Maximal Common Substructure component must contain the proportion of the cluster molecules. The proportion was set to 0.5 to find the largest maximal common substructure that was contained in at least 50% of the cluster molecules. All of the compounds in the center of the cluster contained the carbocyclic structure, similar to the five approved drugs for AD. The cluster sizes varied, with the largest containing 464 compounds and the smallest containing only six compounds. Every cluster had a primary target protein that could better combine with the compounds in the cluster. The binding of the best member of the cluster to the primary target protein and their structures are also shown in Fig. 6.

To see if there were similar or different scaffolds able to bind each target protein, clusters in the sets of ligands that bind to each individual protein were shown in Tables S5 to S35. Each set of ligands was assigned to 10 clusters. The compounds in the center of the cluster do not have the same scaffolds for each individual protein.

Eleven candidate compounds for AD with favorable ADMET and logP properties

Most orally approved drugs have favorable druggability properties. After ADMET and logP filtering, 11 compounds from 1,476 candidate anti-AD compounds possessed favorable properties (Table 3). The 11 compounds all have good HIA and don’t have PPB and hepatotoxicity. Compound 28,468 has optimal aqueous solubility, and compound 16,167 has high BBB. Of the 11 compounds, eight were single-target compounds, and the remaining three were multi-target compounds. For example, compounds 5,868, 9,593, 10,639, 28,814 and 30,713 could only bind with one target protein (GABA(B), MAOB, PDE4, RAR and MAOA, respectively), whereas compounds 5,862, 5,869 and 16,167 shared one common target (AchE). Compound 5,863 was able to bind with three target proteins (AchE, GABA(B) and MGLUR), and compound 26,629 and 28,468 were able to interact with five target proteins. These 11 compound structures and their corresponding plants are shown in Fig. 7. The plants corresponding to seven of these compounds are regarded as anti-AD plants in TCM, including Curcuma kwangsiensis, Poria cocos, Lindera aggregate, Ophiopogon japonicus (L. f.) Ker-Gawl. and Glycyrrhiza glabra. Eleven compounds belong to different organic compound classes. Compounds 5,862, 5,863, 5,868, 5,869 belong to linear diarylheptanoids. Compound 10,639, 16,167, 26,629, 28,468, 28,814, 30,713 belongs to guanidines, pyridines, hydrophenanthrenes, angular pyranocoumarins, naphthyridines, 8-prenylated isoflavanones, respectively. We also checked the 11 comopund structures and existing research in pubchem database. Compounds 5,862, 5,863, 5,868, 5,869 have anti-inflammatory activity and compound 28,814 mainly has immunoregulatory activity. The other 6 compounds have not yet been reported for any biology activity at present.