Abstract

   Alzheimer's disease (AD) is the most known neurodegenerative disease,
   and its prevalence is predicted to increase significantly. Discovering
   novel drugs and treatments for AD is urgently needed. Drugs from
   natural products have been preferred lately due to their high potential
   and low toxicity. Citrus hystrix DC. (kaffir lime; KL) is one such
   herbal plant that is found abundantly in Southeast Asia with many
   biological activities. In this study, the potential of bioactive
   compounds from KL peel, leaf, and essential oil as anti-AD agents was
   explored using network pharmacology. First, the compounds were
   identified with KNApSAcK database and related literature. Subsequently,
   the targets of each corresponding compound were determined with SEA
   Search Server and Swiss Target Prediction, while the proteins
   associated with AD were identified using OMIM, GenCLiP3, and DisGeNET.
   Furthermore, a protein–protein interaction network and a
   compound–target interaction network were constructed to identify the
   most crucial proteins and compounds in the network by employing
   Cytoscape v3.9.1. The study continued with pathway enrichment analysis
   using STRING v1.7.1, molecular docking with PyRx and SwissDock, and
   molecular dynamics simulation with YASARA for further confirmation. Our
   results showed that almost all the secondary metabolites of KL targeted
   AD-associated genes, with oxypeucedanin and citrusoside A showing the
   highest anti-AD potential and targeting essential genes, EGFR and
   MAPK14, respectively. These targets were associated with inflammatory
   and oxidative stress pathways, indicating the potential mechanism of KL
   in attenuating AD clinical manifestation.

   Keywords: Alzheimer's disease, Citrus hystrix, protein network, network
   pharmacology, molecular docking, molecular dynamics

1. Introduction

   Alzheimer's disease (AD) is a progressive neurodegenerative disease
   with a strong correlation with age. One of the earliest symptoms of
   this disease is the impairment of cognitive functions, which may lead
   to the weakening of other vital body functions and even death. The
   prevalence of AD is estimated to increase rapidly in the middle of the
   century, affecting approximately 13,8 million Americans aged >65 years
   old [[29]1].

   AD is a multifactorial disease that could be characterized by various
   cellular and molecular processes; although the actual cause of this
   disease remains unknown, age is believed to be the highest risk factor
   [[30][2], [31][3], [32][4], [33][5], [34][6], [35][7]]. Studies have
   shown a significant correlation between aging and numerous AD
   pathogenesis [[36][10], [37][11], [38][8], [39][9]]. Several hypotheses
   have been proposed to explain the cause and mechanism of AD, such as
   the accumulation of β-amyloid (Aβ) plaques and neurofibrillary tangles
   in neurons, cholinergic neurons, neuroinflammation, and oxidative
   activity in the brain [[40]3,[41]5,[42][12], [43][13], [44][14],
   [45][15]]. Such factors cause synaptic failure and neuronal death,
   leading to brain atrophy and decreased brain functions
   [[46]3,[47]5,[48][12], [49][13], [50][14], [51][15], [52][16]].

   Till date, no pharmacological treatment has been developed that can
   prevent or stop biological activities promoting AD progression. Current
   available treatments only address symptoms and may have mild side
   effects [[53]2,[54]3,[55][13], [56][14], [57][15]]. Hence, an effective
   anti-AD is urgently needed.

   Treatment with natural products is often preferred due to its
   multitarget mechanism which enables it to target different factors that
   promote disease progression [[58][17], [59][18], [60][19], [61][20]].
   Citrus hystrix DC. (kaffir lime; KL) is a citrus plant that is abundant
   in Southeast Asia, including Indonesia. It has several beneficial
   parts, such as leaf, fruit, peel, and essential oil, which are rich in
   bioactive compounds. Previous studies have shown that KL may have some
   potential in managing AD. In an in vitro study using a neuronal
   senescent model on SH-SY5Y cells, Pattarachotanat and Tencomnao (2020)
   [[62]21] reported that the leaf and peel extracts of KL exhibit
   neuroprotective ability against senescence induced by high blood
   glucose levels and might have potential as an anti-AD agent.

   In addition to its leaf and peel, KL essential oil is used as an
   aromatherapy ingredient. Aromatherapy has been proposed as an AD
   treatment [[63][22], [64][23], [65][24], [66][25]], possibly owing to
   its capability to improve cognitive functions through olfactory
   pathways [[67]26]. Moreover, bioactive compounds of KL, such as
   sitosterol, are known to exhibit anti-AD properties. Sitoserol is
   reported to be a good inhibitor of cholinesterase (AChE), capable of
   preventing glutamate and Aβ toxicity, and possesses antioxidant
   activity [[68]21,[69][27], [70][28], [71][29]].

   However, the potential and mechanism of KL in managing AD are yet to be
   explored. For elucidating such multitarget mechanism, network
   pharmacology can be utilized to understand the system and process of
   how KL may cure AD through various pathologic pathways. The network
   pharmacology approach aims to understand how molecules act on a system
   level, impacting multiple targets and pathways simultaneously rather
   than concentrating on just one. It involves analyzing biological
   networks to gain insight into the molecular mechanisms of drugs and
   diseases [[72]30,[73]31]. This method shows promise for drug discovery,
   particularly in complex, multifactorial diseases like AD.

   This study attempts to explore the potential mechanism of KL in
   attenuating AD by applying network pharmacology to its bioactive
   compounds against AD targets. In addition, molecular docking and
   molecular dynamics simulations are performed to validate the network
   analysis results. Our findings indicated that oxypeucedanin and
   citrusoside A have potential as inhibitors of EGFR and MAPK14, which
   could potentially improve symptoms of AD by suppressing related
   inflammatory and oxidative stress pathways.

