In silico studies of natural product-like caffeine derivatives as potential MAO-B inhibitors/AA2AR antagonists for the treatment of Parkinson's disease

Abstract Parkinson’s disease is considered the second most frequent neurodegenerative disease. It is described by the loss of dopaminergic neurons in the mid-brain. For many decades, L-DOPA has been considered as the gold standard for treating Parkinson’s disease motor symptoms, however, due to the decrease of efficacy, in the long run, there is an urgent need for novel antiparkinsonian drugs. Caffeine derivatives have been reported several times for their neuroprotective properties and dual blockade of monoamine oxidase (MAO) and adenosine A2A receptors (AA2AR). Natural products are currently attracting more focus due to structural diversity and safety in contrast to synthetic drugs. In the present work, computational studies were conducted on natural product-like caffeine derivatives to search for novel potent candidates acting as dual MAO-B inhibitors/AA2AR antagonists for Parkinson’s disease. Our findings revealed two natural products among the top hits: CNP0202316 and CNP0365210 fulfill the requirements of drugs acting on the brain. The selected lead compounds were further studied using molecular dynamics simulation to assess their stability with MAO-B. Current findings might shift the interest towards natural-based compounds and could be exploited to further optimize caffeine derivatives into a successful dual-target-directed drug for managing and halting the neuronal damage in Parkinson’s disease patients.


