Tomatidine and Patchouli Alcohol as Inhibitors of SARS-CoV-2 Enzymes (3CLpro, PLpro and NSP15) by Molecular Docking and Molecular Dynamics Simulations

Considering the current dramatic and fatal situation due to the high spreading of SARS-CoV-2 infection, there is an urgent unmet medical need to identify novel and effective approaches for prevention and treatment of Coronavirus disease (COVID 19) by re-evaluating and repurposing of known drugs. For this, tomatidine and patchouli alcohol have been selected as potential drugs for combating the virus. The hit compounds were subsequently docked into the active site and molecular docking analyses revealed that both drugs can bind the active site of SARS-CoV-2 3CLpro, PLpro, NSP15, COX-2 and PLA2 targets with a number of important binding interactions. To further validate the interactions of promising compound tomatidine, Molecular dynamics study of 100 ns was carried out towards 3CLpro, NSP15 and COX-2. This indicated that the protein-ligand complex was stable throughout the simulation period, and minimal backbone fluctuations have ensued in the system. Post dynamic MM-GBSA analysis of molecular dynamics data showed promising mean binding free energy 47.4633 ± 9.28, 51.8064 ± 8.91 and 54.8918 ± 7.55 kcal/mol, respectively. Likewise, in silico ADMET studies of the selected ligands showed excellent pharmacokinetic properties with good absorption, bioavailability and devoid of toxicity. Therefore, patchouli alcohol and especially, tomatidine may provide prospect treatment options against SARS-CoV-2 infection by potentially inhibiting virus duplication though more research is guaranteed and secured.


