In Silico Neuroprotective Effects of Specific Rheum palmatum Metabolites on Parkinson’s Disease Targets

Parkinson’s disease (PD) is one of the large-scale health issues detrimental to human quality of life, and current treatments are only focused on neuroprotection and easing symptoms. This study evaluated in silico binding activity and estimated the stability of major metabolites in the roots of R. palmatum (RP) with main protein targets in Parkinson’s disease and their ADMET properties. The major metabolites of RP were subjected to molecular docking and QSAR with α-synuclein, monoamine oxidase isoform B, catechol o-methyltransferase, and A2A adenosine receptor. From this, emodin had the greatest binding activity with Parkinson’s disease targets. The chemical stability of the selected compounds was estimated using density functional theory analyses. The docked compounds showed good stability for inhibitory action compared to dopamine and levodopa. According to their structure–activity relationship, aloe-emodin, chrysophanol, emodin, and rhein exhibited good inhibitory activity to specific targets. Finally, mediocre pharmacokinetic properties were observed due to unexceptional blood–brain barrier penetration and safety profile. It was revealed that the major metabolites of RP may have good neuroprotective activity as an additional hit for PD drug development. Also, an association between redox-mediating and activities with PD-relevant protein targets was observed, potentially opening discussion on electrochemical mechanisms with biological functions.