Introduction
Neurodegenerative diseases and brain-associated diseases are major concerns among aging populations across the world [1]. Neurodegenerative diseases such as Parkinson's and Alzheimer's diseases have a multifactorial nature that is characterized by the progressive loss of neurons in the brain [2]. Parkinson's disease (PD) is defined especially by the progressive loss of dopaminergic neurons in the substantia nigra pars compacta (SNpc) of the midbrain [3]. More than six million people in the world are affected today with a prevalence of 150 in every 100,000 people which is further increasing with age and affects 1% of the population over 60 years [4]. Current pharmaceutical treatments for PD include levodopa or levodopa plus dopa-decarboxylase inhibitors, dopamine agonists, and catechol-O-methyl transferase (COMT)/monoamine oxidase B (MAO-B) inhibitors [5]. Recently, other non-dopaminergic drugs have shown promising efficacy to relieve PD symptoms such as adenosine A 2A receptor (AA 2A R) antagonists [6].
Monoamine Oxidase (MAO) (EC 1.4.3.4) belongs to a family of flavin adenine dinucleotide (FAD)dependant enzymes that are expressed in the outer mitochondrial membrane of neuronal cells. The MAO enzymes are responsible for the oxidative deamination of monoamine neurotransmitters such as dopamine, adrenaline, and noradrenaline in the central nervous system (CNS) [5,7]. The MAO enzymes exist in two isoforms, MAO-A and MAO-B that share sequence similarities of 70% but differ in tissue distribution, substrate, and inhibitor preferences [5]. The development of the first MAO inhibitors was abandoned due to side effects related to the metabolism of tyramine, which causes a cardiovascular crisis [8]. However, a new class of selective MAO-B inhibitors has been proven to be efficient in treating PD symptoms. It was also shown that this new class of selective MAO-B inhibitors is devoid of tyramine-related side effects. Furthermore, the selective MAO-B inhibitors may act as neuroprotective agents by limiting the release of free radical species and hence may decrease the progression of the disease [5,9].
During aging, the expression of MAO-B increases in the brain and is connected with an enhanced dopamine metabolism which results in an increased reactive oxygen species (ROS) production such as hydrogen peroxide (H 2 O 2 ) inducing oxidative damage and apoptotic signaling events [11].
Previously approved MAO-B inhibitors are selegiline and rasagiline which irreversibly inhibit MAO-B with an IC 50 value of 6.8 and 14 nM respectively [12]. The latest approved MAO-B inhibitor is safinamide which reversibly inhibits MAO-B with an IC 50 value of 450 nM [13]. Istradefylline, a caffeine-based inhibitor that was approved in Japan in 2013 and also approved for medical use in the United States in 2019 acts as a dual inhibitor of MAO-B and AA 2A R [14,15]. However, istradefylline was found to be a weak inhibitor of MAO-B (IC 50 = 28 μM) which encourages further research on new substitutions to the caffeine core [16]. The chemical structures of MAO-B inhibitors are shown in Figure 1.
The crystal structure of MAO-A (PDB ID: 2Z5Y) has a monopartite substrate cavity of ∼550 Å 3 volume while the crystal structure of MAO-B contains a bipartite cavity structure with an entrance cavity of ∼290 Å 3 and a substrate cavity of ∼400 Å 3 [17]. ILE-199 and TYR-326 separate these two cavities in MAO-B serving as "gating" residues and a structural determinant for substrate and inhibitor recognition by MAO-B [18,19]. The superposition of MAO-A and MAO-B and their active site residues are shown in Figure 2.
The structural study revealed that MAO-B (PDB ID: 2V5Z) is formed by two monomers consisting of a globular domain anchored to the membrane through a C-terminal helix [20]. MAO-B active site residues that share similarities to MAO-A active site are TYR-60, LEU-164, PHE-168, GLN-206, ILE-198, ILE-316, PHE-343,  TYR-398, and TYR-435. Meanwhile, the amino acids that are specific to MAO-B are located in the hydrophobic pocket which is formed by LEU-171, CYS-172, ILE-199, TYR-326 [13].
There is a great deal of literature supporting the use and efficacy of natural products (NP) in PD such as flavonoids, xanthones, phenolic derivatives, alkaloids, and caffeine [21,22]. These natural resources and their derivatives have been reported for their potential to selectively inhibit MAO-B and may offer a safer alternative compared to conventional drugs [23]. Furthermore, caffeine has been used in several studies as a scaffold for the design of dual MAO inhibitors/AA 2A R antagonists. Pretorius et al. synthesized a series of C-8 substituted caffeinyl analogues and it was found that the compound bearing a 4-phenylbutadien moiety is the most potent candidate for MAO-B and AA 2A R [24]. On the other hand, Azam et al. explored numerous caffeine derivatives from the literature bearing multiple substitutions through molecular docking and structure-activity relationship studies, it was found that the placement of hydrophobic moieties at C8 is essential for both MAO-B inhibition and AA 2A R antagonism, whereas replacements occurring at C1 and C3 are optimal for AA 2A R but not detrimental for MAO-B [25]. Although research on caffeine is underway for decades, its naturally occurring derivatives are yet to be investigated in detail [26].
NPs and NP-based compounds are an ideal choice for scientists and researchers due to the broad-spectrum activity of NPs with their minimal or no toxic effect on human health [27]. The literature has indicated that caffeine among other NPs is a potent compound that has neuroprotective properties [28]. Considering the link between neurodegeneration and oxidative stress due to the mitochondrial imbalance and the accumulation of reactive oxygen species (ROS), MAO-B was and still is, considered a valid therapeutic target for slowing down the progression of Parkinson's disease.
In the present study, a substructure search was conducted on natural products databases to retrieve caffeine-containing natural products since it is known for its neuroprotective properties and its potency to act as an antagonist of AA 2A R, a validated target for PD [29]. Structure-based virtual screening was employed to evaluate the affinity of the selected natural compounds towards MAO-B and AA 2A R. ADMET properties were evaluated using in silico methods. Finally, molecular dynamics simulations were performed to study the interactions and the stability between the selected compounds and MAO-B over the simulation time.