Introduction
The current contagious pandemic of the novel severe acute respiratory syndrome coronavirus 2 named SARS-CoV-2, is a viral pneumonia and a threat to global public health [1,2]. According to the current situational report (https://www.worldometers.info/ coronavirus/), released on 22 July 2021, the SARS-CoV-2 has already affected over 179 million people across the world, including 396,758 confirmed cases in Kingdom Saudi Arabia with 7691 deaths. The SARS-CoV2 is a single-stranded positive-sense RNA genome of size 29.7 kb [3]. Structurally, it was composed of Spike glycoprotein (S), Envelope protein (E), Membrane protein (M) and Nucleocapsid protein (N). Envelope (E) protein is essential to release the virus, Membrane (M) protein by increasing the membrane curvature, promotes the viral assembly; Nucleocapsid (N) protein is interferon (IFN) antagonistic and supports viral replication. Hemagglutinin esterase (HE) and Helicase (H) proteins, and nonstructural proteins (NSP) which include proteases, papain-like proteases (PL pro or NSP3) and 3C-like protease (M pro or NSP5), and Replicase proteins [4][5][6]. The envelope and nucleocapsid proteins of SARS-CoV-2 are two evolutionarily conserved regions, with sequence identities of 96% and 89.6%, respectively [7]. M pro and PL pro play significant roles in the post-translational modification of the two replicase polyproteins, pp1a and pp1ab [8,9]. M pro display a central role in proteolysis, viral replication processes, transcriptional processes and is essential for the life cycle of the virus [10]. PL pro is essential for replication/transcription complex (RTC) formation and is responsible for release of NSP1, NSP2 and NSP3 from the N-terminal region of pp1a and 1ab [11]. Their inhibition can block these processes and prevent the viral infection. NSP15 protein also acts as an endoribonuclease and preferentially cleaves 3 of uridylates through a ribonuclease A (RNase A)-like mechanism and also facilitates viral replication and transcription [10]. As a part of the body's immune response, inflammation is a crucial parameter that may be taken into account to reduce infection, limit viral replication and transmission. NF-κB-inducing kinase (NIK), Cyclooxygenase-2 (COX-2), phospholipase A2 (PLA2) and interleukin-1 receptor associated kinase 4 (IRAK-4) are important druggable targets involved in SARS-CoV-2 induced inflammatory response and can be used to screen anti-inflammatory molecules. In fact, NIK activates NF-κB2 by promoting proteolytic processing and the generation of NF-κB transcription of the targeted gene, also known to regulates both inflammation-induced and tumor-associated angiogenesis [12]. COX-2 catalyze the biosynthesis of prostaglandins, prostacylins and thromboxanes, from arachidonic acid mostly involved in pathological conditions such as inflammation, pain and fever [11]. PLA2 enzymes are required to increase the level of arachidonic acid for metabolism and biosynthesis of eicosanoid under physiological condition as well as in inflammatory cell activation [11]. Interleukin-1 receptor-associated kinase 4 (IRAK-4) plays a pivotal role in signaling cascades associated with the immune and inflammatory diseases, and may be an effective therapeutic target for various diseases associated with deregulated inflammation [11].
Alongside the development of new, often high-budget drugs, different strategies were applied in parallel with drug repurposing which remains the best approach to rapidly identify antiviral plausible therapeutic agents. Accordingly, computational approaches to identify suitable and adequate therapeutic options in the short-term and loss-cost, to explore possible interactions with viral proteins have been well investigated [13,14]. The molecular docking and dynamic simulation process were largely applied in structure-based drug design owing to their ability to predict, with reliable accuracy which are intended to decipher the mechanism of binding interactions [15,16].
Extensive research on the therapeutic effect of traditional medicinal plant remedies, have been highlighted against various health ailments as well as an alternative medicine [13,[17][18][19][20]. A particular attention was offered a phytocomponents showing scientific evidence for host protection against SARS-CoV-2 infections. Their effectiveness effect is employed in strengthening of host immune system with less adverse side effects, in contrast with synthetically generated chemical substances to cure many diseases.
Tomatidine as a steroidal alkaloid widely found in skins and leaves of eggplant, potatoes and tomatoes [21]. This steroidal was reported to inhibit the replication of Staphylococcus aureus as well as the invasion activity of A549 human lung adenocarcinoma cells and induces HL60 human myeloid leukemia cell apoptosis [22][23][24]. In addition, tomatidine suppresses iNOS and COX-2 expression by blocking the NF-κB pathway in LPS stimulated RAW 264.7 cells [25]. Recently tomatidine was proved to possess anti-inflammatory effects in macrophages [26]. Other study reported that tomatidine exhibited antiviral activity towards chikungunya virus [27]. Patchouli alcohol, as naturally occurring tricyclic sesquiterpene, has been widely demonstrated for its multi-beneficial pharmacological properties such as immunomodulatory, antitumor, anti-inflammatory, antimicrobial, antioxidative, insecticidal, antiatherogenic and antiemetic [28].
Repurposing existing drugs offers the fastest opportunity to identify new indications for existing drugs as a stable solution against coronavirus disease 2019 (COVID-19). Since, an effective drug against SARS-CoV-2 still needs to be discovered, in this study, we recruited an amalgam of docking-based virtual screening, molecular dynamics (MD) simulations and binding-free energy approaches to identify suitable existing drugs for the treatment of COVID-19. We believe that our findings will help pharmaceutical chemists to optimize suitable drugs for the clinical management of COVID-19 patients. Therefore, based on the above facts, in this study, we performed the molecular docking and dynamic simulation of tomatidine and patchouli alcohol in order to discover a novel potential inhibitor to 3CL pro , PL pro and NSP15 as well as pro-inflammatory mediators such as COX-2 and PLA-2 enzymes. Additionally, pharmacokinetics, drug-likeness and target prediction were assessed.

Molecular Docking Analysis
Patchouli alcohol is one of the important compounds isolated from the essential oil of Pogostemon cablin (patchouli), in parallel with tomatidine as a glycoalkaloid and an agly-cone metabolite of tomatidine which is largely present in high amount in unripe tomato and were tested here, for the first time for their potential effect as anti-SARS-CoV-2, using in silico methods. Proteases 3CLpro and PLpro were targeted to virtually identify the potential drug specifically blocking its catalytic dyad and triad, respectively. Their inhibition plays a pivotal role either in stopping the production of functional SARS-CoV-2 polypeptides or in reducing the viral replication on the host.
The docking method involves searching for binding sites on the whole macromolecular surface. Therefore, (SARS-CoV-2) inhibitors were used for blind docking analysis of the known drugs, tomatidine and patchouli alcohol. The protein-ligand interactions of the stable docked complexes were investigated. According to the docking score and binding mode analyses, both leads were found to tightly bind to the different selected targets.

Docking Analysis of Tomatidine and Patchouli Alcohol with Target 3CLpro
As shown (Figure 1), tomatidine interacts to 3CLpro (PDB ID: 6LU7) via van der Waals interactions with Thr25, Leu27, His41, Met49, Leu141, Asn142, Gly143, Ser144, His164, Glu166, Leu167, Arg188, Gln189 and Thr190, unfavorable donor-donor interactions with Thr26 and Alkyl interactions with Cys145, Met165 and Pro168, amino-acids, however patchouli alcohol established van der Waals interactions with His41, Met49, Phe140, Leu141, Asn142, Gly143, Ser144, Cys145, His164, Glu166 and Gln189, and Alkyl interactions with His163, Met165 and His172 residues. In such a case, all the residues involved in the interaction sites of the tomatidine-3CLpro complex, and patchouli alcohol-3CLpro are located in the preferred high-volume pocket of 3CLpro having Thr24, Thr25, Thr26, Leu27, His41, Cys44, Thr45, Ser46, Met49, Pro52, Tyr54, Phe140, Leu141, Asn142, Gly143, Ser144, Cys145, His163, His164, Met165, Glu166, Leu167, Pro168, His172, Asp187, Arg188, Gln189, Thr190 and Gln192 residues. We note also the presence of both Cys145 and His41 catalytic dyad in tomatidine-3CLpro complex, and patchouli alcohol-3CLpro active center known for their major role in substrate binding and the enzyme activity. Regarding the binding mode of tomatidine and patchouli alcohol with target 3Clpro, our results are in good agreement with what has been proven for SARS-CoV-2 NSP15 indicating that His235, Gln245, Gly248, Gln294 and Thr341 were the main residues in NSP15 active site. We note also that both tomatidine and patchouli alcohol do interact with the important catalytic residue Lys290 via van der Waals and H-bonds interactions, respectively, which is involved preferentially in the hydrolysis of a nucleoside as observed from SARS and MERS NSP15 [7]. Interestingly, the residue His235 which is pivotal for hydrolysis was found in. Additionally, hydroxyl groups of Ser294 and Tyr343 residues which are close to the catalytic residues and may be responsible for the binding of the substrate also make hydrogen bonds with ligands in the models. Moreover, in the case of MERS (NSP15), the importance of the residue Tyr343 for ribonuclease activity has been demonstrated by mutational studies [29].

Molecular Dynamic Simulation and Post Dynamic MMGBSA Binding Free Energy Analysis
Since molecular docking studies were carried out using the protein rigid crystal structure, we have studied target receptor and lead compound interactions in the dynamic behavior of both receptor and ligand using molecular dynamic simulation to investigate the stability of bound conformation after binding of lead compound within the binding cavity of SARS CoV-2 3CLpro, NSP15 Endoribonuclease and human COX-2. To study the conformational stability of the Tomatidine-protein complex and its changes, we have simulated the systems up to 100 ns. In this study, we employed a period of 100 ns, which is sufficient time for the configurations of Cα atoms in complex with lead compound in 3CLpro, NSP15 and human COX-2.
The Root Mean Square Deviation (RMSD) and Root Mean Square Fluctuation (RMSF) values of the Cα atoms of proteins were calculated to indicate thermodynamic conformational stability during a 100 ns time period. The RMSD analysis indicates the simulation's general status and whether it has equilibrated or not. It is believed that the lower the RMSD value throughout the simulation, shows higher stability of the protein-ligand complex. Conversely, greater the RMSD value, reveals the less stable the protein-ligand complex. For protein biomolecules, fluctuations in the lower RMSD range are perfectly acceptable. However, fluctuations in the wide range may indicate that the protein is likely to undergo substantial or significant conformational change during the simulation [8,9].
The RMSD plot of the tomatidine-3CLpro complex is displayed in Figure 6A. For the RMSD study, the backbone of the protein structures enumerated throughout the MD simulation was aligned to the initial structure. After the initial fluctuation due to the equilibration, the RMSD of the tomatidine in complex with SARS CoV-2 3CLpro system gradually increased for 18ns, after which the RMSD remained between 1.2 Å to 2.5 Å till the end of the simulation. In this complex highest fluctuation up to 5.8 Å was observed at a 16-18 ns time span. Targeted protein 3CLpro revealed a maximum RMSD value of 2.7 Å, which indicates that the tomatidine-3CLpro complex was maintained constantly during the simulation time. In SARS CoV-2 NSP15-tomatidine complex, for the initial phase of the simulation showed higher fluctuations up to 5 Å because of the equilibration. For the initial 61 ns of the simulation, the tomatidine RMSD plot gradually increases, later 40 ns the tomatidine fluctuations remained within the range of 2.0 Å, indicating stabilization in protein complex structure. A similar RMSD pattern was observed with NSP15 protein as portrayed in Figure 7A.  The tomatidine in complex with human COX-2 exhibited an average RMSD and RMSF value of 5.19 Å and 1.97 Å, respectively. Initially, the RMSD of the ligand increases up to 13 ns, after which a promising result was observed until 100 ns, having a constant RMSD value of 6.0 Å. In the RMSD plot of protein, initially RMSD increases up to a 23 ns time span, after that minor RMSD variation was noticed that suggests lead compound tomatidine bound tightly within the cavity of human COX-2 ( Figure 8A).
The binding efficiency of lead compound tomatidine was also examined using RMSF values for Cα atoms of all the protein residues based on 100 ns trajectory data. The RMSF value reveals the mobility and flexibility of each amino acid in a protein throughout the simulation time. Larger RMSF values indicate more flexibility during simulation, while lower RMSF values indicate good system stability. In this plot, secondary structural elements, such as the α-helical and β-strand regions, are shown in red and blue backgrounds, respectively, while the loop area is shown in a white background. α-helical and β-strands are often stiffer than the unstructured portion of the protein and so fluctuate less than loop regions. If the active site and main chain atoms fluctuated little, it showed that the conformational change was minimal, suggesting that the reported lead compound was snugly bound within the cavity of the target protein binding pocket 5,6.  Figure 1B). The residues Thr26, Thr25, Val 42 and Gln166, exhibited hydrogen bond interactions with the tomatidine. The Met165, Leu167 and Pro168 showed hydrophobic interactions during the MD simulation. The larger majority of the interactions observed between the 3CLpro and tomatidine during the docking were retained during the MD simulation, indicating a stable binding mode prediction during the docking ( Figure 6C,D). The RMSF plot of tomatidine -NSP15complex is displayed in Figure 2B. The fluctuation of amino acid residues during the interaction was observed to be below 2.4 Å, which is perfectly acceptable. The highest fluctuation was identified with Val36, Asp273, Val39, Ser242, Ser242, Ile270, Asp37, Pro271, Gln347, Ser244, Gly38, His243 and Met272, which did not interact with tomatidine. The terminal hydroxyl group of tomatidine established three hydrogen bonds with Thr275, Lys71 and Asp297. The active site amino acid residues of the SARS CoV-2 NSP15 contribute to binding interactions of 99%, 43% and 40%, respectively. In tomatidine-NSP15 complex, strong hydrogen bonding was noticed with amino acid Lys71, Gln202, Thr275 and Asp297 ( Figure 7C,D).
The Radius of gyration (rGyr) property was also examined to illustrate the stability of the tomatidine in the SARS CoV-2 3CLpro, NSP15 and human COX-2 binding pockets during the simulation of 100 ns (Figure 9). The rGyr parameter is used to measure how extended a ligand is, and is equivalent to its principal moment of inertia. The tomatidine in complexes with 3CLpro, NSP15 and human COX-2 exhibited an average rGyr value of 4.76 ± 0.036 Å, 4.77± 0.034 Å and 4.77± 0.034 Å, respectively (Table 1). No major fluctuation was observed in the rGyr. These constant values showed steady behavior. The post dynamic MM-GBSA analysis of free binding energy (∆G Bind) calculation was carried out with the generation of 0-1001 frames having a 101-step sampling size. A total of 10 frames were processed and analyzed throughout the post dynamic MM-GBSA calculation of 100 ns MD data of tomatidine in complex with the SARS CoV-2 3CLpro, NSP15 and human COX-2 revealed by the dynamic's studies. The computed post dynamic-MMGBSA based binding free energy for the protein ligand complexes is depicted in Table  1. The calculated average ∆G Bind of the complex tomatidine in complex with the SARS CoV-2 3CLpro, NSP15 and human COX-2 was found to be −47.4633 ± 9.28 kcal/mol, −51.8064 ± 8.91 kcal/mol and −54.8918 ± 7.55 kcal/mol, respectively. A more negative value shows stronger binding.

Pharmacokinetics, Drug-Likeliness and Toxicity Studies
A good drug must meet the conditions of absorption, distribution, metabolism and excretion (ADME) as well as toxicity and should entirely and quickly absorb from the gastrointestinal tract, specifically distributed to its target, metabolized in a manner that helps its activity, and excreted without any harm. Based on the predicted results summarized in Table 2, drug-likeness properties indicate that both tomatidine and patchouli alcohol obeyed to Lipinski's, Veber's and Egan's rule assessing their flexibility as well as their surface area, with bioavailability score of 0.55 and consensus log Po/w of 4.90 and 3.57, respectively. Lower hydrophobicity results in a decrease in the metabolism of compounds and high absorption. Pharmacokinetic data showed that both leads were permeators of the blood brain barrier and skin with high gastrointestinal absorption meaning their high absorbance potential. It is well known that cytochrome P450s can regulate the metabolism of various drugs. Tomatidine and patchouli alcohol were found to be no inhibitors of all cytochrome P450 isoforms CYP1A2, CYP2C19, CYP2C9 (except for PA), CYP2D6 and CYP3A4 which are the main players in Phase I metabolism. CYP3A4 as major enzyme system known to acts on lipophilic substrates and recognized to be responsible for the metabolism of about 60% of xenobiotics including drugs, carcinogens, steroids and eicosanoids. Both the leads exhibited negative results on Ames mutagenesis and therefore cannot be considered as a mutagenic agent. In addition, no toxicity was shown with the human ether-a-go-gorelated gene (hERG I) inhibitor, and no hepatotoxicity and skin sensitization, which allows them to be good safety drugs.
The bioavailability radar ( Table 2) of examined drugs represented by the pink area regrouping the optimal range for the followings (lipophilicity: XLOGP3 between −0.7 and +5.0, size: MW between 150 and 500 g/mol, polarity: TPSA between 20 and 130 Å, solubility: log S not higher than 6, saturation: fraction of carbons in the sp3 hybridization not less than 0.25, and flexibility: no more than 9 rotatable bonds with the colored zone defined the desired physicochemical space for good oral bioavailability suggesting that they fall entirely the pink area and therefore exhibiting good drug-likeness properties.

Target Prediction
To perform the molecular mechanisms underlying a given phenotype or bioactivity, and to rationalize possible side-effects, molecular target studies were predicted based on their resemblance with known drugs, to estimate their targets. The top 50 results of the closely associated receptors based on Target, Common Name, Uniprot ID, ChEMBL-ID, Target Class, Probability and Known actives in 2D/3D were illustrated as a pie-chart ( Figure 10). As can be seen, tomatidine predicts 42% family AG coupled receptor, 16% enzyme, 16% kinase and 2% protease while, patchouli alcohol has 16% enzyme and 8% protease, as targets. Overall, the obtained results highlighted the effectiveness of both tomatidine and patchouli alcohol as anti-SARS-CoV-2 agents targeting papain-like proteases (PLpro), 3Clike protease (Mpro) and NSP15 proteins. In fact, tomatidine has been demonstrated to inhibit the production of new Dengue virus particles [31]. More recently, Troost and colleagues [27] reported that this akkaloid exhibited antiviral activity against chikungunya virus by decreasing the number of infected cells and acted at a post-entry step of the virus replication cycle. It has also been elucidated that tomatidine inhibits porcine epidemic diarrhea virus (PEDV) by blocking 3CL protease activity in infected cells, and consequently PEDV replication [32]. Additionally, this molecule is active against Sunnhemp Rossette virus and Tobacco mosaic, and do not affect the replication of herpes simplex virus, human respiratory syncytial and influenza virus [33][34][35][36]. In 2021, Vergoten and Baily [37] studied the in silico interaction of four compounds including tomatidine, with Coronavirus 3Clike protease of PEDV (3CLpro) and (SARS-CoV-2)-main protease (Mpro). They reported that the calculated potential energy of interaction (∆E) and free energy of hydration (∆G) for the interaction of the two viral proteases (3CLpro and Mpro) with tomatidine were: −48.3, −21.6, −44.3, −19.1 (kcal/mol), respectively. Their results highlighted the potential inhibition of PEDV and SARS-CoV-2 main proteases by tomatidine.
Similarly, patchouli alcohol was reported to possess antiviral activity against influenza virus [38,39] with an IC 50% about 2.635 µM. This tricyclic sesquiterpene was active against influenza A (H2N2) virus and especially after penetration of the virus into the cell, and strongly bind to neuraminidase protein with an interaction energy of −40.38 kcal mol −1 [40].

Molecular Docking
Three-dimensional structure of the target proteins was retrieved from RCSB Protein Data Bank [51] with PDB ID: 6LUC, 6W9C, 6VWW, 5F1A and 4UY1 ligands-tomatidine and patchouli alcohol was selected for docking to identify the binding mode and binding affinity of the ligands in the binding pocket of the target protein [52].
Docking of protein ligands was carried out using AutoDock Vina [53], GOLD [54,55] and LibDock from Discovery Studio Client v20. 1.0.19295 (DS) [Dassault Systemes, BIOVIA Corp., San Diego, CA, USA, v 20.1]. The docking softwares used employ different algorithms to improve binding accuracy. It has been demonstrated that these docking programmes are capable of predicting experimental poses with root mean squared deviations (RMSDs) that range between 1.5 and 2Å. Nonetheless, flexible receptor docking, particularly backbone flexibility in receptors, continues to pose a significant challenge to the currently available docking methods [56].
The binding affinity and gold fitness scores were obtained from AutoDock Vina and GOLD, respectively, for obtaining the best orientation and conformation of the ligands.
To acquire accurate results, all the docking experiments were performed with the default parameters.

Docking Using AutoDock Vina
PDBQT files of receptor protein and ligands were prepared using Graphical User Interface program AutoDock Tools (ADT). The bound ligand was removed, and the grid box was created around the ligand bind site of each of the protein. Dimensions of the grid around the binding site vary with each of the proteins. For 6LU7, the grid box was created with size 50 × 50 × 50 xyz points, grid spacing of 0.375 Å and grid center of x, y and z dimensions of −10.729, 12 While for 4UY1, the grid box was created with size 40 × 40 × 40 xyz points with grid spacing of 0.375 Å and grid center at −0.058, −2.230 and 0.718 in x, y and z coordinates, respectively. Protein and ligands were set to the rigid mode during the docking procedure and a configuration file consisting for protein and ligand information along with grid box properties was prepared for executing docking using AutoDock Vina. The ligand pose with highest binding energy/binding affinity was selected for exploring close intra-molecular interactions with the receptor.

Docking Using GOLD (Genetic Optimization for Ligand Docking)
In GOLD suite, wizard was used for docking protein and ligands with default parameters. The active site was defined by selecting the bound ligand within the protein. 10 solutions for each ligand were obtained by applying default Genetic Algorithm settings. The best ligand was selected based on the highest GoldScore fitness function. The ligand and the protein docked complex was further analyzed for close intra-molecular interactions.
The molecular docking and visualization studies were also carried out with the help of commercially available site-features directed docking (LibDock) program in Discovery Studio. Each of the target protein was prepared and binding spheres was defined.. The docking preferences were set to "high quality", and "Best" Conformation method with maximum conformations of 255 were selected. The ligand pose with highest LibDock Score was selected to form docked complex with the receptor for further analysis.

Intra-Molecular Interactions in Docked Complex
Docked complexes of the target proteins with tomatidine and patchouli alcohol obtained using AutoDock Vina, GOLD and LibDock were further analyzed for intra-molecular interactions using View Interaction tool from Discovery Studio Client v20.1.0.19295. Interacting residues of protein and ligand were visualized in 3D and 2D view.

Molecular Dynamic Simulation and Post Dynamic MMGBSA Binding free Energy Analysis
The molecular dynamics simulation was run using the "Academic version of Desmond" (Desmond 2020-3) to study the change in protein structure within the solvent system [57]. Desmond's System Builder panel was used to design the water-soaked solvated system. For the simulations, the complex was centered in an orthorhombic cubic box with periodic boundary conditions and filled with single-point charge (SPC) water molecules buffered at a distance of a minimum of 10 Å between a protein atom and box edges [58,59]. The system was neutralized by adding counter-ions (Na + and Cl -) at random, and an isosmotic state was maintained by adding 0.15 M NaCl. Using OPLS 2005 force field parameters as the default protocol associated with Desmond, the solvated built system was minimised and relaxed [60][61][62]. MD simulations were performed using an isothermal, isobaric ensemble (NPT) with a temperature of 300 K, pressure of 1 atm and a 200 ps thermostat relaxation period. A total of 100 ns simulations were run, during which 1000 trajectories were recorded. Finally, the Simulation Interaction Diagram (SID) tool was used to examine the MD simulation trajectory [63][64][65]. The interaction energies between the protein and the ligand poses were computed using the MM-GBSA (molecular merchandised generalized-born/surface area) method implemented by Schrodinger [23][24][25]. The average post dynamic binding free energy (∆G bind) based on MM-GBSA was calculated using the thermal_mmgbsa.py script.

Molecular Target Predictions
Molecular target predictions are important to find the phenotypical side effects or potential cross reactivity caused by the action of small biomolecules were obtained by using the web tool (http://www.swisstargetprediction.ch/, accessed on 1 May 2021) and entering the smile formats of the desired drugs to obtain the targets. The prediction concerns the putative targets of the given molecule by utilizing 2D and 3D similarity index with known ligands.

Conclusions
In this study, we used docking-based virtual screening to identify potential inhibitors SARS-CoV-2 3CLpro, PLpro, NSP15, COX-2 and PLA2. Both tomatidine and patchouli alcohol modulate interactions that are required for infection by the SARS-CoV-2 virus. These two molecules possessed good pharmacokinetic and druglikness properties. In addition, Tomatidine was found to be potentially bind sufficiently to 3CLpro, NSP15 and human COX-2. To further validate the stability of a potential inhibitor tomatidine, a MD simulation study was performed for 100 ns. The MD simulation result revealed significant stability of test compound tomatidine in the active site of SARS-CoV-2 3CLpro, NSP15 and human COX-2, and exhibited good interactions with the surrounded amino acid residues. The post dynamic MM-GBSA analysis of the compound tomatidine showed significant binding affinity with its selected target. Moreover, much research is still required for the effectiveness of tomatidine in human clinical trials.

Conflicts of Interest:
The authors declare no conflict of interest.