Academia.eduAcademia.edu
pharmaceuticals Article Exploring Toxins for Hunting SARS-CoV-2 Main Protease Inhibitors: Molecular Docking, Molecular Dynamics, Pharmacokinetic Properties, and Reactome Study Mahmoud A. A. Ibrahim 1, * , Alaa H. M. Abdelrahman 1 , Laila A. Jaragh-Alhadad 2,3, *, Mohamed A. M. Atia 4 , Othman R. Alzahrani 5 , Muhammad Naeem Ahmed 6 , Moustafa Sherief Moustafa 2 , Mahmoud E. S. Soliman 7 , Ahmed M. Shawky 8 , Paul W. Paré 9 , Mohamed-Elamir F. Hegazy 10 and Peter A. Sidhom 11 1 2 3 4   Citation: Ibrahim, M.A.A.; Abdelrahman, A.H.M.; Jaragh-Alhadad, L.A.; Atia, M.A.M.; Alzahrani, O.R.; Ahmed, M.N.; Moustafa, M.S.; Soliman, M.E.S.; Shawky, A.M.; Paré, P.W.; et al. 5 6 7 8 9 Exploring Toxins for Hunting SARS-CoV-2 Main Protease 10 Inhibitors: Molecular Docking, Molecular Dynamics, 11 Pharmacokinetic Properties, and Reactome Study. Pharmaceuticals 2022, 15, 153. https://doi.org/10.3390/ ph15020153 Academic Editor: Paweł Kafarski Received: 4 January 2022 Accepted: 23 January 2022 Published: 27 January 2022 Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations. Copyright: © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and * Computational Chemistry Laboratory, Chemistry Department, Faculty of Science, Minia University, Minia 61519, Egypt; a.abdelrahman@compchem.net Department of Chemistry, Faculty of Science, Kuwait University, Kuwait City 13060, Kuwait; mostafa_msm@hotmail.com Cardiovascular and Metabolic Sciences Department, Lerner Research Institute, Cleveland Clinic, Cleveland, OH 44195, USA Molecular Genetics and Genome Mapping Laboratory, Genome Mapping Department, Agricultural Genetic Engineering Research Institute (AGERI), Agricultural Research Center (ARC), Giza 12619, Egypt; matia@ageri.sci.eg Department of Biology, Faculty of Science, University of Tabuk, Tabuk 71491, Saudi Arabia; o-alzahrani@ut.edu.sa Department of Chemistry, The University of Azad Jammu and Kashmir, Muzaffarabad 13100, Pakistan; drnaeem@ajku.edu.pk Molecular Modelling and Drug Design Research Group, School of Health Sciences, University of KwaZulu-Natal, Westville, Durban 4000, South Africa; soliman@ukzn.ac.za Science and Technology Unit (STU), Umm Al-Qura University, Makkah 21955, Saudi Arabia; amesmail@uqu.edu.sa Department of Chemistry & Biochemistry, Texas Tech University, Lubbock, TX 79409, USA; paul.pare@ttu.edu Chemistry of Medicinal Plants Department, National Research Centre, 33 El-Bohouth St., Dokki, Giza 12622, Egypt; elamir77@live.com Department of Pharmaceutical Chemistry, Faculty of Pharmacy, Tanta University, Tanta 31527, Egypt; peter.ayoub@pharm.tanta.edu.eg Correspondence: m.ibrahim@compchem.net (M.A.A.I.); laila.alhadad@ku.edu.kw (L.A.J.-A.) Abstract: The main protease (Mpro ) is a potential druggable target in SARS-CoV-2 replication. Herein, an in silico study was conducted to mine for Mpro inhibitors from toxin sources. A toxin and toxin-target database (T3DB) was virtually screened for inhibitor activity towards the Mpro enzyme utilizing molecular docking calculations. Promising toxins were subsequently characterized using a combination of molecular dynamics (MD) simulations and molecular mechanics-generalized Born surface area (MM-GBSA) binding energy estimations. According to the MM-GBSA binding energies over 200 ns MD simulations, three toxins—namely philanthotoxin (T3D2489), azaspiracid (T3D2672), and taziprinone (T3D2378)—demonstrated higher binding affinities against SARS-CoV-2 Mpro than the co-crystalized inhibitor XF7 with MM-GBSA binding energies of −58.9, −55.9, −50.1, and −43.7 kcal/mol, respectively. The molecular network analyses showed that philanthotoxin provides a ligand lead using the STRING database, which includes the biochemical top 20 signaling genes CTSB, CTSL, and CTSK. Ultimately, pathway enrichment analysis (PEA) and Reactome mining results revealed that philanthotoxin could prevent severe lung injury in COVID-19 patients through the remodeling of interleukins (IL-4 and IL-13) and the matrix metalloproteinases (MMPs). These findings have identified that philanthotoxin—a venom of the Egyptian solitary wasp—holds promise as a potential Mpro inhibitor and warrants further in vitro/in vivo validation. conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/). Pharmaceuticals 2022, 15, 153. https://doi.org/10.3390/ph15020153 https://www.mdpi.com/journal/pharmaceuticals Pharmaceuticals 2022, 15, 153 2 of 22 Keywords: toxins; SARS-CoV-2 Mpro ; in silico screening; molecular docking calculations; molecular dynamics (MD) simulations; reactome 1. Introduction The causative factor in COVID-19 infection is Severe Acute Respiratory Syndrome Coronavirus-2 (SARS-CoV-2), a novel β-coronavirus of the positive-stranded RNA virus that results in gastrointestinal, respiratory, and neurological symptoms in humans [1,2]. From December 2019, a gigantic economic epidemic has been disseminated globally because of COVID-19 disease [3,4]. As of 29 December 2021, more than 281 million confirmed cases and over 5.4 million international deaths had been reported [5]. A small number of vaccines have currently been approved under emergency use authorization [6]. Notwithstanding the weak vaccination rate, the deficiency of specific therapies, and the development of numerous viral variants, the pandemic goes on to distribute quickly and intricately. As a consequence, more outstanding efforts are required to discover safe and potent drugs against SARS-CoV-2. SARS-CoV-2 main protease (Mpro ) is a crucial enzyme for viral gene replication, expression, and transcription [7–9]. Therefore, inhibition of the viral Mpro enzyme is a putative strategy towards SARS-CoV-2 antiviral drug development. Several in silico and experimental attempts have been made to repurpose approved drugs as prospective curative candidates for the remediation of COVID-19 [10–13]. Further, a combination of virtual screening and molecular dynamics (MD) simulations of chemical libraries towards SARS-CoV-2 targets has been executed [14–17]. Sundry small compounds have been identified during and after the first and second coronavirus prevalence waves. Among these, many edible plant-derived natural products and their related synthetic derivatives have attracted considerable interest as prospective COVID-19 drug candidates. Furthermore, natural products from marine species have recently demonstrated substantial antiviral characteristics [18–21]. Very recently, an emergency use authorization of PAXLOVID (PF07321332), a covalent Mpro inhibitor with submicromolar activity developed by Pfizer, has been granted by the U.S. Food and Drug Administration (FDA) for treatment of patients with mild-to-moderate COVID-19 [22]. Toxin or toxin subunits have been used as therapeutic agents to treat an enormous number of diseases when they are not capable of causing damage or death to humanity [10,23]. Among the approved drugs, eleven drugs were recognized as toxins, such as exanta, ziconotide, exenatide, and lixisenatide [24]. Toxin and Toxin-Target Database (T3DB) database is an unparalleled bioinformatics resource that collects overall information about popular or omnipresent toxins and their toxin-targets into a single electronic storehouse [25]. The database includes more than 2900 small compounds and peptide toxins, over 33,000 toxintarget associations, and 1300 toxin-targets [25]. Exploring toxins to hunt potential inhibitors towards SARS-CoV-2 Mpro has not been conducted. So, in the current study, the T3DB database was virtually screened as Mpro specific drug candidates. Based on the predicted docking scores, the most potent toxins were submitted to molecular dynamics (MD) simulations combined with binding energy calculations using the molecular mechanics-generalized Born surface area (MM-GBSA) approach. Pathway enrichment analysis (PEA) and Reactome mining was performed to dissect biological aspects of the inhibitor hits on drug-target interactions with an interactive layout [16,26]. Such in silico screening can provide worthy insights with respect to the appropriateness of the obtained hits as future development of prospective clinical candidates. 2. Results and Discussion In searching for small compounds to prohibit viral replication and transcription, in silico techniques were utilized to explore a chemical library containing more than Pharmaceuticals 2022, 15, 153 3 of 22 3678 toxins as potential SARS-CoV-2 Mpro inhibitors. The employed protocol was first validated based on available experimental data. 2.1. Validation of In Silico Protocol The performance of the utilized in silico protocol to predict the binding mode of SARS-CoV-2 Mpro inhibitor was evaluated. The co-crystallized (5S)-5-(3-{3-chloro-5-[(2chlorophenyl)methoxy]phenyl2-oxo[2H-[1,3′ -bipyridine]]-5-yl)pyrimidine-2,4(3H,5H)dione (XF7) inhibitor was redocked against the SARS-CoV-2 Mpro , as well as the anticipated binding mode was compared to the experimental binding mode (PDB code: 7L13) [27]. As shown in Figure 1, the predicted binding mode was very similar to the resolved experimental binding mode with an RMSD of 0.20 Å and a binding affinity of −9.5 kcal/mol. This data comparison revealed the superior performance of AutoDock4.2.6 software in anticipating the experimental binding mode of Mpro inhibitors. The robust binding of XF7 with Mpro is attributed to the NH and two CO groups of pyrimidine-2,4(1H,3H)-dione ring to form hydrogen bonds with a backbone CO, NH group of THR26 and the backbone NH of GLY143 with bond lengths of 2.49, 2.28, and 2.07 Å, respectively (Figure 1). Besides, a nitrogen atom of pyridine rings interacts with the backbone CO group of SER144 and the imidazole ring of HIS163 with bond lengths of 3.29 and 1.91 Å, respectively (Figure 1). At the same time, the carbonyl group of pyridin-2(1H)-one ring exhibits a hydrogen bond with the backbone NH group of GLU166 with a bond length of 1.81 Å (Figure 1). Therefore, the docking protocol confirms the outperformance of this approach in identifying potent inhibitors as prospective SARS-CoV-2 Mpro inhibitors. 2.2. T3DB Database Virtual Screening To identify Mpro inhibitors from toxins, AutoDock4.2.6 software was employed to virtually screen the T3DB database. At the outset, the T3DB database was filtered towards Mpro with conventional docking parameters. According to the portended binding affinity, 200 toxins demonstrated docking scores less than −8.0 kcal/mol towards Mpro . Therefore, those 200 toxins were submitted to more elaborate molecular docking calculations with costly docking parameters. The evaluated docking scores for the top 200 hits are summarized in Table S1. Thirty-two toxins displayed docking scores less than the co-crystalized ligand (XF7 = −9.5 kcal/mol). 2D docking poses showing key amino acids inside the Mpro ’s binding site are illustrated in Figure S1. Most of the scrutinized toxins demonstrated similar binding modes inside the Mpro ’s active site, forming a fundamental hydrogen bond with GLN189 and GLU166 (Figure S1). 2D chemical structures and calculated docking scores for those toxins are listed in Table 1. Philanthotoxin (T3D2489), a component of the venom of the Egyptian solitary wasp, demonstrated the highest binding affinity with a docking score of −11.7 kcal/mol, displaying a total of six hydrogen bonds with the key amino acid residues of Mpro (Table 1). Inspecting its binding mode showed that the ammonium group (NH3 + ) participates in a hydrogen bond with the backbone CO group of GLY170 with a bond length of 2.25 Å (Figure 2). Furthermore, two dimethylaminium groups participate in two hydrogen bonds with the backbone CO groups of PHE140 and ASN142 with bond lengths of 1.98 and 1.75 Å, respectively (Figure 2). The CO and NH of two N-methylacetamide groups display two hydrogen bonds with the backbone NH of GLU166 and the backbone carbonyl of GLN189 with bond lengths of 2.97 and 2.19 Å, respectively (Figure 2). The hydroxy group of the phenol ring interacts with the backbone NH of THR26 with a bond length of 1.99 Å (Figure 2). It is worth mentioning that gram-scale synthesis of philanthotoxin analogs has been obtained [28]. Azaspiracid (T3D2672), an alkaloid from Mytilus edulis (blue mussel) [29], exhibited the second-greatest binding affinity towards Mpro with a docking score of −11.6 kcal/mol (Table 1). Inspecting the binding mode of T3D2672 within the Mpro ’s binding pocket disclosed that the OH group of the tetrahydro-3,5-dimethyl-2H-pyran-2-ol ring and NH of 3,5-dimethylpiperidine ring exhibit two hydrogen bonds with the backbone CO group Pharmaceuticals 2022, 15, 153 4 of 22 of GLN189 with bond lengths of 2.1 and 2.18 Å, respectively (Figure 2). Furthermore, CO of the carboxylate group contributes a hydrogen bond with the backbone NH of ALA191 with a bond length of 2.23 Å (Figure 2). Taziprinone (T3D2378), a β-amino acid derivative, also showed a strong binding affinity against Mpro with an average value of −11.2 kcal/mol. The interaction was based in part on two hydrogen bonds with GLN189 and GLU166 with bond lengths of 2.32 and 1.85 Å, respectively (Figure 2). Figure 1. (a) 3D representation of the anticipated docking pose of XF7 (in pink) and experimental structure (in mauve) of XF7, (b) 3D, in addition to (c) 2D representations of the predicted binding mode of XF7 complexed with SARS-CoV-2 main protease (Mpro ). Pharmaceuticals 2022, 15, 153 5 of 22 Table 1. Estimated conventional and expansive docking scores (in kcal/mol), 2D chemical structures, and origin/usage for most promising potent toxins towards SARS-CoV-2 Mpro a . Compound Origin/ Name/Code Usage b Insect in toxin an (Egyptian solitary wasp) Marinetoxin Marine toxin T3D2672 (Mytilus uedulis) T3D2378 T3D2807 Industrial/ workplace toxin place toxin Synthetic Synthetic compound nd (anti- ety anxiety agent) Synthetic T3D2825 Conv. c Exp. d − 9.1 − − −9.5 − − − − − − XF7 XF7 T3D2489 Docking Score (kcal/mol) Chemical Structure Synthetic com- nd pound (antihyper(antihypertensive agent) T3D2874 Synthetic Synthetic compound nd (antiic allergic agent) T3D2938 Synthetic Synthetic compound nd (antiic allergic agent) − − − − − 11.7− − −11.7 − − − − − − − − −11.6 − − − − − − − − − − − − − − − − − − − − − − −11.2 − − − − − − − − − − − − − − −−10.9− − − − − − − − − − −10.9 − − − − − − − − − − − − − − − − − 10.9− − − − − − − − − −10.9 − − − − − − − − − − − − − − − − − − 10.8− − − − − − − − −10.9 − − − − − − − − − − − − − − − − − − − 10.5− − − − − − − −10.8 − − − − − − − − − − − − − − − − − − − − − − − − − − − − − − − − 11.3− −−11.2− − − − − − − Compound Origin/ Name/Code Usage b Conv. − − − − − − − − − − − − − − − − − − c Exp. d T3D2324 Industrial/ workplace toxin − − −9.2 − − − − − − −9.9 T3D2680 Synthetic com-nd pound (anticholesteremic agent) − − − − − − −9.8 − − −9.9 Synthetic compoundnd T3D2884 (antineoplastic plastic agent) agent) T3D4082 Plant Plant toxin toxin (Veratrum licalifornicum) ) Food toxin Food toxin (antihy(antihypoparT3D2694 poparathyathyroid roid agent) T3D2871 Food toxin (antipsychotic chotic agent) Plant toxin Plant T3D4050 T3D2536 toxin (Solanum chacoense) Animal Animal toxin toxin (B. rubescens, s, B. marinus) ) O OO O O O O OO O OO OO O O N N O OO O O N N O N N N N O O HO HO HO HO HO HO − −9.8 − − − − − − − − − − − − − − − − − − − − − − − − − −−9.8 − − − − − − − − − − − − − − − − − − − − − − − − − −−9.7 − − − − − − − − − − − − − − − − − − − − − − − − − − −9.7 − − − − − − − − − − − − − − − − − − − − − − − − − − − −9.6 − − − − − − − − − − − − − − − − − − − − − − − − − −9.6 − − − − − − − − − − − − − − − − − − − − − − − − − − − − − − − − − − − − − − − − − − − − − − − − − − − − − − − N N O − − − − − − N N − − − − − − Docking Score kcal/mol) Chemical Structure O O O O O O O O O −9.9 −9.8 −9.8 −9.8 −9.8 −9.8 Pharmaceuticals 2022, 15, 153 − − −− −− − − − − − − − −− − − Table 1. Cont. − Compound Origin/ Name/Code Usage b T3D2913 T3D4084 Chemical Structure − − −− − Plant Planttoxin toxin (genus − −− c −10.5 −− − − − Exp. − − − − − −− −−10.2 −− − − Veratrum) − −−− − − −Docking Score −(kcal/mol) − Conv. Synthetic Synthetic compound nd (antipsychotic chotic agent) agent) − d −−10.8 T3D2933 − − − −10.6 − Chemical Structure− − − −− − − Synthetic compound nd (antipsychotic chotic agent) O toxinO (Tiostrea chilensis) − O OO O O O N O − − −− − −Docking − Score − kcal/mol) − HO HO O HO O Exp. −9.6 − − −− − − −− − − − − −− −− −− − 9.5 − − − O HO c d −9.7 N N O HO 6 of 22 N N N O −− − Conv. O Marine toxin Marine T3D4051 −− − − Compound Origin/ Name/Code Usage b −− − O −9.7 O O HO Animal T3D2727 Synthetic Synthetic compound nd (antineoplastic plastic agent)agent) T3D2460 Synthetic com- nd pound (vasoconstrictor agent) T3D2750 Synthetic com- nd pound (vasoconstrictor agent) T3D2801 Synthetic Synthetic compound nd (psychotropic chotropic drug) − − − − − − − −10.0 − − − − −− − −− − − −10.1 − − − Planttoxin Plant toxin (genus −− − −−−− −− − −−9.9− −− − −− − −− − −10.1 − −− −−− −− −−9.8− − − − −− − − −10.0 −− −−− −− −− − −9.8− − −− − − − −− − −− − 10.1 − − − −10.5 − − −−10.1 − − − − −− − −− − −− − − −10.2 − −− − − − − − −−10.0 − −− − −− − −− − − − −10.2 −− − Animal toxintoxin − T3D2527 − − −− − − − − − 9.5 − − − (genus Dendroand bates and genus PhylloPhyllobates) −9.7 Synthetic Synthetic T3D4083 Veratrum) Synthetic T3D2939 T3D2143 − −− − − Synthetic com- nd pound (vasocon(vasoconstrictor agent) − −− − − Industrial/ workplace toxin place toxin T3D0233 Synthetic compoundnd (pesticide) Synthetic com- nd pound T3D2910 (cholinester(cholinesterase aseinhibitor) inhibitor) T3D2863 Synthetic com- nd pound (anti-HIV agent) Animal toxin Animal T3D2535 T3D4232 toxin (B. gargarizans) Bacterial toxin yano(cyanobacteria) − − − − −− − − −− − − − − − −− − − −− − − − − −− − − − − − − − − − −− − − − − − − − − − 9.5 − − − − −− −− − − − −− −− − − −9.7 − − 9.5 − −9.6 −9.5 − − − − − −9.6 − −− − − − − − −− − −9.4 − −− − −9.6 − −−− − − − − − −− − − − − 9.3 − − − − − −− −9. 6 −10.0 a Data sorted according to the expensive docking scores. b Taken from T3DB website 25 . c Conv. stands for the conventional docking computation. d Exp. stands for the expensive docking estimation. − − − − − − −− − − − Pharmaceuticals 2022, 15, 153 7 of 22 Figure 2. 3D and 2D molecular interactions of the predicted docking poses of toxins (i) T3D2489, (ii) T3D2672, and (iii) T3D2378 towards SARS-CoV-2 Mpro . 2.3. Molecular Dynamics (MD) Simulations Molecular dynamics (MD) simulations investigate the stabilization of the inhibitorreceptor complexes, conformational pliabilities, the reliability of inhibitor-receptor binding energies, and structural details [30]. Therefore, the most potent toxins with docking scores lower than −9.5 kcal/mol were subjected to MD simulations, followed by binding affinity estimations. The simulations were executed for 5 ns to reduce time and computational costs. The corresponding MM-GBSA binding affinities were computed and summarized in Table S2. From these data, it is apparent that seven toxins showed lower binding energies (∆Gbinding ) compared to −the co-crystallized XF7 inhibitor (calc. −40.1 kcal/mol)−(Figure − − 3). As a result, these toxins were selected and subjected to a 50 ns MD simulation to obtain more meticulous Mpro binding affinities. In addition, MM-GBSA binding energies were evaluated (Figure 3). Three out of these seven toxins—namely T3D2489, T3D2672, and T3D2378— displayed lower binding energies (∆Gbinding ) than the co-crystallized XF7 inhibitor (calc. −43.4 kcal/mol). The evaluated MM-GBSA binding affinities for T3D2489, T3D2672, and − − − − T3D2378 towards Mpro were −58.2, −54.7, and −48.7 kcal/mol over 50 ns MD simulations, respectively. To obtain more reliable binding affinities, longer MD simulations of 200 ns were executed for those three potent toxins in complex with Mpro , and the corresponding binding affinities were computed (Figure 3). − Pharmaceuticals 2022, 15, 153 Δ − − − 8 of 22 Figure 3. Calculated MM-GBSA binding energies for the native XF7 inhibitor and the most seven potent toxins complexed with SARS-CoV-2 Mpro throughout 5 ns, 50 ns, and 200 ns MD simulations. There is no noteworthy disparity between the computed MM-GBSA binding affinities throughout 50 ns and 200 ns MD simulations for T3D2489-, T3D2672-, and T3D2378-Mpro complexes (Figure 3). Compared to the binding energy of the co-crystallized XF7 (calc. −43.7 −kcal/mol), T3D2489, T3D2672, and T3D2378 revealed a greater binding affinity against Mpro throughout the 200 ns MD simulations with an average ∆GΔbinding of −−58.9, −−55.9, and − − 50.1 kcal/mol, respectively. To identify the principal driving forces in binding the identified toxins with Mpro , decomposition of the MM-GBSA binding affinities was carried out (Figure 4). Evdw was a significant contributor in the XF7-Mpro binding affinity with an average value ffi − value of of −54.5 kcal/mol; besides, Eele contribution was favorable with an average − −21.7 kcal/mol (Figure 4). For compounds T3D2489 and T3D2378, a predominance with a −value of −442.9 and −115.8 kcal/mol, respectively of Eele forces was observed − (Figure 4). Evdw was also favorable, with − an average − value of −53.3 and −51.2 kcal/mol for the T3D2489- and T3D2378-Mpro complexes, respectively (Figure 4). On the other hand, Evdw and Eele had approximately the same contribution in T3D2672-Mpro binding affinity with− average values of −68.9 and −65.4 kcal/mol (Figure 4). − To more fully scrutinize enzyme-inhibitor interactions and the participation of proximal amino acids in the inhibitor-enzyme complexes, total ∆Gbinding values were separated to individual residues with the assistance of an MM-GBSA approach (Figure 5); only residues with ∆Gbinding values lower than −0.50 kcal/mol were considered. GLY143, HIS164, GLU166, and GLN189 interact with T3D2489, T3D2672, T3D2378, and XF7. GLN189 contributed to the total binding affinity significantly with values of −3.7, −4.7, −4.7, and −3.2 kcal/mol for T3D2489-, T3D2672-, T3D2378-, and XF7-Mpro complexes, respectively (Figure 5). GLU166 was the second-highest contributor to total binding free with values of −0.7, −3.1, −1.2, and −2.3 kcal/mol for T3D2489-, T3D2672-, T3D2378-, and XF7-Mpro complexes, respectively (Figure 5). It is worth mentioning that all investigated complexes have similar interaction patterns with key amino acid residues, which signifies a resemblance in the binding mode in these complexes. Pharmaceuticals 2022, 15, 153 9 of 22 Δ Δ − ffi − − − − − − − − Figure 4. Components of the MM-GBSA binding energies for (a) XF7, (b) T3D2489, (c) T3D2672, and (d) T3D2378 complexed with SARS-CoV-2 Mpro throughout the simulation time of 200 ns. Δ Δ − ffi − − − − − − − − Figure 5. Energy participation of the proximal residues to the total binding free energy (kcal/mol) of T3D2489, T3D2672, T3D2378, and XF7 complexed with SARS-CoV-2 main protease (Mpro ). 2.4. Post-MD Analyses To further examine the constancy of T3D2489, T3D2672, and T3D2378 complexed with Mpro , structural and energetical analyses were executed over 200 ns MD simulations and compared to those of the co-crystalized XF7 inhibitor. Monitoring the structural steadiness of the scrutinized complexes was effectuated via inspecting hydrogen bond length, root-mean-square deviation (RMSD), binding energy per frame, and center-of-mass (CoM) distance. Pharmaceuticals 2022, 15, 153 10 of 22 2.4.1. Binding Energy per Frame The structural steadiness of T3D2489, T3D2672, T3D2378, and XF7 complexed with Mpro was comprehensively evaluated over the 200 ns MD simulations via mensuration correlations between binding affinity and time (Figure 6). An interesting aspect of binding energy per frame was the overall stabilities for T3D2489, T3D2672, T3D2378, and XF7 with average ∆G Δ binding of−−58.9,− −55.9, − −50.1, and − −43.7 kcal/mol, respectively. Based on this analysis, all investigated systems preserved constancy over the 200 ns MD simulations. Figure 6. Computed MM-GBSA binding energy per frame for XF7 (in black), T3D2489 (in red), T3D2672 (in blue), in addition to T3D2378 (in cyan) in complex with SARS-CoV-2 main protease (Mpro ) throughout 200 ns MD simulations. 2.4.2. Intermolecular Hydrogen Bonds Hydrogen bond analysis was used to estimate the constancy of hydrogen bonding between identified toxins and Mpro throughout a 200 ns MD simulation. The number of hydrogen bonds per frame was evaluated and depicted in Figure 7. The number of hydrogen bonds oscillated throughout the 200 ns MD simulations, and the average number of hydrogen bonds was four, two, and two for T3D2489-, T3D2672-, and T3D2378Mpro complexes. It is worth noting that XF7 showed the lowest number of hydrogen bonds with the proximal amino acids within Mpro ’s binding pocket, while the outstanding Δ − binding affinity of XF7 with average ∆Gbinding of −43.7 kcal/mol may be ascribed to other interactions like van der Waals and hydrophobic interactions. The dominance of van der Waals feature of interactions of XF7 with Mpro is in agreement with the MM-GBSA binding energies decomposition (Figure 4). The results obtained from the intermolecular bonds assured the presentence of a great stationary for the T3D2489, T3D2672, and T3D2378 complexed with Mpro compared to the XF7-Mpro complex. 2.4.3. Center-of-Mass Distance To obtain further in-depth insight into the steadiness of toxin-Mpro complexes throughout the 200 ns MD simulations, the center-of-mass (CoM) distance was estimated between the toxin and GLN189 (Figure 8). From the CoM graph, it is apparent that the measured CoM distance was more stable for T3D2672 and T3D2378 in complex with Mpro than T3D2489 and XF7 with average values of 5.1, 3.5, 10.2, and 5.6 Å, respectively. These findings demonstrate that the most identified toxins bind more tightly with the SARS-CoV-2 Mpro than XF7. Pharmaceuticals 2022, 15, 153 11 of 22 Figure 7. Number of hydrogen bonds formed between XF7 (in black), T3D2489 (in red), T3D2672 (in blue), and T3D2378 (in cyan) and SARS-CoV-2 Mpro throughout 200 ns MD simulations. Figure 8. Distance between the center-of-mass (CoM) (in Å) of XF7 (in black), T3D2489 (in red), T3D2672 (in blue), and T3D2378 (in cyan) and GLN189 of SARS-CoV-2 Mpro over the 200 ns MD simulations. 2.4.4. Root-Mean-Square Deviation The principal objective of the MD simulations is to inspect the positional and conformational changes of ligands upon binding to the binding pocket, which supplies insight into the binding steadiness. For ease of comparison, the root-mean-square deviation (RMSD) of the entire complex backbone atoms were evaluated to examine the structural stability of the T3D2489, T3D2672, T3D2378, and XF7 in complex with Mpro (Figure 9). Unambiguously, the evaluated RMSD values for the identified toxins complexed with Mpro stayed below 0.25 nm throughout the MD simulations (Figure 9). It is worth noting that RMSD analysis shows that the complexes stabilized after 10 ns and conserved their stabilities up to the end of the simulations. In general, the current findings confirmed that T3D2489, T3D2672, and T3D2378 are tightly bonded and do not leverage the structural constancy of the Mpro , in addition to keeping structural integrity. Pharmaceuticals 2022, 15, 153 12 of 22 constancy of the M , in addition to keeping structural integrity. Figure 9. Root-mean-square deviation (RMSD) of the backbone atoms from the initial structure of XF7 (in black), T3D2489 (in red), T3D2672 (in blue), and T3D2378 (in cyan) with SARS-CoV-2 Mpro during the simulation time of 200 ns. 2.5. Drug-like Properties The efficacy of therapeutic medicaments is substantially based on the molecular characteristics and bioactivity of the chemical compounds [31]. To predict the drug-like and bioactivity of the identified toxins as SARS-CoV-2 inhibitors, a Molinspiration tool was employed to estimate the in silico molecular characteristics. The anticipated physiochemical properties are summarized in Table 2. Promising drug-likeness properties were observed, except for T3D2672, which showed violations in some parameters such as mi logP, TPSA, number of hydrogen bonds acceptors (nON), and molecular weight. The mi logP values of T3D2489, T3D2378, and XF7 were promising, with values less than 5 [32]. The TPSA values of the T3D2489, T3D2378, and XF7 were less than 140 Å, revealing that the compounds have eminent oral absorption or membrane permeability [33]. In addition, the number of hydrogen bond acceptors (nON) was lower than 10, while the number of hydrogen bond donors (nOHNH) was lower than 5, except for T3D2489. The molecular weights for T3D2489, T3D2378, and XF7 were 435.6, 385.5, and 492.5 daltons, respectively, proposing that most investigated toxins have good absorption and/or permeation across the cell membrane. Table 2. Portended physiochemical parameters and structural descriptors of T3D2489, T3D2672, T3D2378, and XF7 as SARS-CoV-2 main protease (Mpro ) inhibitors. Compound Name/Code mi logP TPSA nON nOHNH Nrotb MWt %ABS T3D2489 T3D2672 T3D2378 XF7 0.03 6.7 1.7 4.3 128.5 163.7 61.9 109.9 8 13 6 8 7 4 1 2 18 9 4 6 435.6 842.1 385.5 492.5 64.7% 52.5% 87.6% 71.1% 2.6. In Silico ADMET Analysis ADMET analysis is especially helpful in simplifying clinical trials, especially in the early stage of drug design [34]. GI absorption, skin, Caco2 permeability, and aqueous solubility are absorption properties needed to be considered in any drug discovery process [35]. It is signified that a GI absorption value greater than 30% indicates perfect absorbance. T3D2489 (47.6%), T3D2672 (56.7%), T3D2378 (95.0%), and XF7 (93.5%) manifested good absorbance rates (Table 3). The identified toxins revealed appropriate skin permeability, Pharmaceuticals 2022, 15, 153 13 of 22 with a skin permeability value higher than −2.5 cm/h. The identified toxins also showed low Caco2 permeability (<0.9 cm/s). Another substantial agent through ADMET analysis was to anticipate the P-glycoprotein non-substrate. All toxins were identified as a substrate for P-glycoprotein (Table 3). Table 3. Anticipated ADMET characteristics of the top potent inhibitors. ADME Parameters XF7 T3D2489 Water solubility Caco2 permeability Intestinal absorption (human) Skin Permeability P-glycoprotein substrate P-glycoprotein I inhibitor P-glycoprotein II inhibitor −3.5 0.5 93.5 −2.7 Yes Yes Yes T3D2672 T3D2378 −3.0 −0.2 47.6 −2.7 Yes No No −3.4 −0.1 56.7 −2.7 Yes Yes Yes −2.5 0.3 95.0 −3.2 Yes No No 1.6 −0.7 −4.1 0.6 −1.7 −3.4 1.3 0.2 −2.6 No No No No No No No No Yes No No No No No No Yes No No No No No 1.4 −0.1 0.8 No 0.3 No Yes 2.7 3.1 Yes No 0.3 2.2 No −0.3 No No 2.8 2.7 Yes No 0.3 0.9 No −0.4 No Yes 2.4 0.7 Yes No 0.6 4.7 Absorption Distribution VDss (human) BBB permeability CNS permeability 0.3 −1.1 −2.5 Metabolism CYP2D6 substrate CYP3A4 substrate CYP1A2 inhibitior CYP2C19 inhibitior CYP2C9 inhibitior CYP2D6 inhibitior CYP3A4 inhibitior No Yes No Yes Yes No Yes Excretion Total Clearance 0.8 AMES toxicity Max. tolerated dose (human) hERG I inhibitor hERG II inhibitor Oral Rat Acute Toxicity (LD50) Oral Rat Chronic Toxicity (LOAEL) Hepatotoxicity Skin Sensitisation T. Pyriformis toxicity Minnow toxicity No 0.4 No Yes 3.3 1.1 Yes No 0.3 2.4 Toxicity In order to examine the drug distribution, the VDss, BBB membrane permeability, and CNS were assessed [32]. Higher distribution volumes were observed for T3D2489, T3D2672, and T3D2378 with log VDss values of 1.6, 0.6, and 1.3, respectively (Table 3). For BBB membrane permeability, log BB values in the range of −1.0 to 0.3 implied that the drug candidates passed the BBB membrane. For CNS permeability, log PS values ranged from −3 to −2, pointing out impenetrability. Most of the investigated inhibitors were forecasted to be neither able to permeate the CNS nor pass the BBB membrane (Table 3). CYP450 has a substantial role in drug metabolism in the liver system [36]. The metabolism scores revealed that all the inspected inhibitors did not prohibit CYP2D6 enzymes and did not perform as inhibitors for CYP3A4, CYP2C9, and CYP2C19 enzymes, except for the reference compound. The total drug clearance was inspected via a collection of hepatic and renal clearance. Total clearance construes the drug concentration in the body utilizing its removal rate. The data indicated that compound excretion rates range from −0.1 to 1.4 mL/min/kg (Table 3). In drug design, toxicity is a significant criterion and plays a remarkable role in selecting of sufficient drug candidates [32]. All the drug candidates in this analysis have not passed any skin allergic action and hepatotoxic influence (Table 3). hERG inhibition (I and II) is a Pharmaceuticals 2022, 15, 153 14 of 22 fundamental agent for toxicity analysis in addition to it also including cardiotoxicity. None of the inhibitors displayed inhibitory behaviors for hERG I. T3D2489, T3D2378, and XF7 were foretold to be hERG II inhibitors. All the investigated inhibitors had not crossed any Tetrahymena pyriformis as well as AMES toxicities. The maximum tolerated dosage range, lowest-observed-adverse-effect level (LOAEL), and LD50 were expected via the toxicity analysis server and the anticipated scores as summarized in Table 3. Consequently, this study concluded that these bioactive identified toxins could be utilized as possible Mpro inhibitor candidates. 2.7. Molecular Target Prediction and Network Analysis The lysosomal cathepsins, primarily cathepsin L (CTSL) and cathepsin B (CTSB) cleave and activate S proteins, which then merge with host cells [37,38]. SARS-CoV-2 upregulates CTSL expression both in vivo and in vitro [39]. This increases pseudo-virus infection in human cells. CTSL may be a therapeutic target to remedy COVID-19 disease [40,41]. CatL is a secreted lysosomal protease that is also known as a lysosomal protease. Excessive production of cysteine cathepsins has been related to various clinical conditions, including inflammation. As a result, high quantities of CatL might be activated in inflammatory cells under inflammatory circumstances. Moreover, because this enzyme is involved in the genesis, progression, and metastasis of cancer, its inhibition might be beneficial in treatment [42]. The plasma CatL has been recommended as a marker for pancreatic cancer. Philanthotoxin (T3D2489) was identified as an inhibitor against Mpro and provided a ligand lead using the STRING database, including the biochemical top 20 signaling genes CTSB, CTSL, and CTSK (Figure 10 and Table S3). Figure 10. Bioinformatic analysis: (A) The PPI network of top genes; (B) the top 20 gene comparison in lung adenocarcinoma and lung squamous cell carcinoma; (C,D) Gene Expression Profiling Interactive Analysis (GEPIA), the expression level of CTSB was increased in lung adenocarcinoma and squamous cell lung carcinoma. 2.8. Pathway Enrichment Analysis (PEA) and Reactome Mining For extensive and biological-wide mining of philanthotoxin (T3D2489) target-function interactions, a Reacfoam map was constructed based on PEA analysis and Reactome mining/modeling. The Reacfoam map tree was constructed to illustrate the top enriched Pharmaceuticals 2022, 15, 153 15 of 22 pathways affected by 20 top gene targets in response to philanthotoxin (Figure 11). Moreover, an illustrative graphical map was built to show the top enriched pathway as a response to T3D2489 treatment (Figure 12). Prominently, the top enriched pathways involved: (1) signaling by the interleukin pathway, (2) interleukin-4 and interleukin-13 signaling pathway, and (3) immune system pathway. These pathways were determined to be the most significantly enriched pathways targeted by T3D2489 (Table S3). Figure 11. The Reacfoam map shows the top enriched pathway (Interleukin-4 and Interleukin-13 signaling) influenced by the top 20 gene targets in response to philanthotoxin (T3D2489). The SARS-CoV-2 infection causes an induced inflammation and generates subsequent higher levels of immune cell infiltration and cytokines that trigger matrix metalloproteinases (MMPs) activation. Interestingly, the PEA-Reactome mining coupled approach revealed higher enrichment of signaling by Interleukins (particularly: Interleukin-4 and Interleukin-13 signaling), activation of matrix metalloproteinases signaling pathway, and immune system as a result of philanthotoxin induction among all other human biological pathways. Interestingly, remodeling of interleukins (IL-4 and IL-13) can transform growth factor-beta (TGF-β) levels, and consequently, the number of M2 macrophages in SARSCoV-2 patients, which ultimately can prevent severe lung injury [43]. Furthermore, a recent report emphasized that gathering and remodeling the immune cells, matrix metalloproteinases, secreted cytokines (especially interleukin-4 and Interleukin-13), and several other mediators is proposed as a feasible option for treating COVID-patients [44]. Pharmaceuticals 2022, 15, 153 16 of 22 Figure 12. Graphic representation of the Interleukin-4 and Interleukin-13 signaling Reactome pathway influenced as a response to philanthotoxin (T3D2489) in the human genome. 3. Computational Methodology 3.1. Target Preparation The X-ray resolved three-dimensional structure of SARS-CoV-2 main protease (Mpro ) (PDB ID: 7L13, resolution: 2.17 Å) in complex with a noncovalent inhibitor (XF7) was retrieved and utilized as a template for all in silico calculations [27]. The viral target was prepared by eradicating all heteroatoms involving ligands, water molecules, and ions. The protonation states of the titratable amino acids were investigated and assigned using the H++ web server [45]. In addition, all missing hydrogen atoms are topologically added. In H++ estimations, physiologic conditions of 10, 0.15, 80, and 6.5 for internal dielectric constant, salinity, external dielectric constant, and PH, respectively. 3.2. Database Preparation The Toxin and Toxin Target Database (T3DB) was downloaded and prepared for virtual screening [25]. All the molecules were retrieved in a 2D structural data format (SDF). Omega2 software was then utilized to construct the 3D chemical structures [46,47]. A conformational search was executed to generate all conformers with an energy window of 10 kcal/mol. The lowest energy conformer was then minimized using an MMFF94S force field available within SZYBKI software [48,49]. Tautomer and fixpka applications implemented inside QUACPAC software were applied to investigate the tautomer enumeration and the protonation state of the toxins [50]. The partial atomic charges of each compound within the T3DB database were determined using a Gasteiger-Marsili method [51]. Duplicated molecules with congruent international chemical identifier keys (InChIKey) were removed [52]. Prepared T3DB data files are available through www.compchem.net/ccdb. An illustrative diagram of the employed computational approaches for the screening process of the T3DB database is illustrated in Figure S2. 3.3. Molecular Docking AutoDock4.2.6 software was utilized to conduct all molecular docking calculations [53]. On the basis of the AutoDock protocol [54], the MGL tools (version 1.5.7) were used to generate the pdbqt file for the SARS-CoV-2 Mpro . For molecular docking calculations, the number of generations and population size were set to 27,000 and 300, respectively. The maximum number of energy evaluations (eval) and the number of genetic algorithms (GA) run variables were adjusted to 5,000,000, and 25,000,000, and 50, and 250 for con- Pharmaceuticals 2022, 15, 153 17 of 22 ventional and expensive molecular docking calculations, respectively. The rest of the docking parameters were maintained at the default settings. The grid box size was set at 60 Å × 60 Å × 60 Å, which is capable of accommodating the entire Mpro ’s binding site. Grid map files with a spacing value of 0.375 Å were established utilizing AutoGrid 4.2.6. The grid was centered at the coordinates X = −13.069, Y = 9.740, and Z = 68.490. A genetic algorithm inside the AutoDock software was employed to evaluate the various inhibitor conformers. Conformations were clustered via the root-mean-square deviation (RMSD) tolerance of 1.0 Å and were sorted based on the docking score [55]. Besides, the lowest docking score within the largest cluster was considered for opting as the representative pose. 3.4. Molecular Dynamics Simulations AMBER16 software was used to conduct all molecular dynamics (MD) simulations for the most potent toxins in complex with SARS-CoV-2 Mpro [56]. A general AMBER force field (GAFF2) [57] was employed to characterize the investigated toxins, whereas the AMBER force field of 14SB [58] was adopted for the characterization of the viral enzyme. The restrained electrostatic potential (RESP) fitting approach was used to estimate the atomic partial charge of the studied toxins at the HF/6-31G* level with the assistance of Gaussian09 software [59,60]. A solvated octahedron box of TIP3P water model with 12 Å distances between the box edge and atoms of the toxin-Mpro complexes was constructed. The sodium (Na+ ) and chloride (Cl− ) counter-ions were inserted to neutralize all solvated systems, as well as to preserve the isosmotic condition (0.15 M NaCl). The prepared systems were initially minimized through combined steepest and conjugate gradient algorithms for 5000 steps to remove any steric clashes or inappropriate geometries. The minimized systems were thereafter gradually heated to 300 K over a brief interval of 50 ps with a weak constraint of 10 kcal mol−1 Å−1 on the amino acid residues. To guarantee a reasonable initial structure, an equilibration stage was executed over a total duration of 1000 ps under a constant number of particles, pressure (1 atm), and temperature (NPT) ensemble. Eventually, the production runs were conducted over simulation times of 5 ns, 50 ns, and 200 ns. Snapshots were recorded and saved each 10 ps for post-MD analyses. A cutoff distance of 12 Å was used for the non-bonded interactions. The Particle-Mesh Ewald (PME) algorithm was utilized to estimate the long-range electrostatic interactions [61]. The collision frequency (gamma_ln) 1.0 ps−1 was applied to maintain the temperature at 298 K [62]. The Berendsen barostat with a pressure relaxation time of 2 ps was utilized to keep pressure constant [63]. The SHAKE algorithm with 2 fs integration step was employed to restrict all bonds including hydrogen atoms [64]. All MD simulations were carried out utilizing the CUDA version of pmemd in AMBER 16 software. All molecular docking calculations, MD simulations, as well as quantum mechanics (QM) evaluations were performed on the CompChem GPU/CPU hybrid cluster (hpc.compchem.net). All molecular interactions were visualized via the Discovery Studio module of Biovia software [65]. 3.5. MM-GBSA Binding Energy A molecular mechanic generalized Born surface area (MM-GBSA) approach was adopted to calculate binding free energies of the investigated toxins in complex with SARS-CoV-2 Mpro . The MM-GBSA (∆Gbinding ) binding energies were computed based on uncorrelated snapshots collected over the MD simulations as follows: ∆Gbinding = Gcomplex − (Gtoxin + GMpro ) (1) where the energy term (G) is calculated as: G = Evdw + Eele + GGB + GSA (2) Eele and Evdw stand for electrostatic and van der Waals energies, respectively. GSA implies the nonpolar solvation-free energy, generally computed with a linear relation to Pharmaceuticals 2022, 15, 153 18 of 22 the solvent-accessible surface area (SASA). GGB indicates the electrostatic solvation free energy estimated from the generalized Born equation. In this study, the modified GB model developed via Onufriev et al. (igb = 2) was used [66]. A single-trajectory method was utilized, in which the coordinates of every toxin-Mpro , toxin, and Mpro were calculated from a single trajectory. Due to the expensive computational demand, entropy estimations were neglected [67,68]. 3.6. Drug-Likeness Properties For the most promising toxins, the physicochemical properties were estimated utilizing an online Molinspiration cheminformatics package (http://www.molinspiration.com) according to Lipinski Rules. Molinspiration supports for prediction of significant physicochemical properties such as molecular weight (MW), topological polar surface area (TPSA), octanol/water partition coefficient (mi logP), hydrogen bond donor (HBD), hydrogen bond acceptor (HBA), rotatable bond count (RB), and percent absorption (%ABS). %ABS was computed as follows [69]: %ABS = 109 − [0.345 × TPSA] (3) 3.7. In Silico ADMET Analysis A freely accessible web server, pkCSM (http://biosig.unimelb.edu.au/pkcsm/prediction) is an in silico tool for anticipating substantial pharmacokinetic properties [34]. ADMET properties involve absorption (A): human intestinal absorption (HIA), water-solubility, P-glycoprotein I and II inhibitors, Caco-2 permeability, P-glycoprotein substrate, skin permeability. Distribution (D) is anticipated according to blood-brain barrier (BBB) permeability, steady-state volume of distribution (VDss), fraction unbound, and central nervous system (CNS) permeability. CYP2D6/CYP3A4 substrate is used to detect the metabolism (M). Excretion (E) is determined via drug total clearance. Subsequently, the toxicity (T) of the drugs is predicated on the Human Ether-a-go-go-related gene inhibition, carcinogenic status, mutagenic status, as well as acute oral toxicity [31]. 3.8. Protein Interaction Network Analyses To identify the probable targets for each ligand, the target molecules were tested using SwissTargetPrediction (http://www.swisstargetprediction.ch), an internet website-based program. We obtained the top one hundred genes for the possible metabolite. After that, a functional STRING database for the most likely targets was utilized to build proteinprotein interactions (PPI). According to network architecture, Cytoscape 3.8.2 was utilized to analyze all possible receptor-function relationships. Pathway enrichment analysis was also performed utilizing Cytoscape 3.8.2 to scrutinize all prospective receptor-function connections for the top 20 targeted genes. The top 20 genes were investigated further using a newly developed interactive webserver (Gene Expression Profiling Interactive Analysis, GEPIA, http://gepia.cancer-pku.cn/index.html). Besides, to investigate all potential targetfunction relationships for the top most 20 targeted genes and their biological influencedpathways/networks, the pathway enrichment analysis (PEA) approach was conducted with the assistance of Cytoscape 3.8.2 software [70]. Subsequently, the Reactome mining analysis and visualization were achieved using the ReactomeFIViz tool towards modeling and annotating all the philanthotoxin (T3D2489)-target interactions [71]. 4. Conclusions Herein, a toxin and toxin-target database (T3DB) was mined to identify potential SARSCoV-2 Mpro inhibitors utilizing combined molecular docking and molecular dynamics simulations. On the basis of molecular docking calculations and MD simulations combined with molecular mechanics-generalized born surface area binding energy calculations, three compounds, T3D2489, T3D2672, and T3D2378, demonstrated promising binding affinity with ∆Gbinding < −50.0 kcal/mol towards Mpro . The energetic and structural analyses Pharmaceuticals 2022, 15, 153 19 of 22 during 200 ns MD simulations pointed to great constancy for these compounds in complex with Mpro . Additionally, the identified toxins demonstrated favorable pharmacokinetic and pharmacodynamic properties. The results obtained from PEA analysis combined Reactome-mining presented an interesting primary role of philanthotoxin in remodeling the interleukins (especially IL-4, IL -13) and matrix metalloproteinases (MMPs), which could alleviate lung injury in COVID-19 patients. In vitro and in vivo evaluations are planned to elucidate further the role of these compounds as prospective drug candidates and validate the computational findings. Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/ph15020153/s1, Figure S1: 2D representations of the binding modes of the thirty-two potent toxins complexed with SARS-CoV-2 main protease (Mpro ); Figure S2: Schematic representation of the utilized in silico techniques and the filtration process; Table S1: Estimated conventional and expansive docking scores (in kcal/mol), for top 200 toxins towards SARS-CoV-2 main protease (Mpro ); Table S2: Estimated conventional and expansive docking scores (in kcal/mol), and MM-GBSA binding energies (in kcal/mol) over 5 ns MD simulations for XF7 and the top 32 potent toxins towards SARS-CoV-2 main protease (Mpro ); Table S3: Top 20 enriched pathways influenced by philanthotoxin (T3D2489) targets resulted from PEA analysis. Author Contributions: Conceptualization, M.A.A.I. and L.A.J.-A.; Data curation, A.H.M.A.; Formal analysis, A.H.M.A. and M.A.M.A.; Investigation, A.H.M.A., M.A.M.A. and P.A.S.; Methodology, M.A.A.I., O.R.A., M.-E.F.H. and P.A.S.; Project administration, M.A.A.I. and L.A.J.-A.; Resources, M.A.A.I.; Software, M.A.A.I.; Supervision, M.A.A.I.; Visualization, A.H.M.A., M.A.M.A. and P.A.S.; Writing—original draft, A.H.M.A. and P.A.S.; Writing—review & editing, M.A.A.I., L.A.J.-A., M.A.M.A., O.R.A., M.N.A., M.S.M., M.E.S.S., A.M.S., P.W.P. and M.-E.F.H. All authors have read and agreed to the published version of the manuscript. Funding: This research received no external funding. Institutional Review Board Statement: Not applicable. Informed Consent Statement: Not applicable. Data Availability Statement: The data presented in this study are available in the supplementary material. Acknowledgments: Ahmed M. Shawky would like to thank the Deanship of Scientific Research at Umm Al-Qura University for supporting this work by Grant: 22UQU433174DSR01. Laila A. JaraghAlhadad acknowledges and appreciates both Research Sector at Kuwait University and Learner Research Institute at Cleveland Clinic for their support. The computational work was completed with resources supported by the Science and Technology Development Fund, STDF, Egypt, Grants No. 5480 & 7972 (Granted to Mahmoud A. A. Ibrahim). Conflicts of Interest: The authors declare that they have no conflict of interest. References 1. 2. 3. 4. 5. 6. 7. 8. Wu, F.; Zhao, S.; Yu, B.; Chen, Y.M.; Wang, W.; Song, Z.G.; Hu, Y.; Tao, Z.W.; Tian, J.H.; Pei, Y.Y.; et al. A new coronavirus associated with human respiratory disease in China. Nature 2020, 579, 265–269. [CrossRef] [PubMed] Gorbalenya, A.E.; Baker, S.C.; Baric, R.S.; de Groot, R.J.; Drosten, C.; Gulyaeva, A.A.; Haagmans, B.L.; Lauber, C.; Leontovich, A.M.; Neuman, B.W.; et al. The species Severe acute respiratory syndrome-related coronavirus: Classifying 2019-nCoV and naming it SARS-CoV-2. Nat. Microbiol. 2020, 5, 536–544. Huang, C.; Wang, Y.; Li, X.; Ren, L.; Zhao, J.; Hu, Y.; Zhang, L.; Fan, G.; Xu, J.; Gu, X.; et al. Clinical features of patients infected with 2019 novel coronavirus in Wuhan, China. Lancet 2020, 395, 497–506. [CrossRef] Zhu, N.; Zhang, D.; Wang, W.; Li, X.; Yang, B.; Song, J.; Zhao, X.; Huang, B.; Shi, W.; Lu, R.; et al. A Novel Coronavirus from Patients with Pneumonia in China, 2019. N. Engl. J. Med. 2020, 382, 727–733. [CrossRef] WHO Coronavirus Disease (COVID-19) Dashboard. Available online: https://covid19.who.int/ (accessed on 29 December 2021). WHO COVID-19 Vaccines. Available online: https://www.who.int/emergencies/diseases/novel-coronavirus-2019/covid-19 -vaccines (accessed on 29 December 2021). Li, Q.; Kang, C. Progress in Developing Inhibitors of SARS-CoV-2 3C-Like Protease. Microorganisms 2020, 8, 1250. [CrossRef] Suarez, D.; Diaz, N. SARS-CoV-2 Main Protease: A Molecular Dynamics Study. J. Chem. Inf. Model. 2020, 60, 5815–5831. [CrossRef] Pharmaceuticals 2022, 15, 153 9. 10. 11. 12. 13. 14. 15. 16. 17. 18. 19. 20. 21. 22. 23. 24. 25. 26. 27. 28. 29. 30. 31. 32. 33. 20 of 22 Verma, S.; Twilley, D.; Esmear, T.; Oosthuizen, C.B.; Reid, A.M.; Nel, M.; Lall, N. Anti-SARS-CoV Natural Products With the Potential to Inhibit SARS-CoV-2 (COVID-19). Front. Pharmacol. 2020, 11, 561334. [CrossRef] Ibrahim, M.A.A.; Abdelrahman, A.H.M.; Allemailem, K.S.; Almatroudi, A.; Moustafa, M.F.; Hegazy, M.F. In Silico Evaluation of Prospective Anti-COVID-19 Drug Candidates as Potential SARS-CoV-2 Main Protease Inhibitors. Protein J. 2021, 40, 296–309. [CrossRef] Ibrahim, M.A.A.; Abdelrahman, A.H.M.; Hegazy, M.F. In-silico drug repurposing and molecular dynamics puzzled out potential SARS-CoV-2 main protease inhibitors. J. Biomol. Struct. Dyn. 2021, 39, 5756–5767. [CrossRef] Caly, L.; Druce, J.D.; Catton, M.G.; Jans, D.A.; Wagstaff, K.M. The FDA-approved drug ivermectin inhibits the replication of SARS-CoV-2 in vitro. Antivir. Res. 2020, 178, 104787. [CrossRef] Liu, J.; Cao, R.Y.; Xu, M.Y.; Wang, X.; Zhang, H.Y.; Hu, H.R.; Li, Y.F.; Hu, Z.H.; Zhong, W.; Wang, M.L. Hydroxychloroquine, a less toxic derivative of chloroquine, is effective in inhibiting SARS-CoV-2 infection in vitro. Cell Discov. 2020, 6, 16. [CrossRef] [PubMed] Ibrahim, M.A.A.; Abdelrahman, A.H.M.; Hussien, T.A.; Badr, E.A.A.; Mohamed, T.A.; El-Seedi, H.R.; Pare, P.W.; Efferth, T.; Hegazy, M.F. In silico drug discovery of major metabolites from spices as SARS-CoV-2 main protease inhibitors. Comput. Biol. Med. 2020, 126, 104046. [CrossRef] Ibrahim, M.A.A.; Mohamed, E.A.R.; Abdelrahman, A.H.M.; Allemailem, K.S.; Moustafa, M.F.; Shawky, A.M.; Mahzari, A.; Hakami, A.R.; Abdeljawaad, K.A.A.; Atia, M.A.M. Rutin and flavone analogs as prospective SARS-CoV-2 main protease inhibitors: In silico drug discovery study. J. Mol. Graph. Model. 2021, 105, 107904. [CrossRef] [PubMed] Ibrahim, M.A.A.; Abdelrahman, A.H.M.; Atia, M.A.M.; Mohamed, T.A.; Moustafa, M.F.; Hakami, A.R.; Khalifa, S.A.M.; Alhumaydhi, F.A.; Alrumaihi, F.; Abidi, S.H.; et al. Blue Biotechnology: Computational Screening of Sarcophyton Cembranoid Diterpenes for SARS-CoV-2 Main Protease Inhibition. Mar Drugs 2021, 19, 391. [CrossRef] [PubMed] Ibrahim, M.A.A.; Abdelrahman, A.H.M.; Mohamed, T.A.; Atia, M.A.M.; Al-Hammady, M.A.M.; Abdeljawaad, K.A.A.; Elkady, E.M.; Moustafa, M.F.; Alrumaihi, F.; Allemailem, K.S.; et al. In Silico Mining of Terpenes from Red-Sea Invertebrates for SARS-CoV-2 Main Protease (M(pro)) Inhibitors. Molecules 2021, 26, 2082. [CrossRef] Ibrahim, M.A.A.; Abdeljawaad, K.A.A.; Abdelrahman, A.H.M.; Hegazy, M.F. Natural-like products as potential SARS-CoV-2 M(pro) inhibitors: In-silico drug discovery. J. Biomol. Struct. Dyn. 2021, 39, 5722–5734. [CrossRef] Zakaryan, H.; Arabyan, E.; Oo, A.; Zandi, K. Flavonoids: Promising natural compounds against viral infections. Arch. Virol. 2017, 162, 2539–2551. [CrossRef] Cherrak, S.A.; Merzouk, H.; Mokhtari-Soulimane, N. Potential bioactive glycosylated flavonoids as SARS-CoV-2 main protease inhibitors: A molecular docking and simulation studies. PLoS ONE 2020, 15, e0240653–e0240666. [CrossRef] Jo, S.; Kim, S.; Kim, D.Y.; Kim, M.-S.; Shin, D.H. Flavonoids with inhibitory activity against SARS-CoV-2 3CLpro. J. Enzym. Inhib. Med. Chem. 2020, 35, 1539–1544. [CrossRef] Fact Sheet for Healthcare Providers: Emergency Use Authorization for PAXLOVID. Available online: https://www.fda.gov/ media/155050/download (accessed on 20 December 2021). Harvey, A.L. Toxins and drug discovery. Toxicon 2014, 92, 193–200. [CrossRef] Bordon, K.C.F.; Cologna, C.T.; Fornari-Baldo, E.C.; Pinheiro-Junior, E.L.; Cerni, F.A.; Amorim, F.G.; Anjolette, F.A.P.; Cordeiro, F.A.; Wiezel, G.A.; Cardoso, I.A.; et al. From Animal Poisons and Venoms to Medicines: Achievements, Challenges and Perspectives in Drug Discovery. Front. Pharmacol. 2020, 11, 1132–1160. [CrossRef] [PubMed] Lim, E.; Pon, A.; Djoumbou, Y.; Knox, C.; Shrivastava, S.; Guo, A.C.; Neveu, V.; Wishart, D.S. T3DB: A comprehensively annotated database of common toxins and their targets. Nucleic Acids Res. 2010, 38, D781–D786. [CrossRef] [PubMed] Mohamed, T.A.; Elshamy, A.I.; Ibrahim, M.A.A.; Atia, M.A.M.; Ahmed, R.F.; Ali, S.K.; Mahdy, K.A.; Alshammari, S.O.; Al-Abd, A.M.; Moustafa, M.F.; et al. Gastroprotection against Rat Ulcers by Nephthea Sterol Derivative. Biomolecules 2021, 11, 1247. [CrossRef] Zhang, C.H.; Stone, E.A.; Deshmukh, M.; Ippolito, J.A.; Ghahremanpour, M.M.; Tirado-Rives, J.; Spasov, K.A.; Zhang, S.; Takeo, Y.; Kudalkar, S.N.; et al. Potent Noncovalent Inhibitors of the Main Protease of SARS-CoV-2 from Molecular Sculpting of the Drug Perampanel Guided by Free Energy Perturbation Calculations. ACS Cent. Sci. 2021, 7, 467–475. [CrossRef] Wellendorph, P.; Jaroszewski, J.W.; Hansen, S.H.; Franzyk, H. A sequential high-yielding large-scale solution-method for synthesis of philanthotoxin analogues. Eur. J. Med. Chem. 2003, 38, 117–122. [CrossRef] Satake, M.; Ofuji, K.; Naoki, H.; James, K.J.; Furey, A.; McMahon, T.; Silke, J.; Yasumoto, T. Azaspiracid, a new marine toxin having unique spiro ring assemblies, isolated from Irish mussels, Mytilus edulis. J. Am. Chem. Soc. 1998, 120, 9967–9968. [CrossRef] De Vivo, M.; Masetti, M.; Bottegoni, G.; Cavalli, A. Role of molecular dynamics and related methods in drug discovery. J. Med. Chem. 2016, 59, 4035–4061. [CrossRef] Shen, M.; Tian, S.; Li, Y.; Li, Q.; Xu, X.; Wang, J.; Hou, T. Drug-likeness analysis of traditional Chinese medicines: 1. property distributions of drug-like compounds, non-drug-like compounds and natural compounds from traditional Chinese medicines. J. Cheminform. 2012, 4, 31. [CrossRef] Han, Y.; Zhang, J.; Hu, C.Q.; Zhang, X.; Ma, B.; Zhang, P. In silico ADME and Toxicity Prediction of Ceftazidime and Its Impurities. Front. Pharmacol. 2019, 10, 434. [CrossRef] Bakht, M.A.; Yar, M.S.; Abdel-Hamid, S.G.; Al Qasoumi, S.I.; Samad, A. Molecular properties prediction, synthesis and antimicrobial activity of some newer oxadiazole derivatives. Eur. J. Med. Chem. 2010, 45, 5862–5869. [CrossRef] Pharmaceuticals 2022, 15, 153 34. 35. 36. 37. 38. 39. 40. 41. 42. 43. 44. 45. 46. 47. 48. 49. 50. 51. 52. 53. 54. 55. 56. 57. 58. 59. 60. 61. 62. 63. 21 of 22 Pires, D.E.; Blundell, T.L.; Ascher, D.B. pkCSM: Predicting Small-Molecule Pharmacokinetic and Toxicity Properties Using Graph-Based Signatures. J. Med. Chem. 2015, 58, 4066–4072. [CrossRef] [PubMed] Dahlgren, D.; Lennernas, H. Intestinal Permeability and Drug Absorption: Predictive Experimental, Computational and In Vivo Approaches. Pharmaceutics 2019, 11, 411. [CrossRef] [PubMed] Zanger, U.M.; Schwab, M. Cytochrome P450 enzymes in drug metabolism: Regulation of gene expression, enzyme activities, and impact of genetic variation. Pharmacol. Ther. 2013, 138, 103–141. [CrossRef] [PubMed] Zhou, Y.; Vedantham, P.; Lu, K.; Agudelo, J.; Carrion, R., Jr.; Nunneley, J.W.; Barnard, D.; Pohlmann, S.; McKerrow, J.H.; Renslo, A.R.; et al. Protease inhibitors targeting coronavirus and filovirus entry. Antiviral. Res. 2015, 116, 76–84. [CrossRef] Ou, X.; Liu, Y.; Lei, X.; Li, P.; Mi, D.; Ren, L.; Guo, L.; Guo, R.; Chen, T.; Hu, J.; et al. Characterization of spike glycoprotein of SARS-CoV-2 on virus entry and its immune cross-reactivity with SARS-CoV. Nat. Commun. 2020, 11, 1620. [CrossRef] Zhao, M.M.; Yang, W.L.; Yang, F.Y.; Zhang, L.; Huang, W.J.; Hou, W.; Fan, C.F.; Jin, R.H.; Feng, Y.M.; Wang, Y.C.; et al. Cathepsin L plays a key role in SARS-CoV-2 infection in humans and humanized mice and is a promising target for new drug development. Signal Transduct. Target. Ther. 2021, 6, 134. [CrossRef] [PubMed] Liu, T.; Luo, S.; Libby, P.; Shi, G.P. Cathepsin L-selective inhibitors: A potentially promising treatment for COVID-19 patients. Pharmacol. Ther. 2020, 213, 107587. [CrossRef] Smieszek, S.P.; Przychodzen, B.P.; Polymeropoulos, M.H. Amantadine disrupts lysosomal gene expression: A hypothesis for COVID19 treatment. Int. J. Antimicrob. Agents 2020, 55, 106004. [CrossRef] Gomes, C.P.; Fernandes, D.E.; Casimiro, F.; da Mata, G.F.; Passos, M.T.; Varela, P.; Mastroianni-Kirsztajn, G.; Pesquero, J.B. Cathepsin L in COVID-19: From Pharmacological Evidences to Genetics. Front. Cell. Infect. Microbiol. 2020, 10, 589505. [CrossRef] Vaz de Paula, C.B.; de Azevedo, M.L.V.; Nagashima, S.; Martins, A.P.C.; Malaquias, M.A.S.; Miggiolaro, A.; da Silva Motta Junior, J.; Avelino, G.; do Carmo, L.A.P.; Carstens, L.B.; et al. IL-4/IL-13 remodeling pathway of COVID-19 lung injury. Sci. Rep. 2020, 10, 18689. [CrossRef] Guizani, I.; Fourti, N.; Zidi, W.; Feki, M.; Allal-Elasmi, M. SARS-CoV-2 and pathological matrix remodeling mediators. Inflamm. Res. 2021, 70, 847–858. [CrossRef] [PubMed] Gordon, J.C.; Myers, J.B.; Folta, T.; Shoja, V.; Heath, L.S.; Onufriev, A. H++: A server for estimating pKas and adding missing hydrogens to macromolecules. Nucleic Acids Res. 2005, 33, W368–W371. [CrossRef] [PubMed] OMEGA, 2.5.1.4; OpenEye Scientific Software: Santa Fe, NM, USA, 2013. Hawkins, P.C.; Skillman, A.G.; Warren, G.L.; Ellingson, B.A.; Stahl, M.T. Conformer generation with OMEGA: Algorithm and validation using high quality structures from the Protein Databank and Cambridge Structural Database. J. Chem. Inf. Model. 2010, 50, 572–584. [CrossRef] [PubMed] SZYBKI 1.9.0.3, 1.9.0.3; OpenEye Scientific Software: Santa Fe, NM, USA, 2016. Halgren, T.A. MMFF VI. MMFF94s option for energy minimization studies. J. Comput. Chem. 1999, 20, 720–729. [CrossRef] QUACPAC, 1.7.0.2; OpenEye Scientific Software: Santa Fe, NM, USA, 2016. Gasteiger, J.; Marsili, M. Iterative partial equalization of orbital electronegativity—A rapid access to atomic charges. Tetrahedron 1980, 36, 3219–3228. [CrossRef] Heller, S.R.; McNaught, A.; Pletnev, I.; Stein, S.; Tchekhovskoi, D. InChI, the IUPAC International Chemical Identifier. J. Cheminform. 2015, 7, 23. [CrossRef] Morris, G.M.; Huey, R.; Lindstrom, W.; Sanner, M.F.; Belew, R.K.; Goodsell, D.S.; Olson, A.J. AutoDock4 and AutoDockTools4: Automated docking with selective receptor flexibility. J. Comput. Chem. 2009, 30, 2785–2791. [CrossRef] Forli, S.; Huey, R.; Pique, M.E.; Sanner, M.F.; Goodsell, D.S.; Olson, A.J. Computational protein-ligand docking and virtual drug screening with the AutoDock suite. Nat. Protoc. 2016, 11, 905–919. [CrossRef] Mansourian, M.; Fassihi, A.; Saghaie, L.; Madadkar-Sobhani, A.; Mahnam, K.; Abbasi, M. QSAR and docking analysis of A2B adenosine receptor antagonists based on non-xanthine scaffold. Med. Chem. Res. 2014, 24, 394–407. [CrossRef] Case, D.A.; Betz, R.M.; Cerutti, D.S.; Cheatham, T.E.; Darden, T.A.; Duke, R.E.; Giese, T.J.; Gohlke, H.; Goetz, A.W.; Homeyer, N.; et al. AMBER 2016; University of California: San Francisco, CA, USA, 2016. Wang, J.; Wolf, R.M.; Caldwell, J.W.; Kollman, P.A.; Case, D.A. Development and testing of a general amber force field. J. Comput. Chem. 2004, 25, 1157–1174. [CrossRef] Maier, J.A.; Martinez, C.; Kasavajhala, K.; Wickstrom, L.; Hauser, K.E.; Simmerling, C. ff14SB: Improving the accuracy of protein side chain and backbone parameters from ff99SB. J. Chem. Theory Comput. 2015, 11, 3696–3713. [CrossRef] [PubMed] Bayly, C.I.; Cieplak, P.; Cornell, W.D.; Kollman, P.A. A Well-Behaved Electrostatic Potential Based Method Using Charge Restraints for Deriving Atomic Charges—The Resp Model. J. Phys. Chem. 1993, 97, 10269–10280. [CrossRef] Frisch, M.J.; Trucks, G.W.; Schlegel, H.B.; Scuseria, G.E.; Robb, M.A.; Cheeseman, J.R.; Scalmani, G.; Barone, V.; Mennucci, B.; Petersson, G.A.; et al. Gaussian 09, Revision E01, Gaussian Inc.: Wallingford, CT, USA, 2009. Darden, T.; York, D.; Pedersen, L. Particle mesh Ewald: AnN·log(N) method for Ewald sums in large systems. J. Chem. Phys. 1993, 98, 10089–10092. [CrossRef] Izaguirre, J.A.; Catarello, D.P.; Wozniak, J.M.; Skeel, R.D. Langevin stabilization of molecular dynamics. J. Chem. Phys. 2001, 114, 2090–2098. [CrossRef] Berendsen, H.J.C.; Postma, J.P.M.; Vangunsteren, W.F.; Dinola, A.; Haak, J.R. Molecular-dynamics with coupling to an external bath. J. Chem. Phys. 1984, 81, 3684–3690. [CrossRef] Pharmaceuticals 2022, 15, 153 64. 65. 66. 67. 68. 69. 70. 71. 22 of 22 Miyamoto, S.; Kollman, P.A. Settle—An Analytical Version of the Shake and Rattle Algorithm for Rigid Water Models. J. Comput. Chem. 1992, 13, 952–962. [CrossRef] Fish, R.H.; Jaouen, G. Bioorganometallic Chemistry: Structural Diversity of Organometallic Complexes with Bioligands and Molecular Recognition Studies of Several Supramolecular Hosts with Biomolecules, Alkali-Metal Ions, and Organometallic Pharmaceuticals. Organometallics 2003, 22, 2166–2177. [CrossRef] Onufriev, A.; Bashford, D.; Case, D.A. Exploring protein native states and large-scale conformational changes with a modified generalized born model. Proteins 2004, 55, 383–394. [CrossRef] Hou, T.; Wang, J.; Li, Y.; Wang, W. Assessing the performance of the molecular mechanics/Poisson Boltzmann surface area and molecular mechanics/generalized Born surface area methods. II. The accuracy of ranking poses generated from docking. J. Comput. Chem. 2011, 32, 866–877. [CrossRef] Wang, E.; Sun, H.; Wang, J.; Wang, Z.; Liu, H.; Zhang, J.Z.H.; Hou, T. End-Point Binding Free Energy Calculation with MM/PBSA and MM/GBSA: Strategies and Applications in Drug Design. Chem. Rev. 2019, 119, 9478–9508. [CrossRef] Zhao, Y.H.; Abraham, M.H.; Le, J.; Hersey, A.; Luscombe, C.N.; Beck, G.; Sherborne, B.; Cooper, I. Rate-limited steps of human oral absorption and QSAR studies. Pharm. Res. 2002, 19, 1446–1457. [CrossRef] [PubMed] Shannon, P.; Markiel, A.; Ozier, O.; Baliga, N.S.; Wang, J.T.; Ramage, D.; Amin, N.; Schwikowski, B.; Ideker, T. Cytoscape: A software environment for integrated models of biomolecular interaction networks. Genome Res. 2003, 13, 2498–2504. [CrossRef] [PubMed] Blucher, A.S.; McWeeney, S.K.; Stein, L.; Wu, G. Visualization of drug target interactions in the contexts of pathways and networks with ReactomeFIViz. F1000Research 2019, 8, 908. [CrossRef] [PubMed]