Data sources
To retrieve all the available natural compounds based on the caffeine scaffold, we used the COCONUT database (https://coconut.naturalproducts.net/); the largest open-source natural products (NP) database to date containing more than 400 000 unique NP from over 50 sources [30]. The search was conducted using the Ulmann algorithm for the substructure search with the caffeine scaffold as a pharmacophore [31]. The search results revealed 144 caffeine-containing natural products. These compounds were downloaded in SDF format for further analysis.
Residues with missing atoms were fixed using the CHARMM-GUI web server [32]. Water molecules were removed since they are not involved in the ligand binding. Since MAO-B is expressed as a dimer, only one chain was kept along with the FAD cofactor for the molecular modeling studies to ease the computational cost [33]. Finally, polar hydrogens and Kollman charges were added using AutoDockTools 1.5.6 [34]. The grid box was placed near the FAD cofactor with a spacing of 1 Å. Grid dimensions were chosen large enough (24 × 24 × 24 Å in x, y, and z directions, respectively) to fit all the residues forming both cavities of the active site in the protein. The grid box was positioned in a way to cover the entire binding site and to allow larger molecules to dock properly: 53 × 155 × 27 Å for MAO-B and −21.6 × 6.1 × 17.5 for AA 2A R in x, y, and z directions, respectively. Lastly, the generated coordinates for the grid box were saved in a text file.

Preparation of ligands
The selected caffeine derivatives were split into multiple files, with each file containing a single ligand. The 3D conformations were generated for all the compounds, geometrical optimization was performed using Merck molecular force field (MMFF94) implemented in the Open Babel chemical toolbox [35]. The minimized ligands were then prepared for the molecular docking study using the prepare_ligand4.py package of AutoDockTools 1.5.6. Partial charges, atomic types, and polar hydrogens were added to all compounds and then converted to PDBQT format.

Structure-based virtual screening workflow
Structure-based virtual screening was performed using a Perl script for the automated execution of AutoDock Vina 1.1.2 [36]. The proposed methodology is detailed in Figure 3, a text file containing all the names of the prepared ligands was created to serve as a single input file for the docking screens. To facilitate the analysis of the virtual screening results, all the generated output log files were concatenated into a single output text file. All procedures were performed in respect of good practices using state-of-the-art virtual screening approaches for natural products bioprospecting as shown in Figure 4. The standard virtual screening consists mainly of target identification, selection of the chemical library, molecular docking studies, ADMET evaluation, molecular dynamics simulations, and finally experimental validation of the lead compounds [37].

Visualization and analysis
The ligands were ranked by their binding affinities. The compounds displaying a binding score of −8.0 kcal/mol or less were subject to further analysis. The conformations of the selected compounds were visualized using UCSF Chimera visualization software and superposed to the ligand of reference [38]. Discovery Studio Visualizer program was used to identify hydrogen bonds and other hydrophobic interactions [39].

In silico ADMET prediction
The profiling of compound pharmacokinetics is very essential in drug discovery. As of today, many online tools are available to predict the ADMET profiles of drugs based on their chemical properties [40]. In silico ADMET profiling can be useful to speed up the drug development process by limiting the number of compounds for experimental testing. In this study, physicochemical properties and pharmacokinetic parameters were evaluated using SwissADME online calculation toolkit (http://www.swissadme.ch/) [41]. Lipinski's rule of five was taken into account to assess the ability of the compounds to be active for orally administrated drugs [42]. Other parameters such as water solubility, gastrointestinal absorption, and blood-brain barrier permeability were predicted. Pain-assay interference compounds (PAINS) are chemical molecules that often give false-positive results in high-throughput  screens due to the presence of several disruptive functional groups that interact nonspecifically with various biological targets rather than selectively affecting the therapeutic target of interest [43]. Hepatotoxicity was predicted using the ProTox-II web server (https://tox-new.charite.de/protox_II/) [44].

Molecular dynamics simulations
AMBER16 software [45] was employed to conduct molecular dynamics simulations on the most potent compounds in complex with MAO-B. AMBER force field 14SB [46] and the general AMBER force field (GAFF2) [47] were used to parametrize the protein and the identified inhibitors, respectively. The TIP3P water model with a margin of 15.0 Å (1.5 nm) in each direction from the solute was used to construct a water-solvated cubic box. The specifics of the used MD simulations are elucidated in Ref. [48][49][50][51][52][53][54][55]. In synopsis, energy minimization was initially used on the investigated inhibitors in complex with MAO-B for 5000 steps using the combined steepest and conjugate gradient algorithms. Thereafter, the minimized systems were progressively heated from 0 k to 300 k over 50 ps. The complexes were equilibrated to a free simulation for 1000 ps. Ultimately, a production run for 100 ns was subsequently carried out utilizing an NPT ensemble at 300 K with 1.0 atm pressure. All the periodic boundary PME (Particle Mesh Ewald) simulations were conducted using the "pmemd.cuda" implementation in AMBER for GPU-accelerated simulations on the CompChem hybrid GPU/CPU cluster.

MM-GBSA binding energy
The molecular mechanics-generalized Born surface area (MM-GBSA) approach was applied to estimate the binding free energies (ΔG binding ) of the investigated inhibitors in complex with MAO-B [56]. Thus, the total binding free energies were evaluated according to IGB value of 2. For each system, the binding free energy calculations were executed for 10,000 snapshots recorded throughout 100 ns MD simulations. For each snapshot, the MM-GBSA (ΔG binding ) binding energy was calculated by the standard formula:

Results
Developing efficient therapies against neurodegenerative diseases such as Parkinson's disease remains a great challenge. The use of natural products has been known for a long time to offer great promise and they're often a safer alternative compared to synthetized drugs. Currently, in silico studies are providing much-needed preliminary data about potential drugs, which can be a great help in conducting additional in vitro and in vivo studies [57].

Validation of molecular docking accuracy
Molecular docking protocol was first validated by cross-docking the co-crystallized ligands using the same parameters that were applied for the studied compounds against different crystallographic structures (PDB ID: 2V61 for MAO-B and PDB ID: 5IU4 for AA 2A R) [13,58]. Computed root-mean-square deviation (RMSD) was calculated by mean of superposition; the obtained values are below 2 Å which indicates a good quality of the docking program ( Figure 5).
Additionally, molecular docking accuracy was further validated using two datasets of 10 caffeine derivatives with reported half-maximal inhibitory concentrations (IC 50 ) for MAO-B and dissociation constants (Ki) for AA 2A R respectively [24,59]. A good correlation was established between the docking results and the experimental values for MAO-B and AA 2A R, which confirms the reliability of the molecular docking approach to study the natural product-like caffeine derivatives with MAO-B. The correlation charts, correlation coefficients, and slopes are shown in Figure 6.

Natural product-like caffeine derivatives screening
In the present study, we screened 144 natural product-like caffeine derivatives against MAO-B using structurebased virtual screening. The compounds were ranked by their binding affinities (kcal/mol). The highestranking molecules displaying a docking score of −8.0 kcal/mol or less were further analyzed based on their interactions with the MAO-B active site cavity. Molecular docking results and protein-ligand interactions are shown in Table 1.

In silico ADMET prediction results
In silico pharmacokinetics, toxicity, and drug-likeness prediction results are shown in Table 2. All compounds were predicted as either soluble or highly soluble. Furthermore, most of the molecules are showing high gastro-intestinal absorption which indicates a good oral bioavailability. However, the blood-brain barrier permeability parameter revealed only two compounds besides safinamide that may readily cross the bloodbrain barrier and act on the central nervous system. Moreover, all compounds, excluding safinamide, were identified as non-inhibitors of CYP2D6, which is particularly necessary for drugs acting on the brain since the expression of CYP2D6 is higher in the brain and is involved in metabolizing endogenous neural compounds that suggest its neuroprotective effects [60]. Moreover, the inhibition of CYP enzymes can decrease drug efficacy leading to therapeutic failure or increased drug side effects and toxicity [61][62][63]. Organ toxicity predicted using the ProTox-II webserver revealed that all the compounds are safe for the liver and do not disrupt its normal function. Physicochemical properties profiling of the selected compounds revealed that all the compounds are drug-like according to Lipinski's rule of five. Finally, Pain-assay interference compounds (PAINS) alerts calculations indicated that all the compounds do not contain any disruptive functional groups except CNP0074614 displaying one PAINS alert due to the catecholamine group.

Interaction analysis of lead compounds with MAO-B and AA 2A R
According to the molecular docking and ADMET-based screening as summarized in Figure 7, two compounds were identified as potential drug candidates that possess the desired pharmacokinetics properties for drugs acting on the central nervous system: CNP0202316 and CNP0365210 superposed to the reference inhibitor, safinamide in complex with MAO-B are illustrated in Figure 8. The binding scores of these molecules were −10.1 and −9.7 kcal/mol for MAO-B respectively and are comparable to the reference inhibitor, safinamide Water solubility, insoluble < −10 < poorly < −6 < moderately < −4 < soluble < −2 < very < 0 < highly; GI absorption, gastrointestinal absorption; BBB, blood-brain barrier permeability; CYP2D6 inhibitor, Likeliness of a drug to act as inhibitor of cytochrome P450 CYP2D6; Hepatotoxicity, prediction of drug-induced liver injury; Lipinski violation, number of violations to the rule of five (log P o/w ≤ 5; MW ≤ 500 g/mol; HBA ≤ 10; HBD ≤ 5; RB ≤ 10); PAINS alert, number of disruptive functional groups shared by many PAINS (Pan-assay interference compounds).

Figure 7:
Step-wise structure and ADMET-based screening of the selected natural product-like caffeine derivatives.
which displayed a binding score of −9.9 kcal/mol. However, molecular docking of istradefylline revealed a low binding affinity (−9.3 kcal/mol) compared to the aforementioned compounds. Safinamide was found to interact through hydrogen bonding with GLN-206 which is known to be a hydrogen bond acceptor for the majority of MAO-B inhibitors and is by the literature [13]. Structural analysis shows that both compounds share a phenyl ring linked to the caffeine scaffold with a pentane group. The phenyl ring seems to be favorable for the stability of the ligands within the MAO-B active site by establishing numerous hydrophobic interactions with the nearby residues of the entrance cavity. Meanwhile, the caffeine scaffold is directed towards the FAD cofactor and interacts with CYS-172 in CNP0202316 and TYR-435 of the aromatic cage in CNP0365210 through hydrogen bonding.
It has been shown in previous studies that hydrophobic interactions through the phenyl ring are vital for establishing MAO-B binding and are more favorable than all other interactions such as hydrogen or halogen bonds [64]. However, CNP0202316 where the phenyl ring is placed at C7 of the xanthine core seems to be more stable than the second compound implying that the presence of hydrophilic interaction with CYS-172 contributes more to the anchoring and the stability of this compound in the active site cavity of MAO-B.
Alternatively, the selected compounds were analyzed through molecular docking studies with AA 2A R to assess their binding affinities with the AA 2A R active site and compare their binding conformations to istradefylline. Molecular docking results and protein-ligand interactions are shown in Table 3.
The binding poses of the selected compounds and istradefylline were superposed to the co-crystallized structure of caffeine ( Figure 9). The structural analysis indicates that both compounds were able to interact with a key polar residue, ASN-253 through the formation of hydrophilic hydrogen bonds, similarly to istradefylline and other potent AA 2A R antagonists [65]. A large network of hydrophobic interactions was also  observed, where key residues namely PHE-168 and GLU-169 were found to make a consistent appearance. Since all compounds share the same core that characterizes the caffeine molecule, this would entail a somewhat similar disposition inside the binding pocket. This holds especially for CNP0202316 where the xanthine core was found to be positioned similarly to istradefylline. Moreover, the propylbenzene moiety at position C7 in CNP0202316 might be more favorable to the hydrophobic pocket of the receptor suggesting its affinity potential which may be on par or better than istradefylline.

Molecular dynamics simulations and binding energy calculations
Molecular dynamics (MD) simulations were applied to probe the stability of the selected ligand-protein complexes, structural specifics, conformational flexibilities, and realize reliable inhibitor-enzyme binding affinities [66,67]. The most promising compounds in complex with the MAO-B enzyme were further inspected via MD throughout 100 ns simulation time. According to the gathered inhibitor-enzyme snapshots over the production run of 100 ns, the MM-GBSA approach was used to calculate the binding free energies (ΔG binding ) and are illustrated in Figure 10. From the data in Figure

Post-MD simulations analysis
Molecular docking calculations, and MD simulations combined with MM-GBSA binding energy calculations, unveiled the most potent molecules as potential MAO-B inhibitors. The MD-based analysis could be required  to demonstrate structural and energetic stabilities for the scrutinized inhibitors in complex with MAO-B. The structural and energetical analysis included root-mean-square deviation (RMSD), Center-of-Mass (CoM) distance, and binding energy per-frame.

Root-mean-square deviation (RMSD)
The root-mean-square deviation (RMSD) values of the backbone atoms within the whole complex throughout the simulation time were estimated to monitor the structural stability of the CNP0202316, CNP0365210, and safinamide in complex with MAO-B. The RMSD of the backbone atoms as a function of time following the initial structure of the three investigated systems is displayed in Figure 12. The platform in RMSD curves emphasizes that all three inspected systems attain an equilibrium within 1000-10,000 ps throughout MD simulations, exposing that the three investigated systems are converged over the simulation window. These findings suggest that all the compounds are tightly bound and not influenced by the topology of the protein. This is especially true for CNP0202316, which displayed very slight deviations somewhat similar to safinamide suggesting its high stability compared to CNP0365210.

Center-of-mass distance
To get a more in-depth insight into the stability of the selected compounds throughout the MD simulation time, center-of-mass (CoM) distances were evaluated ( Figure 13). The most interesting aspect of this graph is  that CoM distances were consistent for the CNP0202316-and CNP0365210 in complex with MAO-B compared to safinamide-MAO-B complex, with average values of 5.7, 5.8, and 6.1 Å, respectively. The most obvious finding to emerge from this analysis is that CNP0202316 and CNP0365210 bound more tightly to the MAO-B complex than the reference inhibitor, safinamide.

Binding energy per frame
The comprehensive structural stability of CNP0202316, CNP0365210, and safinamide complexed with MAO-B was evaluated throughout a 100 ns MD simulation via inspecting the correlation between the binding energy per frame and time ( Figure 14). Overall stabilities for CNP0202316, CNP0365210, and safinamide were noticed with average binding energies (ΔG binding ) of −36.6, −34.5, and −37.9 kcal/mol, respectively. Based on this analysis, all investigated complexes preserved their stability over the 100 ns MD simulations.

Discussion
The results obtained from the molecular docking study with MAO-B were further analyzed through molecular dynamics simulations and binding free energy calculations. MM-GBSA binding energies further confirmed that van der Waals (ΔE vdw ) energy is a considerable contributor to the stability of MAO-B inhibitors. This finding confirms that hydrophobic interactions outweigh any other interactions in terms of MAO-B inhibition, and is also supported by the literature [65]. Post-MD simulations analysis confirmed the stability of the three compounds and revealed that all the complexes achieve equilibrium within 1-10 ns throughout the simulation time. These results suggest that the proposed natural products may be on par or better than the inhibitor of reference, safinamide regarding MAO-B inhibition. The key residues involved in MAO-B inhibition were found to be LEU-171, ILE-199, TYR-326, and GLN-206, these residues contribute the most to the stability of the inhibitors when bound to MAO-B. LEU-171, ILE-199, and TYR-326 are specific residues to MAO-B isoform, they are located in the entrance cavity and play a role in substrate and inhibitor specificity [68,69]. Meanwhile, GLN-206 is recognized as a hydrogen bond acceptor for most MAO-B inhibitors and is responsible for their stability in the substrate cavity [17]. Moreover, the identified compounds may confer neuroprotective properties linked to the xanthine core as reported in the literature [70]. For MAO-B, the orientation of the phenyl ring linked to the xanthine core was found to be similar to safinamide especially in CNP0202316, and is favorable to the entrance hydrophobic cavity. On the other hand, this molecule adopted a similar conformation to istradefylline when bound to AA 2A R. The propylbenzene moiety attached to the imidazole of the caffeine in CNP0202316 at position C7 is more favorable to the hydrophobic pocket of AA 2A R, whereas the oxygen atom of the xanthine core maintained a hydrogen bond with ASN-253 that is deemed crucial to the binding of most AA 2A R antagonists [71]. Thus, the identified compounds might offer a dual-target activity in the context of a polypharmacological approach and might represent a more efficient alternative for treating and slowing down the neuronal damage in Parkinson's disease patients.

Conclusions
The present study aimed to find novel natural product-like caffeine derivatives as potential dual MAO-B inhibitors/AA 2A R antagonists. Structure-based virtual screening and ADMET analysis revealed two natural products that fulfill the requirements for drugs acting on the brain. The selected compounds in complex with MAO-B were subject to molecular dynamics simulations to assess their stability over the simulation time along with the inhibitor of reference, safinamide. Our findings show that the presence of the phenyl ring in the selected compounds is crucial for the ligands to fill the long-shaped cavity of the MAO-B active site and is a major contributor to various van der Waals interactions responsible for the stability and the tight-binding of these compounds to MAO-B. Similarly, the propylbenzene moiety was found to be more favorable for the hydrophobic pocket of AA 2A R especially when linked at position C7 of the xanthine core which allowed the caffeine core to adopt a similar conformation to istradefylline suggesting the dual-target properties of the identified natural products. In conclusion, the structure-based virtual screening helped provide valuable insight on the studied natural product-like caffeine derivatives and our findings may attract more focus for the development of novel antiparkinsonian drugs with dual-targeting properties. However, in vitro experiments such as bioactivity assays for MAO-B and AA 2A R, membrane permeability and cell viability assays remains necessary to further validate these findings.