2. Materials and methods

2.1. Identification of bioactive compounds

   The bioactive compounds of KL leaf, peel, and essential oil (from
   blossom, leaf, branch, and fruit peel) were identified from the
   KNApSAck database ([74]http://www.knapsackfamily.com/KNApSAcK/)
   [[75]32], Dr. Duke database
   ([76]https://phytochem.nal.usda.gov/phytochem/search), and previous
   reports [[77]21,[78]33]. A total of 69 compounds were combined and
   deduplicated, and the secondary metabolites were selected. Finally, 52
   secondary metabolites were included in the analysis.

2.2. Identification of compound-targeted genes and AD-associated genes

   The SMILES of each identified secondary metabolite was collected from
   PubChem database ([79]https://pubchem.ncbi.nlm.nih.gov/) [[80]34] and
   used for target identification using SEA Search Server
   ([81]https://sea.bkslab.org/) [[82]35] and Swiss Target Prediction
   (STP) ([83]http://www.swisstargetprediction.ch/) [[84]36] with the Homo
   sapiens filter. The top three targets with the highest scores were
   selected from each database. The top 10 genes in cases where all the
   identified genes had the same score were chosen. All the results were
   combined and deduplicated. AD-associated genes were also collected from
   DisGeNET database ([85]https://www.disgenet.org/) [[86]37], Online
   Mendelian Inheritance in Man (OMIM) ([87]https://omim.org/) [[88]38],
   and GenCLiP3 database ([89]http://ci.smu.edu.cn/genclip3/analysis.php)
   [[90]39] with the Homo sapiens filter, combined, and deduplicated. Both
   lists of genes were entered to Venny 2.1. to find the overlapping
   genes, generating 64 KL target genes that were also AD-associated
   genes.

2.3. Construction of a Protein–Protein Interaction (PPI) network

   The overlapping genes were inserted to STRING v1.7.1. integrated with
   Cytoscape v3.9.1 [[91]40]. as organism Homo sapiens to visualize the
   interaction of proteins. The top 10 most essential genes in the network
   were identified using another application integrated with Cytoscape,
   CytoHubba v0.1 [[92]41] through the maximal clique centrality (MCC)
   algorithm. The redder the color of the node, the higher its MCC scores.

2.4. Construction of a Compound–Target Interaction (CTI) network

   A CTI network was constructed manually from the target identification
   results using Cytoscape v3.9.1. to map the interaction between protein
   targets and their corresponding targeting compounds. The target nodes
   were colored as yellow to red according to their MCC score, and the
   compound nodes were yellow. The edges show the association between the
   compound and its targets.

2.5. Analysis of KEGG pathway enrichment

   The previously constructed PPI network was further analyzed by Kyoto
   Encyclopedia of Genes and Genomes (KEGG) pathway enrichment using
   STRING Enrichment App integrated with Cytoscape. When retrieving the
   functional enrichment, the background network was set to be the genome,
   providing a list of pathways along with their corresponding genes,
   number of background genes, P-value, and false discovery rate (FDR; Q)
   value. The pathways were selected and filtered with Q-value <0.001.
   Additionally, the extent of pathway enrichment was determined by
   computing the rich factor, which was obtained by dividing the number of
   genes in a specific pathway within the network by the number of
   background genes. Subsequently, data from these 14 highly enriched
   pathways was consolidated and presented visually using Datylon
   ([93]https://www.datylon.com).

2.6. Preparation of proteins and ligands

   Molecular docking was conducted to validate the results of the network
   analysis. Epidermal growth factor receptor (EGFR) (PDB ID: [94]7B85;
   deposited by [[95]42], [[96]43]), mitogen-activated protein kinase 14
   (MAPK14) (PDB ID: [97]6HWU; deposited by [[98]44], [[99]45]), and
   peroxisome proliferator activated receptor alpha (PPARA) (PDB ID:
   [100]6KXY; deposited by [[101]46], [[102]47]) were selected as the
   proteins. The selection was made based on the results of target
   identification, MCC scores, and involvement in the most enriched
   pathway. Their 3D structures were searched in the Uniprot
   ([103]https://www.uniprot.org) [[104]48] and RCSB PDB database
   ([105]https://www.rcsb.org/). The protein's crystal structure was
   chosen based on various factors including the role of the native ligand
   bound to the protein, methodology, resolution, completeness of the
   protein chain or domain, deposition dates, and presence of any absent
   residues or mutations. The structural integrity was additionally
   validated using Ramachandran Plot statistics by PROCHECK accessed
   through PDBsum
   ([106]https://www.ebi.ac.uk/thornton-srv/databases/pdbsum/), with only
   structures exhibiting a value exceeding 90 % being ultimately selected
   for docking. The selected structures were then prepared with PyMOL v2.0
   [[107]49] by removing water and nonstandard (ligand) molecules. If any
   missing residues were found, the structure was repaired using Dock Prep
   in Chimera v1.1.6 [[108]50].

   The bioactive compounds targeting these three proteins based on the
   previous identification target results were selected as the ligands,
   and the 3D structures of these bioactive compounds were obtained from
   PubChem as.sdf files. On the other hand, the controls for each protein
   were their respective native ligands obtained from the RCSB database.
   The chosen controls were propan-2-yl
   2-[[4-[2-(dimethylamino)ethyl-methyl-amino]-2-methoxy-5-(propanoylamino
   )phenyl]amino]-4-(1-methylindol-3-yl)pyrimidine-5-carboxylate (R28) for
   EGFR,
   3-(2,5-dimethoxyphenyl)-∼{N}-[4-[4-(4-fluorophenyl)-2-[(∼{E})-phenyldia
   zenyl]-1,3-thiazol-5-yl]pyridin-2-yl]propanamide (GE5) for MAPK14, and
   6-ethyl-1-(4-fluorophenyl)-3-pentan-3-yl-pyrazolo[3,4-b]pyridine-4-carb
   oxylic acid (T06) for PPARA.

   The protein and ligand conformations were additionally prepared using
   PyRx and SwissDock, correspondingly. In the case of PyRx, the ligand
   structures were prepared in OpenBabel integrated with PyRx v0.9.8
   [[109]51]. by minimizing the energy with the Universal Force Field
   (UFF) using the conjugate gradient algorithm. The total number of steps
   was defined as 200 and the number of steps for update was set to 1. In
   addition, the procedure was configured to stop if the energy difference
   dropped below 0.1 kcal/mol. After performing ligand energy
   minimization, the proteins and compounds were inserted into the
   AutoDock Wizard integrated with PyRx. This involved designating them as
   macromolecules and ligands respectively. This step automatically added
   hydrogen molecules, merged non-polar hydrogens, and added Gasteiger
   charges to the structures, converting them into.pdbqt files ready for
   docking.

   In the case of SwissDock ([110]https://www.swissdock.ch/) [[111]52],
   Attracting Cavities 2.0 [[112]53] engine was chosen for docking. Ligand
   data was transformed into Mol2 files through PyMOL v2.0 and then
   uploaded onto the platform for processing. SwissParam with the
   MMFF-based method was utilized for parameterization. PDB IDs were
   provided for the targets, selecting the A chain in cases where multiple
   chains existed, and specifying no heteroatom. The preparation of
   targets was carried out to facilitate their utilization in CHARMM
   simulations.

2.7. Molecular docking: AutoDock Vina (PyRx)

   Targeted molecular docking was performed for each protein in AutoDock
   Vina integrated with PyRx v0.9.8. The exhaustiveness and the mode were
   set to 50 and 9, respectively. Prior to the actual docking using the
   selected ligands, a validation was conducted by redocking the proteins
   with each corresponding native ligand (control). The redocking process
   was carried out repeatedly until the RMSD (Root Mean Square Deviation)
   result was less than 2 Å, demonstrating its conformational similarity
   with the co-crystallized control and confirming the accuracy of the
   docking procedure [[113][53], [114][54], [115][55], [116][56]]. The
   RMSD values were calculated using DockRMSD
   ([117]https://zhanggroup.org/DockRMSD/) [[118]57].

   Afterward, molecular docking was performed with the prepared proteins
   and ligands using identical grid box coordinates employed for redocking
   as shown in [119]Table S5. Subsequently, the best docking poses were
   combined with their respective protein structures utilizing PyMOL
   v2.5.2. These complexes were further analyzed and visualized with
   Discovery Studio v21.1.0.20298 [[120]58] to determine the binding types
   and residues involved in the interactions.

2.8. Molecular docking: Attracting Cavities 2.0 (SwissDock)

   The same proteins and ligands were also subjected to molecular docking
   using another software known as Attracting Cavities 2.0, which is
   accessible through the SwissDock web-based server. Once the ligand and
   target were prepared, a search space was defined using the identical
   grid box coordinates employed for docking with Vina. Additionally,
   parameters were configured based on default settings, such as setting
   the number of RIC to 1, choosing medium for sampling exhaustivity, and
   choosing buried for cavity prioritization. The scoring function
   comprises the AC score, which considers CHARMM force field energy and
   fast analytical continuum treatment of solvation (FACTS) terms, as well
   as the SwissParam score, a weighted sum of polar and nonpolar terms.
   The AC score was employed for selecting the optimal model, whereas the
   SwissParam score was computed to assess the affinity of various
   compounds towards a specific target. The docking outcomes were accessed
   using ViewDock in Chimera v1.16, and the most favorable samples were
   stored as.pdb files for further examination in Discovery Studio
   v21.1.0.20298.

2.9. Molecular dynamics simulation

   The top docked complexes, EGFR-oxypeucedanin and MAPK14-citrusoside A,
   were further investigated by performing molecular dynamics simulation
   in YASARA v21.12.19 [[121]59] for a period of 10.00 ns with a timestep
   of 2.5 fs. The md_runfast.mcr script was used to define parameters such
   as the AMBER14 force field, temperature of 310K, pH of 7.4, NaCl
   concentration of 0.9 %, and pressure at 1 atm. The MD simulation was
   also performed on the control complexes for comparison purposes. The
   findings were assessed through the execution of the md_analyze.mcr
   script to gain insights into the RMSD Cα, RMSD ligand conformation, and
   RMSD ligand movement.

3. Results

3.1. Identification of compound-targeted genes and AD-associated genes

   The KL compounds were gathered from KNApSAck, Dr. Duke, and a few
   research papers [[122]21,[123]33]. The KNApSAck database yielded 45
   compounds from volatiles of blossom extract [[124]60] and fruit peel
   extract [[125]61]. In contrast, the Dr. Duke database did not yield any
   results. Warsito et al. (2017) [[126]33] reported 28 compounds derived
   from KL leaf, branch, and fruit peel essential oils. Finally,
   Pattarachotanat & Tencomnao (2020) [[127]21] documented a total of 27
   compounds combined from the leaf and fruit peel extracts. These
   compounds were combined and deduplicated, resulting in a total of 69
   compounds as displayed in [128]Table S1. From this pool of compounds,
   52 secondary metabolites were derived due to their diverse
   pharmacological activities.

   The targets of each compound were then identified using SEA Search
   Server and STP, which generated 127 target genes ([129]Table S2).
   AD-associated genes were collected from three databases (DisGeNET,
   OMIM, and GenCLiP3), resulting in 4647 genes after deduplication. A
   Venn diagram was created using Venny 2.1. to find the intersection of
   both lists of genes, resulting in 64 KL-targeted genes that were also
   associated with AD ([130]Fig. 1.).

Fig. 1.

   [131]Fig. 1
   [132]Open in a new tab

   Venn diagram of compound-targeted genes (collected from SEA Search
   Server and Swiss Target Prediction) and Alzheimer-associated genes
   (collected from DisGeNET, OMIM, and GenCLiP3).

3.2. Protein–Protein Interaction (PPI) network

   The 64 intersecting genes were inputted to STRING v1.7.1. integrated
   with Cytoscape v3.9.1. to construct the PPI network. The network
   consisted of 64 nodes, 135 edges, and 4.426 average neighbors. All
   proteins had interactions with each other except for NLRP1, TACR2, and
   MSR1. The MCC algorithm in CytoHubba v0.1., which is also integrated
   with Cytoscape, was then used to find the most crucial genes in the
   network. The redder the color of the node, the higher the MCC score and
   the greater the importance of the node in the network. The top 10 most
   essential genes were identified based on MCC as follows: EGFR, STAT3,
   TRPV1, PRKCA, CNR1, MAPK14, PPARA, PRKCE, AR, and FAAH ([133]Fig. 2.).
   The exact scores are listed in [134]Table S3.

Fig. 2.

   [135]Fig. 2
   [136]Open in a new tab

   PPI network constructed with Cytoscape v3.9.1. Each node represents the
   overlapping proteins from the previous step, and the edges show the
   interaction of the nodes. The thick red color of the node reveals the
   node's importance in the network based on its MCC score. (For
   interpretation of the references to color in this figure legend, the