Introduction
Parkinson's disease (PD) is the second most prevalent neurodegenerative disease after Alzheimer's disease [1,2].It has been considered the world's fastest-growing neurological disorder and an economic burden for older adults because of medication costs and productivity losses [3].This disease is characterized by bradykinesia, resting tremors, muscle stiffness, speech and writing changes, postural instability, and a wide array of motor and non-motor symptoms [4][5][6].While the exact pathogenesis of PD is still unclear, primary observations reveal that such progressive disease develops from the degradation of dopaminergic neurons in the substantia nigra pars compacta (SNpc) in the basal ganglia, and the development of intraneuronal proteinaceous cytoplasmic inclusions called Lewy bodies [7].Approximately four of five PD patients are idiopathic; the others are presumed of genetic origins [4].Unfortunately, current treatments for PD are for neuroprotection and symptomatic therapy due to a vague understanding of its exact pathogenesis [5,8].Dopamine (a catecholamine) is a crucial target for treating symptoms of PD since it is the primary neurotransmitter for movement coordination and control (in the motor circuitry of the brain in the basal ganglia), as well as in the cognition and reward system [9].The death of dopaminergic neurons in SNpc consequently decreases dopamine production responsible for the motor symptoms of PD, such as its effect on the degradation of the nigrostriatal pathway responsible for sensory stimuli and movement, causing bradykinesia [9,10].Because of the altered motor circuitry in the basal ganglia, automatic postural adjustments are disrupted, causing postural instability and immobility [11].Also, thalamic activity from basal ganglia dysfunction was found to be promoted due to decreased dopamine, which regulates the thalamocortical loop, and led to abnormal cortical oscillations that cause resting tremors [12,13].While the death of dopaminergic neurons is vital in PD indications, no complete conclusions on neuronal death have been reached.However, prior studies hypothesized that neuroinflammation from oxidative stresses, metabolic alterations, mitochondrial dysfunction, the impairment of autophagy functions, and intestinal inflammation from gut dysbiosis affecting the enteric nervous system (ENS) are some of the fundamental mechanisms of neuronal death [5,[14][15][16].
Some specified mechanisms of neuronal death are the loss of mitochondrial complex I functionality responsible for dopamine production and fragmented mitochondria from mutations of coiled-coil-helix-coiled-coil-helix domain containing 2 (CHCHD2) [17][18][19].Also, dysregulations in the tryptophan (Trp)-kynurenine (KYN) metabolic system observed in PD patients potentially led to further neuroinflammation due to the immunomodulatory properties of the metabolites, as well as the accumulation of neuronal death-causing compounds [17,20].Indeed, imbalances in the Trp-KYN pathway were suspected of causing neurodegeneration due to a heightened ratio of neurotoxic to neuroprotective metabolites, such as the reported increase in quinolinic acid (QUIN) in blood plasma and 3hydroxykynurenine (3-HK) in the cerebrospinal fluid of PD patients, which was implicated to increased oxidative stress, excitotoxicity, and neuroinflammation where dopaminergic neurons were particularly vulnerable [17,21,22].Since dopamine is a precursor of other catecholamines, specifically epinephrine and norepinephrine, other non-motor functionssuch as peristalsis, heart rate, and blood pressure regulation-are diminished for PD patients [10,[23][24][25].Gnanaraj et al. [8] summarized the PD-relevant proteins of utmost concern for their symptomatic treatment and neuroprotection.One of PD's first genetic links is the overexpression of SNCA, an α-synuclein (ASN) coding gene, causing the formation of Lewy bodies from ASN aggregations [26].On the other hand, in the standard treatment of PD, monoamine oxidase isoform B (MAOB) and catechol o-methyltransferase (COMT) are most considered for their function in dopamine metabolism [27].Due to the impact of dopamine on human movement and control, its excessive decrease promotes uncontrolled and nonrigid movements.At the same time, the antagonism of the A 2A adenosine receptor (A 2A AR) reduces the wear-off effect of dopaminergic therapies [28].This receptor modulates motor circuits through inhibition of the direct pathway (motion-induction) and promotion of the indirect pathway (motion-arrest) [28,29].Current treatments of PD include, but are not limited to, levodopa and its combination drug (such as carbidopa); dopamine agonists (such as apomorphine, ropinirole, cabergoline); MAOB inhibitors (such as rasagiline and selegiline); COMT inhibitors (such as entacapone and tolcapone); anticholinergics; anti-inflammatory agents; and vitamins A, E, and C [5,30].
Natural products have always been an abundant and diverse source of hit and lead compounds with good pharmacokinetic properties, including for treating neurological diseases [31,32].Interestingly, complex mixtures of these medicines could provide curative potential to PD because of its multiple targeted approach from several potentially bioactive compounds.In traditional Chinese medicine, the roots of Rheum palmatum (RP) or Chinese rhubarb are used for various therapeutic uses, including PD symptom treatment [33].
Previously, it has been revealed that anthraquinones (specifically chrysophanol, rhein, aloe-emodin, emodin, and physcion as shown in Figure 1) are abundant in RP [34].Consequently, neurotransmitters such as dopamine, epinephrine, norepinephrine, and levodopa show redox-mediating activity for microbial fuel cells (MFCs) and quasi-reversible redox ability in cyclic voltammetry measurements [35].Also, PD medications were observed as quasi-reversible redox mediators where the electrochemical mechanism for biological function is still open for follow-up investigations [36].Since RP extracts show good electron-shuttling characteristics, the major metabolites of RP may have similar activity as the treatment mentioned [34].Indeed, prior studies have explored that emodin and chrysophanol exhibited neuroprotective properties in vitro and in vivo for PD [37].
or Chinese rhubarb are used for various therapeutic uses, including PD symptom treatment [33].Previously, it has been revealed that anthraquinones (specifically chrysophanol, rhein, aloe-emodin, emodin, and physcion as shown in Figure 1) are abundant in RP [34].Consequently, neurotransmitters such as dopamine, epinephrine, norepinephrine, and levodopa show redox-mediating activity for microbial fuel cells (MFCs) and quasireversible redox ability in cyclic voltammetry measurements [35].Also, PD medications were observed as quasi-reversible redox mediators where the electrochemical mechanism for biological function is still open for follow-up investigations [36].Since RP extracts show good electron-shuttling characteristics, the major metabolites of RP may have similar activity as the treatment mentioned [34].Indeed, prior studies have explored that emodin and chrysophanol exhibited neuroprotective properties in vitro and in vivo for PD [37].A molecular simulation study was conducted for RP's major metabolites, revealing the potential alternative therapeutic compounds or conjunctive use of RP extracts with conventional medicines for PD.In silico investigations of major compounds of RP were performed, considering neurotransmitters, and PD medicines have good MFC bioenergy amplification.This investigation evaluated their potencies to relevant PD targets and safety profiles with their binding affinity, pharmacokinetic properties, expected reactivity, and quantitative structure-activity relationship for the prediction of their inhibitory action.

Molecular Docking
Primary metabolites from RP were subjected to the ligand preparation protocol in Discovery Studio, which generated the 31 ligands presented in Table S1.Two standard drugs used for PD treatment (dopamine and levodopa) were used as positive controls and were prepared similarly to the mentioned metabolites.Validation of the protein structures of protein targets for PD revealed that the root-mean-square-difference (RMSD) of minimized structures compared to experimentally elucidated structures were no more than 2.0 Å, which indicated valid structures for docking studies (see Table S2) [38,39].Each ligand form was screened using hotspots-based docking, using the LibDock protocol available in Discovery Studio reported in Table 1.The ligand conformers for each metabolite with the highest LibDock score were selected for further analysis to CDOCKER, a CHARMm-based docking protocol.A molecular simulation study was conducted for RP's major metabolites, revealing the potential alternative therapeutic compounds or conjunctive use of RP extracts with conventional medicines for PD.In silico investigations of major compounds of RP were performed, considering neurotransmitters, and PD medicines have good MFC bioenergy amplification.This investigation evaluated their potencies to relevant PD targets and safety profiles with their binding affinity, pharmacokinetic properties, expected reactivity, and quantitative structure-activity relationship for the prediction of their inhibitory action.

Molecular Docking
Primary metabolites from RP were subjected to the ligand preparation protocol in Discovery Studio, which generated the 31 ligands presented in Table S1.Two standard drugs used for PD treatment (dopamine and levodopa) were used as positive controls and were prepared similarly to the mentioned metabolites.Validation of the protein structures of protein targets for PD revealed that the root-mean-square-difference (RMSD) of minimized structures compared to experimentally elucidated structures were no more than 2.0 Å, which indicated valid structures for docking studies (see Table S2) [38,39].Each ligand form was screened using hotspots-based docking, using the LibDock protocol available in Discovery Studio reported in Table 1.The ligand conformers for each metabolite with the highest LibDock score were selected for further analysis to CDOCKER, a CHARMm-based docking protocol.Furthermore, 35 protein-ligand complexes were refined through CDOCKER with the top hits presented in Table 2, where CDOCKER energy combined protein-ligand binding interaction energy and ligand conformational energy.Compounds above or within the dopamine and levodopa standard range were analyzed for their molecular interaction and properties.While binding these metabolites rendered them inferior to levodopa, multiple ligands inhibited many targets.In this condition, emodin exhibited the best binding activity among all protein targets.Specifically, ASN (PDB ID: 1XQ8) interacted the best with emodin (−34.4961kcal/mol), surpassing dopamine (−31.0394kcal/mol).There are only functional assumptions of ASN on dopamine regulation from its ability for synaptic vessel provision in presynaptic terminals and cytoskeletal dynamics from microtubular interaction [40,41].Interestingly, one symptomatic progression of PD concerns alternate folding of ASN forming insoluble fibrils, which further aggregates leading to neurodegeneration, and they are postulated as causative to prion-like cell-to-cell transmission of ASN misfolding [42].Nevertheless, the aggregation of ASN is widespread in numerous brain regions of the central nervous system (CNS) for PD patients.These aggregations impair numerous neural functions, such as proteasomal and lysosomal functions, potentially hindering protein clearance and promoting neuronal death [26,43].Additionally, there are preliminary investigations on ASN aggregation due to CHCHD2 mutations and interaction with 3-HK.Hence, it was suspected that interactions with ASN could inhibit the formation of Lewy bodies for decreased neurotoxicity.Nevertheless, chrysophanol (−39.7554kcal/mol) and physcion (−30.3292kcal/mol) may have a good interaction with MAOB (PDB ID: 2C65).At the same time, the bestinteracting compound for COMT (PDB ID: 3BWM) is emodin (−32.431kcal/mol).Normally, MAOB catalyzes the oxidative deamination of xenobiotic and biogenic amines, including dopamine [44].In contrast, COMT normally regulates levels of catecholamines in the brain and peripheral tissues, especially levodopa, by its inactivation through methyl group transfer from the S-adenosylmethionine to the catechol structure [27,44].For this reason, their inhibition ameliorates the motor dysfunction and fluctuations observed in PD patients and offsets the deterioration of dopaminergic neurons.Also, emodin has superior binding activity (−57.5427kcal/mol) to A 2A AR (PDB ID: 3EML), comparable to levodopa (−61.8193kcal/mol).At the same time, rhein (−39.4557kcal/mol) has a similarity in binding affinity with dopamine (−39.3845kcal/mol) to A 2A AR (PDB ID: 3EML).Due to the inhibitory interaction of A 2A AR with the dopamine D2 receptor, reduced availability of dopamine in specific periods exacerbates motor fluctuations [28,29].Hence, inhibition of A 2A AR lessens the worsening of PD symptoms, such as dyskinesia, making it a non-dopaminergic treatment, usually in combination with dopaminergic treatments [28].
Frequent drug-receptor interactions of known drugs, such as hydrophobic contact, hydrogen bonding (HB), and π-stacked interactions [45], were observed in the selected protein-ligand complexes shown in Figure 2. The ASN-emodin complex formed a salt bridge between LYS45 and the alkoxide of the o-dihydroxy aromatic moiety in emodin.Alkyl and π-alkyl interactions were also a contributor to the stability of this complex.For the MAOB-chrysophanol complex, an attractive charge between LYS296 and the alkoxide of the phenolic structure, as well as HB with SER59 and TYR60, caused stabilization of the formed complex.This complex is further stabilized with stacked π-π, π-alkyl, and alkyl interactions.The superior binding of the MAOB-emodin complex formed an attractive charge between LYS296 and the alkoxide and is further stabilized with π-sulfur, stacked π-π, π-alkyl, and alkyl interactions.At the same time, the MAOB-physcion complex formed HB between CYS172 and TYR435 with π-sulfur, stacked π-π, T-shaped π-π, π-alkyl, and alkyl interactions, suggesting excellent molecular orbital alignment.The COMT-emodin complex formed HB between ASP141 and ASN170 and an attractive charge further supported by an attractive charge, π-cation, π-sulfur, π-alkyl, and alkyl interactions.For the A 2A AR-emodin complex, salt bridges and attractive charges (ARG107, ARG206, LYS1060, ARG1008), as well as HB (ARG107, ARG199, ARG206, and GLU1005), were formed; however, there was an unfavorable donor-donor interaction (ARG1008) which could potentially destabilize the complex structure.Finally, the A 2A AR-rhein complex formed HB (ARG199, LEU202, ALA203, ARG206, and ARG1008), attractive charge (ARG107), and π-alkyl (LEU202 and ALA203).These interactions reveal considerable effects on the PD targets; however, this static simulation could not fully uncover the protein conformation changes that could affect protein-ligand interactions.Hence, further analysis with molecular dynamics simulations could uncover these details.

Density Functional Theory (DFT) Analysis
The energies of frontier molecular orbitals (FMO) can provide helpful quantum mechanical calculations describing compounds' chemical reactivity and stability.According to Koopmans' theorem, the energy of the highest occupied molecular orbital (εHOMO) can be approximated to the required first ionization energy, and the energy of the lowest unoccupied molecular orbital (εLUMO) approximates the electron affinity in the last occupied orbital.In contrast, the energy gap (Δε) relates to the reactivity of the compound [46][47][48][49].Global reactivity descriptors from approximated ionization energies and electron affinity revealed beneficial molecular properties.In a chemical process, the hardness,

Density Functional Theory (DFT) Analysis
The energies of frontier molecular orbitals (FMO) can provide helpful quantum mechanical calculations describing compounds' chemical reactivity and stability.According to Koopmans' theorem, the energy of the highest occupied molecular orbital (ε HOMO ) can be approximated to the required first ionization energy, and the energy of the lowest unoccupied molecular orbital (ε LUMO ) approximates the electron affinity in the last occupied orbital.In contrast, the energy gap (∆ε) relates to the reactivity of the compound [46][47][48][49].Global reactivity descriptors from approximated ionization energies and electron affinity revealed beneficial molecular properties.In a chemical process, the hardness, of molecules describes the resistance of electron cloud deformation to trivial perturbations.
In contrast, the chemical potential, shows the stability of the compounds through their electron transfer capability, where negative values indicate the compound for less propensity for decomposition [46,49,50].
The electrophilic index, is one of the significant descriptors, since it describes potential affinity to a receptor for drugreceptor interactions.This descriptor assesses the inclination of the compound to acquire additional charge and its resistance to electronic charge exchange with the surroundings, indicative of electrophilicity [50].A large energy band gap between the HOMO and LUMO is classified as hard molecules and is consequently less reactive.Conversely, soft molecules are more reactive for smaller energy band gaps [50,51].From this, excellent inhibitors are characterized by their heightened stability.
Compared to the minimized structures (in Table S7), all protein-ligand complexes except rhein in complex with A 2A AR (PDB ID: 3EML) exhibited an increased energy band gap, revealing increased reactivity.Emodin (2.4536 eV) with ASN (PDB ID: 1XQ8) increased from 2.2949 eV.At the same time, chrysophanol (2.6025 eV), emodin (2.3333 eV), and physcion (3.0351 eV) with MAOB (PDB ID: 2C65) increased from 2.5902, 2.1354, and 2.9551 eV.At the same time, emodin (2.3313 eV and 2.3155 eV) with COMT (PDB ID: 3BWM) and A 2A AR (PDB ID: 3EML) increased from 2.1354 eV.However, rhein (2.7222 eV) with A 2A AR (PDB ID: 3EML) indicated a decrease from 2.8997 eV.Compared to the standards, MAOB-physcion and A 2A AR-rhein exhibited better protein-ligand interaction than levodopa.Similarly, ASN-emodin and MAOB-chrysophanol have potentially more stable protein-ligand interactions.The chemical hardness of the ligands shown in Table 3 indicated less polarizability, where the mentioned are again superior in stability.Moreover, chemical potentials reveal chemical stability from decomposition and electron escaping tendency.All ligands were considered to have superior stability compared to the standards.The electrophilicity indices show rhein and physcion as least nucleophilic.

Predicted Inhibitory Activity from 3D QSAR Models
Quantitative structure-activity relationship models were established in this study for the prediction of the inhibitory activity of RP major metabolites on PD targets.Steric (from van der Waals potential of carbon atom probe) and electrostatic (from H + point charge probe) contributions were considered in the 3D grid space of aligned and minimized ligands with known half-maximal inhibitory concentrations.A dataset of 45 known ASN inhibitors, 140 known MAOB inhibitors, 105 known COMT inhibitors, and 120 known A 2A AR inhibitors was considered in the 3D QSAR model building.As shown in Table 4, all models are considered robust and reliable from set significance based on their internal validity.For a further assessment of their accuracy, cross-validation of the models reveals the most accuracy from MAOB (0.488) > A 2A AR (0.400) > COMT (0.303) > ASN (0.225).On the other hand, external validation reveals the most accuracy from MAOB (0.458) > COMT (0.433) > ASN (0.387) > A 2A AR (0.175).Figures S1-S4 indicate predicted and experimental pIC 50 values in the training and test datasets, including their residuals.All predicted values of IC 50 to each protein target are within the dataset range as shown in Table 5. Consistently, emodin (15.957 µM) and rhein (18.457 µM) had the best inhibition of ASN among all metabolites.Compared to epigallocatechin-3-gallate with fibrillogenic inhibitory activity (9.8 µM), the best RP metabolites had 1.63 and 1.88-fold less inhibitory activity [52].This inhibitory activity on ASN could decrease fibril formation and aggregation, thereby decreasing the toxicity of the oligomeric Lewy bodies to the neurons.Similarly, chrysophanol (0.088 µM) and emodin (0.128 µM) had similar outcomes as their binding affinity with MAOB, which are comparable with selegiline (0.0013 µM), a known MAOB inhibitor [53].For this, the inhibition of MAOB slows down dopamine degradation, which restores the motor functions of PD patients.Compared to entacapone (0.386 µM) for COMT inhibition [54] and istradefylline (5.25 µM) for A 2A AR antagonism [55], the comparable inhibition from aloe-emodin (0.330 µM) and good inhibition with A 2A AR of 1.503 µM was revealed, based on its structural features.This interaction could potentially decrease the metabolic activity of levodopa for its increased therapeutic effect and sustained dopamine levels.Also, the decreased blocking activity of A 2A receptors to dopamine D2 receptors could potentially restore the balance of the motor circuitry.Nevertheless, chrysophanol (0.816 µM) was shown to have good inhibitory activity with A 2A AR.Previous literature experimentally measured the inhibitory concentrations of the major metabolite on MAOB, where only emodin was observed to have appreciable inhibitory activity [56].However, these results are open to further validations, specifically from immunochemical assays, protein-ligand binding assays, mutagenesis studies, functional assays, and X-ray crystallographic analyses.

ADMET Properties
Predicted pharmacokinetics and toxicity profiles from structure-activity relationships were performed in ADMETboost and pkCSM, as summarized in Table 6.All major metabolites of RP appeared to satisfy Lipinski's rule of five.Only chrysophanol did not follow the Pfizer rule, which could potentially have a toxic profile.At the same time, both chrysophanol and physcion do not follow the GSK rule, which may provide subpar ADMET properties.Nevertheless, the physicochemical properties of these metabolites follow druglike soft rules related to their ease of metabolism without substantially producing bioactive metabolites or intermediates [57].For this reason, undesired or toxic side effects from metabolic products and drug-drug interactions are minimized.All compounds have desirable n-octanol/water partition coefficients; however, the metabolites are significantly hydrophobic compared to the standards.Chrysophanol and physcion had good Caco-2 permeability with more than −5.15 log cm/s compared to the standards.On the other hand, all major RP metabolites, including the standards, have excellent human intestinal absorption (with more than 70% probable), which indicates a decreased requirement for parenteral administration.Since traditional medicines commonly use the oral route of administration and aqueous solvent, the good n-octanol/water distribution (1-3 log mol/L) observed for all metabolites and standards and poor aqueous solubility (−4-0.5 log mol/L) could potentially lead to decreased potency of RP but with optimum lipophilicity preferred for medicines for optimal lipid bilayer penetration [59,60].Nevertheless, the oral bioavailability of the metabolites and standards are comparable with all compounds possessing a good plasma protein binding rate, since they have less than a 90% probability of combining completely with blood serum proteins.This result means that major RP metabolites could potentially have quick onset action, and they can reach the target site rapidly [61,62].Chrysophanol, physcion, and rhein may cross the blood-brain barrier lightly (between 30-70%) compared to other metabolites and better than levodopa, potentially not requiring advanced drug delivery techniques for therapeutic effect.For these reasons, the metabolites may have good absorption and distribution to affect the required targets appreciably.
The compounds were found to have mediocre CYP2C9 and CYP2D6 inhibitory activities.However, all compounds were revealed to have an excellent probability of CYP2D6 inhibition.The CYP2D6 enzyme has mediocre metabolic activity in these compounds.This enzyme is highly expressed in the SN region of the CNS and is discouraged as a main metabolic route for CNS drugs due to their polymorphic nature [63].While the half-life of the metabolites was superior to the standards, corresponding to a lengthier action, the hepatocyte clearance was faster than the standards, and the microsome clearance was slower.For this, a cumulative effect may be observed, potentially leading to accumulative ability, and their lipophilicity could have a toxic effect.At the same time, the metabolites reveal a toxicity profile akin to the standards.There is an observed small probability for hERG blocking potential, which could lead to cardiotoxicity.However, this probability could be assumed negligible since the RP metabolites were approximately like the standards.Further investigations into the interaction of the metabolites to hERG with immunochemical assays, molecular simulations, and electrophysiological studies could validate the predicted cardiotoxicity profile.However, mediocre mutagenicity and hepatotoxicity were observed for all metabolites.On the other hand, aloe-emodin, chrysophanol, and physcion were observed to have significantly higher tolerable doses than the other metabolites and the standards, and all metabolites had significantly low lethal doses for oral application.In fact, in vitro and in vivo investigations showed that excessive use of emodin manifested hepatotoxicity, renal toxicity, and reproductive toxicity with chrysophanol [64].For this reason, the dosage of these compounds requires careful consideration since accumulation could be lethal.
All in silico pharmacokinetics and toxicity property results were consistent with other tools like pkCSM and ADMETlab 2.0.These findings indicate that RP metabolites may have a favorable therapeutic effect as a neuroprotective agent.

Preparation of Ligands
This study employed the molecular modeling and visualization software BIOVIA Discovery Studio 2022 (Dassault Systèmes, Vélizy-Villacoublay, France) for molecular simulations.Two-dimensional (2D) structure-data files (SDFs) of the identified ligand molecules were collected from the PubChem database (chem.ncbi.nlm.nih.gov,accessed on 10 May 2023) and prepared with the Prepare Ligands protocol in Discovery Studio (see Table S1).Numerous ionization states were generated under pH 7.5 ± 1.0.A structure visualization of the ligands was generated from ChemDraw v.19.

Preparation of Proteins
Experimentally elucidated structures of protein targets available from the Protein Data Bank (PDB, www.rcsb.org,accessed on 10 May 2023) in May 2023 were collected in PDB file format.Afterwards, water molecules and irrelevant heteroatoms for the simulation were removed, and polar hydrogens were added to the protein structure.Minimization of the protein structure with co-crystal ligands with the CHARMm forcefield was performed to verify the accuracy of the results.The protein was prepared with the Prepare Protein protocol in Discovery Studio, which automatically repaired and protonated the protein structure under a pH of 7.4 and 0.145 M of ionic strength with 0.9 kcal/mol cutoff energy.Missing loops were built in the protocol using PDB SEQRES loop definition, limited to 20 residues, with loop refinement through CHARMm forcefield minimization.

Molecular Docking
The initial docking runs were performed using the LibDock algorithm in Discovery Studio to screen ligand forms.Numerous conformations were generated with conjugategradient and quasi-newton minimization with an energy threshold between isomers set to 20.0 kcal/mol.A maximum fraction of steric clashes of 10% and a final clustering radius of 0.5 Å were set for the docking.The fractions of acceptable solvent accessible surface area (SASA) cutoffs were set to 15% for apolar SASA and 5% for polar SASA with 18 grid points.The binding sites were determined from the literature and PDB site records (see Table S2) [8,[65][66][67][68][69]. Docking optimization was performed using a gridbased CDOCKER algorithm in Discovery Studio for the evaluation of the best ligand for each protein target.Random conformations were generated from 1000 dynamics steps at 1000 K target temperature or with consideration of electrostatic interactions.The docked conformations were refined with 800 maximum bad orientations and a vdW energy threshold of 300 kcal/mol.This complex was annealed for 2000 heating steps until 700 K and subsequently cooled for 5000 steps to 300 K. Final minimizations of the refined ligand pose were set to full potential.The docking analysis results were further evaluated through a comparison to standard drugs of Parkinson's disease, i.e., dopamine and levodopa [5,30].

3D QSAR Modeling and Activity Prediction
The QSAR model for each protein target was modeled after IC 50 values of known interacting compounds collected in the ChEMBL database ((ebi.ac.uk/chembl, accessed on 10 May 2023) and enumerated in Tables S3-S6.These compounds were prepared using the Prepare Ligands for QSAR protocol in Discovery Studio, where duplicate data were removed, compound structure and charge were repaired, and numerous ionization states were generated under pH 7.4.The compounds were subjected to CHARMm minimization with Momany-Rone ligand partial charge estimation until no energy gradient was reached.Outlier compounds were removed through software-predefined parameters (atomic logarithmic 1-octanol/water partition coefficient, molecular weight, number of H donors, H acceptors, rotatable bonds, rings, aromatic rings, and fractional polar surface area).Outliers were classified under Poisson distribution considering short-ranged parameters with a p-value of 0.01.Rigid molecular overlay for compound alignment was performed under consensus-based, 50% steric, and 50% electrostatic field fit.The IC 50 values of the compound data were converted into pIC 50 values.These data are separated by diverse split through the predefined parameters into 80% training set and others for the test set.Three-dimensional QSAR models based on shape (steric) and charge (electrostatic) complementarity were generated for each protein target with pIC 50 values as the activity variable, grid spacing of 1.5 Å, and 5-fold cross-validation.The 3D QSAR model correlates the relationship between the activity variable and molecular fields following Equation (4): Variables in this model are n, the number of descriptors which are each grid point; k, the model coefficient for the descriptor; and, V, the descriptor value at the grid point corresponding to the compound.These models are built with partial least squares (PLS) regression.

Density Functional Theory (DFT) Analysis
Molecular properties of ligand structures in the form of frontier molecular orbitals, band gap (∆ε), and electrostatic potential maps were evaluated with the DMol 3 package in Discovery Studio for the structure stability analysis.Reference ligands were optimized with CHARMm forcefield minimization.Single-point energy calculation of the reference and docked ligands was performed with Becke's three-parameter, Lee-Yang-Parr (B3LYP) exchange-correlation functional and aqueous environment using the COSMO continuum solvent models.A double-numeric quality basis set with polarization functions (DNP) and a self-consistent field (SCF) density convergence of 1.0 × 10 −6 was used for the analysis.

ADMET Properties Prediction
SMILES notation records of the prepared ligands were collected from the PubChem database (chem.ncbi.nlm.nih.gov,accessed on 10 May 2023).ADMETboost (ai-druglab.smu.edu/admet, accessed on 10 May 2023) was used for ADMET properties prediction for further analysis.This predictive model utilized MACCS, extended connectivity circular, Mol2Vec, and PubChem fingerprints, as well as Mordred and RDKit descriptors through the extreme gradient boosting machine learning modeling model [58].Parameters include blood-brain barrier (BBB) permeability, metabolism (inhibition and substrate activity on human cytochrome P450 family), half-life, human ether-a-go-go-related gene (hERG) for cardiotoxicity, human hepatotoxicity, and AMES test for mutagenicity.The human maximum tolerable dose was determined from pkCSM (biosig.lab.uq.edu.au/pkcsm,accessed on 10 May 2023) [70].Data visualization was generated from GraphPad Prism v.9.5.0.

Conclusions
With increased patient cases, limited treatment diversity, and the complex and unclear pathogenesis of PD, the current state of research is far beyond a firm and near-complete understanding of the disease and, consequently, remains complicated with regards drug development.Current efforts for PD treatment continue to be symptomatic treatment and neuroprotection to delay its progression.Traditional medicines have been an excellent source of diverse compounds, and because of the complex mechanisms involved in PD, a multiple-targeted approach could potentially uncover interesting mechanisms and interactions.This study initially addressed the interaction of major RP metabolites for neuroprotective activities to Parkinson's disease in silico, and emodin could be attributed as the main effect compound in the extracts.These metabolites inhibited Parkinson's disease targets on par with controls and known inhibitors from its binding activity, and predicted the inhibition activity from structure-activity relationships with good chemical stability.

Future Directions
While there is an interesting outlook on the multi-targeted effect of RP metabolites on key PD targets, experimental analyses are required to validate the computational results from this study.A more thorough computational analysis (such as network pharmacology and molecular dynamics simulations) could be initially performed for the prediction of the mechanism of action.In addition to previously mentioned validation studies, the activity of RP extracts or decoctions, including such extracts with in vitro and in vivo models of PD, may provide additional insights into its treatment.Gene expression studies and protein-ligand interaction experiments are recommended for further interpretation of their mechanism of action.Following this, structure isolation and comparative explorations of these extracts should provide the structural moiety and molecular interactions necessary for PD treatment.Before advanced investigation and clinical trials, these compounds should be subjected to experimental pharmacokinetics and toxicity studies.

Figure 1 .
Figure 1.Chemical structures of RP major metabolites.

Figure 1 .
Figure 1.Chemical structures of RP major metabolites.

Table 1 .
LibDock scores of RP major metabolites to PD targets.

Table 1 .
LibDock scores of RP major metabolites to PD targets.

Table 2 .
CDOCKER energy (in kcal/mol) of RP major metabolites to PD targets.

Table 3 .
Electronic properties of standards and docked conformations of ligands (in eV).
1Minimized structures of standard drugs.

Table 5 .
Predicted IC 50 (in µM) of standards and RP metabolites on PD targets.

Table 6 .
ADMET properties in silico with ADMETboost and pkCSM webserver.