In silico bioprospecting of taraxerol as a main protease inhibitor of SARS-CoV-2 to develop therapy against COVID-19.

COVID-19 was caused by a novel coronavirus known as SARS-CoV-2. The COVID-19 disease outbreak has been avowed as a global pandemic by the World Health Organization at the end of March 2020. It leads to the global economic crash, resulting in the starvation of a large population belonging to economically backward countries. Hence, the development of an alternative medicine along with the vaccine is of the utmost importance for the management of COVID-19. Therefore, screening of several herbal leads was performed to explore their potential against SARS-CoV-2. Furthermore, viral main protease was selected as a key enzyme for performing the study. Various computational approaches, including molecular docking simulation, were used in the current study to find potential inhibitors of viral main protease from a library of 150 herbal leads. Toxicity and ADME prediction of selected molecules were also analysed by Osiris molecular property explorer software. Molecular dynamic simulation of the top 10 docked herbal leads was analysed for stability using 100 ns. Taraxerol (-10.17 kcal/mol), diosgenin (10.12 kcal/mol), amyrin (-9.56 kcal/mol), and asiaticoside (-9.54 kcal/mol) were among the top four herbal leads with the highest binding affinity with the main protease enzyme. Thus, taraxerol was found to be an effective drug candidate against the main protease enzyme for the management of COVID-19. Furthermore, its clinical effect and safety profile need to be established through an in vivo model. Supplementary information: The online version contains supplementary material available at 10.1007/s11224-022-01943-x.

Introduction

The coronavirus disease (COVID-19) outbreak was caused by SARS-CoV-2 and spread from Wuhan, China, at the end of December 2019 [1, 2]. COVID-19 was thought to have emerged from a seafood market from an unknown species, causing the virus’s emergence and symptoms to be comparable with pneumonia [1]. The confirmation of human-to-human transmission of this contagious virus was made by the National Health Commission of China in late January 2020. As SARS-CoV-2 began to spread across international borders, affecting the populations of multiple countries, the WHO designated it as an internationally concerning public health emergency and declared it a global pandemic in March 2020 [3].

Based on the sequence resemblance of the SARS-CoV-2 nucleotide, it is considered a betacoronavirus, having a high similarity to well-known and aggressive strains of human coronaviruses (HCOVs). The nuclear material of SARS-CoV-2 was found to be a long positive-sense single-stranded RNA encoding two discrete types of structural as well as non-structural viral proteins [1, 2]. The 5′-untranslated region (5′UTR) and a replicase complex encode for viral structural proteins like spike, nucleocapsid, matrix, and envelope, while the 3′UTR and several unidentified open reading frames encode for viral non-structural proteins like protease, phosphatase, and polymerases [1, 2].

COVID-19 has initial symptoms of cough, malaise, fatigue, fever, body ache, loss of smell or taste sensation, inadequacy, and shortness of breath followed by respiratory distress [4]. Persistent viral infection affects not only the lungs but also other key organs of the body, eventually leading to organ failure [5-7]. Coronavirus transmits from person to person through coughing and sneezing, which spreads through the nasal mucosa via airborne droplets, where it replicates narrowly in ciliated epithelial cells, causing cell damage and inflammation. The transmission of this virus was also observed from infected surfaces to an individual’s hands and then to their bodies [8].

Globally, more than 216 countries have been affected by the pandemic outbreak of COVID-19 until December 31, 2021. More than 287,022,026 confirmed cases of coronavirus-affected people have been reported, with more than 5,447,886 deaths. According to the WHO, China is the first country to report 102,083 confirmed cases of COVID-19 with 4636 mortalities. In addition, new cases of COVID-19 are being reported day by day, in the USA, with more than 55 million confirmed cases and 845,745 deaths, followed by India with more than 34 million confirmed cases and Brazil with more than 22 million confirmed cases [9]. Even after 18 months of the outbreak of the novel SARS-CoV-2, there is still no approved therapy or vaccine available for the cure of this lethal infectious disease. In light of the foregoing, it is critical to design a novel remedy to combat this worldwide epidemic crisis.

The role of viral main protease enzyme has been identified in causing infection. It plays an important role in pathogenic entry into the host via ACE2 and could be used as a therapeutic antiviral drug target. Inhibiting the viral main protease restricts the pathogen’s entry into the host cell and terminates the viral infection [4, 10]. As a result, an inhibitor of the main protease is considered to be a successful treatment approach for COVID-19 pandemics all over the world. The proven methodology for in silico screening of herbal leads to identify modest inhibitors of the viral main protease enzyme is depicted graphically in Fig. 1.

Fig. 1

An in silico approach for screening herbal leads for potential inhibitors of the viral main protease enzyme in order to discover new antiviral therapies against SARS-CoV-2

Herbal medicine is a holistic medicine, having being used for the management of several health problems for thousands of years. Herbal medicine inspired several drugs, such as artemisinin, paclitaxel, reserpine, morphine, quinine, emetine, aspirin, and many more, which have been discovered for the treatment of numerous diseases [11]. Most of the world’s population relies on herbal medicine as an alternative and complementary medicine for the management of diseases like COVID-19. In this context, several herbal medicines and their bioactive leads including Camellia sinensis (epigallocatechin gallate), Andrographis paniculata (andrographolide), Artemisia annua (artemisinin), Betula sp. (betulinic acid), Citrus sp. (hesperidin), Curcuma longa (curcumin), Ficus benjamina (biorobin), Glycyrrhiza glabra (glycyrrhizin), Mollugo cerviana (vitexin), Myristica fragrans (myricitrin), Piper nigrum (piperine), Radix sophorae (matrine), Stephania tetrandra (tetrandrine), Tinospora cordifolia (berberine), Torreya nucifera (luteolin), etc. have been explored through computational approaches for the treatment of COVID-19 [12-32]. Hence, it inspires us to find such herbal leads through an in silico computational approach for the management of the current pandemic situation. Thus, in the current investigation, putative antagonists of the main protease enzyme of the SARS-CoV-2 were acknowledged using docking-based computational screening of herbal ligands. Moreover, in order to develop novel therapy to combat COVID-19, validation of herbal leads with respect to time by using molecular dynamics (MD) simulation was done.

Material and methods

Preparation of compound library

A molecular ligand library was prepared by considering the diverse compounds already reported for their antiviral properties [33-38]. The already reported antiviral compounds were used to prepare a ligand library in the current study with the intention of validating their effectivity against the SARS-CoV-2 as well as establishing the probable mechanism of action of those compounds which are supposed to be interacting with the main protease enzyme for their antiviral effect against SARS-CoV-2.

Molecular docking simulation-based virtual screening

The main protease enzyme of the SARS-CoV-2 complex with an antagonist N3 was procured from the protein databank (PDBID: 6lu7) [10, 39, 40]. The parted macromolecular model was set for docking by amputation of superfluous water molecules and the addition of polar hydrogens, while the reference ligand N3 was prepared by assigning rotatable, non-rotatable, and un-rotatable bonds [41, 42]. The interacting residues of the viral enzyme with the complex inhibitor were evaluated by using PyMol software to confirm the macromolecular active site and prepare a grid box [42, 43]. The prepared grid box was further utilised by the Autogrid module to evaluate the chemical configuration of the macromolecule as well as the ligand to create maps for various atoms, necessary to accomplish docking simulation [41, 42].

The confirmational search by Autodock for the execution of the docking process was performed by the Lamarckian Genetic Algorithm (LGA). The force-field calculates the ligand’s binding energy by integrating intramolecular energies and assessing energetics for bound and unbound states using a comprehensive thermodynamic model. Each ligand molecule’s docking parameters is stored in the docking parameter file (DPF) [44, 45].

For validation purposes, the reference ligand N3 was docked against the macromolecular target. Validation of the grid parameters as well as protocols used in the current docking procedure was performed by considering the chemical interactions as well as the conformation of the docked ligand N3 with respect to its crystallised bioactive conformation.

Validated docking protocols were further utilised for computational screening of a ligand library consisting of 150 diverse herbal ligands against the viral main protease. Virtual screening of herbal ligands was performed with the intent of identifying prospective lead compounds targeting the viral main protease enzyme [42].

The interactions of the ligands against the macromolecular residues were taken into account while evaluating the docking outcomes for all ligands screened against the viral main protease enzyme. The compounds with the minimum binding energy were selected as the lead molecules based upon the predefined empirical range for the obtained binding energy of −5 to −15 kcal/mol [46-50].

Prediction of pharmacokinetics and toxicological properties

Evaluation of physicochemical, pharmacokinetics, and toxicological parameters of the lead molecules obtained after the docking-based screening were predicted by using the pkCSM-pharmacokinetics server. The online server screens the ligand molecule for smooth kinetics within the body required for its physiological expression without toxicity [51].

Lead molecules had their molecular weight (MW), partition coefficient (LogP), hydrogen bond acceptor/donor (HBA/HBD), rotatable bonds count, topological polar surface area (TPSA), and solubility (LogS) calculated. These parameters regulate the ADMET of the lead candidates. The permeability across the blood–brain barrier (BBB), gastrointestinal absorption, and the inhibitory potential of the cytochrome P450 isoenzymes were also considered for the leads to determine their druggability. These parameters were standardised with the Veber’s rule, the Lipinski’s rules of 5, and druggability [52, 53].

Molecular dynamic simulation

Amyrin, diosgenin, and taraxerol were further selected as potential leads of the viral main protease for further simulation to study the stability of their macromolecular complex with respect to time based on their docking score, observed chemical interactions against the target receptor, pharmacokinetic evaluation, as well as toxicity profiling. MD simulation for a shorter duration of timeframe, lasting 10 ns, was performed initially for all the three shortlisted lead molecules. Based on the results obtained for the 10-ns MD simulation, taraxerol was selected to further magnify the MD simulation for a longer time period of 100 ns.

The protein–ligand macromolecular complex was simulated for 100 ns by using Schrodinger’s Desmond module at a constant temperature of 300 K [54, 55]. MD simulations were performed by solvating the macromolecular complex in an explicit water box of size 10 and modelling the protein–ligand macromolecular complex using the OPLS3e force field [56-58]. Previously, promising and repeatable results were reported using the OPLS3e force field and the SPC water model for the protein targets complexed with the organic ligands [59]. By introducing the ions and subsequently minimising their energy, the macromolecular complex was neutralised. The Nose–Hoover thermostatic algorithm [60] was used to keep the energy-minimised macromolecular complex at 300 K, while the Martina–Tobias–Klein approach was used to keep the pressure constant throughout the simulations [61]. The long-range electrostatic interactions between the ligand and the macromolecule were calculated using the particle mesh Ewald (PME) [62] technique with 0.8 grid spacing and a cutoff radius of 9.0 for Coulomb interactions after the NPT-ensemble MD simulation was performed for 100 ns. The simulation interaction diagram tool in the Desmond was used to examine the ligand’s precise binding interactions with the macromolecular main protease [54, 55].

The root mean square deviation (RMSD) of the atoms of both receptors and the complexed ligand were calculated by comparing their reference frame to determine atomic displacement for a certain timeframe during the complexation process. During the simulation phase, the macromolecular residues’ root mean square fluctuation (RMSF) was calculated in reference to their initial condition in the crystalline structure. Throughout the simulation technique, the residue index was used to plot the distribution of macromolecular secondary structural elements (SSE) such as α-helices and β-strands. The binding interactions of the ligand within the macromolecular site were separated into four groups throughout the simulation procedure: hydrogen bonds, hydrophobic interactions, ionic interactions, and water interaction bridges. The radius of gyration (rGyR), molecular surface area (MolSA), solvent accessible surface area (SASA), and polar surface area of the ligands were all measured (PSA) [63-65]. In comparison to its initial frame, the RMSD value of the ligand molecule was determined throughout the whole simulation duration. The PSA was determined by taking into consideration the total contribution of oxygen and nitrogen atoms, whereas the ligand’s extended length, which correlates to its main moment of inertia, was computed using a 1.4 probe radius.

Results

A three-dimensional structural model of the viral macromolecule with a resolution of 2.16 was obtained using the X-ray diffraction technique. Viral protein is made up of a single chain of 306 amino acid residues. Complex ligand N3 comprises 26 rotatable bonds and 9 aromatic carbons, making it a complex inhibitor. In the current computational investigation, all bonds of the inhibitor molecule were kept flexible, and it was also stored in the AutoDock software’s pdbqt format.

The ligand N3 and active binding residues of the viral main protease were wrapped together to prepare a three-dimensional grid box. The exact grid coordinates were obtained from our previously published study on the same protein [46–50, 63, 66]. The docking results for the complex inhibitor N3 against viral macromolecular targets are tabulated in Table 1.

Table 1

Docking results of ligand N3 against the viral main protease enzyme

Macromolecule	RMSD (Å)	Binding energy (kcal/mol)	Ki (nM)
6lu7	0.88	−8.22	935.62

The lead molecules were chosen based on their affinity for the viral main protease enzyme. The lowest binding energy obtained for preeminent pose for each ligand and their binding interactions with the macromolecule were used to determine the binding affinity of the leads. Table 2 shows the binding energy acquired from docking-based virtual screening of the top 10 shortlisted herbal leads. The interacting residues were explored by analysing the binding site of the target macromolecule with the help of Discovery Studio visualizer software with the intent of identifying the important amino acid residues of the viral main protease enzyme.

Table 2

Binding energy of short-listed herbal leads for the viral main protease enzyme

The movement of drugs within the human body is regulated by physicochemical and pharmacokinetic parameters. These regulatory parameters are crucial for the pharmacological expression via pharmacodynamics as well as toxicological properties of a drug molecule. The pkCSM webserver was used for prediction of physicochemical, pharmacokinetic, as well as toxicological parameters of the lead molecules [51]. The physicochemical, pharmacokinetics, and pharmacodynamics properties of the selected lead molecules are tabulated in Table 3.

Table 3

Physicochemical, pharmacokinetics, and pharmacodynamics properties of the selected lead molecules for the SARS-CoV-2 main protease enzyme

Property	Descriptor	Taraxerol	Diosgenin	Amyrin	Asiaticoside	Momordicin	Hecogenin	Guggulsterone	Andrographolide	Pelargonidin	Lupeol
MW	(g/mol)	412	414.63	426.73	959.13	472.71	430.63	312.45	350.45	271.25	426.73
LogP	-	7.78	5.71	8.025	−1.0328	5.456	4.973	4.644	1.963	3.203	8.025
Rotatable bond	-	0	0	0	9	5	0	0	3	1	1
HBA	-	1	3	1	19	4	4	2	5	4	1
HBD	-	1	1	1	12	3	1	0	3	4	1
TPSA	(Å)2	20.93	38.69	20.23	315.2	77.76	55.76	34.14	86.99	90.15	20.23
Absorption	Water solubility (mol/L)	−6.628	−5.684	−6.906	−2.797	−4.759	−5.453	−5.401	−3.045	−3.393	−5.958
Absorption	Caco2 permeability	1.299	1.227	1.358	−0.637	0.951	1.186	1.477	1.213	0.37	1.356
Absorption	Intestinal absorption (%) (human)	97.705	96.491	97.891	5.812	97.953	95.816	100	94.961	85.871	99.707
Absorption	Skin permeability (Log Kp)	−2.706	−3.449	−2.828	−2.735	−2.546	−3.726	−2.768	−3.158	−2.735	−2.675
Absorption	P-glycoprotein substrate	No	No	No	Yes	Yes	No	No	Yes	Yes	Yes
Absorption	P-glycoprotein I inhibitor	Yes	Yes	Yes	Yes	Yes	Yes	Yes	No	No	Yes
Absorption	P-glycoprotein II inhibitor	Yes	Yes	Yes	No	Yes	Yes	No	No	No	Yes
Distribution	VDss (human)	0.157	0.52	0.321	−0.364	−0.277	0.272	0.343	−0.289	0.17	−0.177
Distribution	Fraction unbound (human)	0	0	0	0.373	0.021	0.001	0	0.225	0.125	0
Distribution	BBB permeability	0.724	0.298	0.713	−2.675	−0.467	−0.17	0.19	−0.766	−1.221	0.721
Distribution	CNS permeability	−1.419	−2.903	−1.925	−5.015	−1.847	−1.484	−2.021	−2.484	−2.088	−1.261
Metabolism	CYP2D6 substrate	No	No	No	No	No	No	No	No	No	No
Metabolism	CYP3A4 substrate	Yes	Yes	Yes	No	No	Yes	Yes	Yes	No	Yes
Metabolism	CYP1A2 inhibitor	No	No	No	No	No	No	No	No	Yes	No
Metabolism	CYP2C19 inhibitor	No	No	No	No	No	No	No	No	Yes	No
Metabolism	CYP2C9 inhibitor	No	No	No	No	No	No	No	No	Yes	No
Metabolism	CYP2D6 inhibitor	No	No	No	No	No	No	No	No	No	No
Metabolism	CYP3A4 inhibitor	No	No	No	No	No	No	Yes	No	Yes	No
Excretion	Total clearance (log ml/min/kg)	0.004	0.328	0.119	0.286	0.457	0.311	0.61	1.183	0.692	0.153
Excretion	Renal OCT2 substrate	No	Yes	No	No	No	Yes	No	No	No	No
Toxicity	AMES toxicity	No	No	No	No	No	No	No	No	No	No
Toxicity	Max. tolerated dose (human) (log mg/kg/day)	−0.194	−0.234	−0.275	−0.847	−1.438	−0.588	−0.355	−1.1	0.787	−0.263
Toxicity	hERG I inhibitor	No	No	No	No	No	No	No	No	No	No
Toxicity	hERG II inhibitor	Yes	No	Yes	Yes	No	No	Yes	No	No	No
Toxicity	Oral rat acute toxicity (LD50) (mol/kg)	2.861	1.812	2.345	2.516	4.003	1.927	1.801	2.786	2.436	3.215
Toxicity	Oral rat chronic toxicity (LOAEL) (mg/kg/day)	1.144	1.562	1.031	8.049	1.765	1.294	1.78	0.389	1.595	1.688
Toxicity	Hepatotoxicity	No	No	No	No	Yes	No	Yes	No	No	No
Toxicity	Skin sensitisation	No	No	No	No	No	No	No	No	No	No
Toxicity	T. Pyriformis toxicity (mg/L)	0.394	0.397	0.418	0.285	0.364	0.376	1.014	0.346	0.435	0.335
Toxicity	Minnow toxicity	−1.552	−0.229	−1.953	10.986	0.095	0.563	−0.028	1.157	0.957	−2.002

Taraxerol has all the parameters well within the predefined range according to Lipinski’s rule of five, i.e., MW < 500, < 10 HBA, < 5 HBD, TPSA 20–130 Å2, except LogP value, which exceeds the limit of 5.0. The demonstrated physicochemical properties suggest that taraxerol has optimised pharmacokinetics for drug-like candidature [51]. Taraxerol has not displayed a substrate-like expression for P-glycoprotein, which acts as a biological barrier for xenobiotics and exogenous toxins.

Taraxerol has not shown substrate-like expression for most of the cytochrome P450 isoenzymes, except for the metabolic CYP3A4. The ligand consistently exhibits a high proclivity for excretion from biological systems and low expression across a variety of toxicity pathways such as AMES, hERG I inhibition, hepatotoxicity, skin sensitization, Tetrahymena pyriformis toxicity, and minnow toxicity. Cumulatively, the predicted physicochemical, pharmacokinetics, pharmacodynamics, and toxicological parameters for taraxerol fall well within the predefined range, clearly indicating its promising potential as a druggable candidate [51-53].

The MD simulation of the viral main protease enzyme for a shorter timeframe of 10 ns found taraxerol to be the most stabilised inhibitor molecule complex in the active cavity of the viral main protease enzyme. The putative herbal inhibitor compound taraxerol was further validated by using the Schrodinger’s Desmond module to run a 100-ns MD simulation. There are 306 residues in the macromolecular receptor, but the ligand atom contains only one rotatable bond and 30 heavy atoms out of a total of 78. Based on structural validation throughout the operation, the RMSD analysis supports the smooth execution of the equilibrated simulation process. The ligand RMSD demonstrates the ligand’s stability in relation to the macromolecular binding residues during the simulation procedure by aligning their heavy metals.

The RMSD value for the macromolecular residues was found to be well within the 2.8 Å range, indicating that the majority of the residues do not shift from their starting position during ligand molecule complexation. The ligand molecule’s RMSD value has been maintained well between 6 and 7.5 Å throughout the simulation run, despite some early swings of up to 9.0 Å during ligand adjustment within the macromolecular binding cavity, indicating the ligand’s adjustment within the macromolecular cavity during the initial 20 ns of the simulation, followed by achieving its stable conformation for the remaining 80-ns duration of the simulation process. The ligand taraxerol undergoes a sequence of vibrations after reaching the active binding site of the viral main protease to achieve the most stable confirmation inside the active binding site. As a result, the early oscillations in the ligand’s RMSD value during 3–10 ns are caused by these continuing vibrations while performing certain manoeuvres to achieve the most stable confirmation inside the active binding site of the viral main protease enzyme. Figure 2 shows the RMSD of the protein and ligand over the 100-ns period of the MD simulation. The RMSF value of macromolecular residues was determined to be well within the allowed range of 3 Å. A few residues changed slightly, with an RMSF value of 2–3 Å. However, the bulk of residues exhibited smaller variations, with an average value of less than 1 Å. Figure 3a, b show the RMSF of the viral main protease and the complex ligand taraxerol during the 100-ns period of the MD simulation. The RMSF observed for the ligand taraxerol complex within the active binding site of the viral main protease enzyme was found to be within the range of 2–3 Å throughout the 100-ns timeframe. This clearly suggests that the ligand was stabilised within the active site of the target macromolecule with slight fluctuations of some functionals required for interacting with the target macromolecule.

Fig. 2

Root mean square deviation (RMSD) of the viral main protease enzyme and ligand taraxerol observed during the 100-ns MD simulation

Fig. 3

Root mean square fluctuation: During the 100-ns timeframe of the MD simulation, the RMSF of the viral main protease enzyme and the complexed ligand taraxerol were measured

During the simulation process, the SSE analysis revealed that it has 14.94% alpha helices and 24.21% beta-sheets, resulting in a total contribution of 39.14% SSE, which may be conserved for the majority of the simulation. The macromolecular ligand interaction study revealed that Pro168, Leu167, Met165, Cys145, Met49, Gln192, Thr190, and Arg188 interact with the ligand over the course of the 100-ns simulation procedure. More than eight macromolecular residues were discovered to be consistently interacting with the complex ligand throughout the simulation run. Figure 4 depicts the detailed protein–ligand connections observed during the whole timeframe of the 100-ns MD simulation. The ligand’s RMSD value was substantially within the range of 1.5 Å, indicating that the ligand fluctuated as little as possible during the simulation process. The ligand’s rGyr value was discovered to be in the range of 4.1–4.3 Å. During the modelling phase, no intramolecular hydrogen bonding in the ligand was observed. Throughout the modelling phase, the ligand’s MolSA was discovered to be in the range of 372–384Å2. After some initial variations, the ligand’s SASA was found to be in the range of 150–300Å2 during the simulation process. The PSA for the complex ligand was discovered to be between 35 and 40 Å2 during the simulation process.

Fig. 4

Protein–ligand contacts: Detailed protein–ligand interactions were found during the 100-ns MD simulation timescale

Discussion

The in silico computational approach is pioneering for drug discovery and development. Several drugs have been identified for the treatment of many diseases, including viral infections. There are several examples of in silico inspired medicines such as chloroquine, hydroxychloroquine, remdesivir, ritonavir, lopinavir, darunavir, arbidol, peptide EK1, cobicistat, amodiaquine, niclosamide, telbivudine, nicotinamide, ribavirin, sofosbuvir, galidesivir, tenofovir, setrobuvir, artemisinin, mefloquine, halofantrine, lumefantrine, amodiaquine, primaquine, doxycycline, atovaquone, sulfonamides didanosine, camptothecin, teicoplanin, velpatasvir, ledipasvir, and favipiravir that have been explored against COVID-19. But their clinical uses are limited due to quality-related safety issues. Most of the drugs are under clinical trials [33].

As we know, the world is facing health and economic problems due to COVID-19. Still, there are no such therapies as some vaccines have been approved for the treatment of COVID-19. There is a continuous quest to find drugs for treating COVID-19. In this context, herbal medicine is very popular throughout the globe for its miraculous effect on health problems. It has been used as an alternative medicine for treating various diseases. As per the WHO, 80 percent of the world’s population depends upon herbal medicine for the cure and treatment of several diseases. Several drugs have been discovered through natural resources, including plants [67]. Therefore, it occupies our mind to identify such potential herbal leads for the management of COVID-19. In the recent past, numerous herbal leads like luteolin, psoralidin, caffeic acid, myricetin, quercetin, vitexin, tryptanthrin, shikonin, silvestrol, scutellarein, biorobin, hesperidin, andrographolide, catechins, reserpine, emetine, artemisinin, quinine, paclitaxel, and morphine have been recognised through computational approaches. They have exhibited remarkable inhibitory activity against SARS-CoV-2. Various research studies at preclinical level (in vitro and in vivo) and clinical trials are going on to ensure their therapeutic efficacy and safety [33].

In the present study, lead molecules with binding energy such as taraxerol (−10.17 kcal/mol), diosgenin (−10.12 kcal/mol), amyrin (−9.56 kcal/mol), asiaticoside (−9.54 kcal/mol), momordicin (−9.51 kcal/mol), hecogenin (−9.42 kcal/mol), guggulsterone (−9.23 kcal/mol), andrographolide (−8.61 kcal/mol), pelargonidin (−8.49 kcal/mol), and lupeol (−8.48 kcal/mol) were identified for inhibiting the viral main protease enzyme. Based on the observed chemical interactions, which are tabulated in Table 2, among the shortlisted leads against the viral macromolecular targets, a SAR was developed that confirms the involvement of macromolecular residues such as Cys145, Pro168, Met165, Leu167, Gln192, Thr190, and Met49 in the active ligand binding. These leads have been reported in medicinal plants, including Clitoria ternatea (taraxerol), Dioscorea sp. (diosgenin), Eclipta alba (amyrin), Centella asiatica (asiaticoside), Momordica charantia (momordicin), Chlorophytum borivilianum (hecogenin), Commiphora sp. (guggulsterone), Andrographis paniculata (andrographolide), Anagallis monelli (pelargonidin), and Carissa spinarum (lupeol) [11]. As a result, these bioactive-containing medicinal plants would be effective in managing the SARS-CoV-2 health issue. In another study, Kar et al. reported the inhibitory potential of Clerodendrum trichotomum containing lead taraxerol, friedelin, and stigmasterol against SARS-CoV-2 through an in silico approach [68].

Taraxerol, among other things, inhibited SARS-CoV-2 enzymes (spike protein (6lzg) [−7.5 0.01 kcal/mol]), main protease (8.4 0.01 kcal/mol), and RNA-dependent RNA polymerase (RdRp) [−7.4 0.02 kcal/mol] [69]. Taraxerol has been reported for several biological activities, including antiplasmodial (Plasmodium falciparum, IC50, 8.5 µM), antiparasitic (Trypanosoma brucei, IC50, 10.5 µM), antioxidant (IC50, 500 µM), etc. [70].

Similarly, Enmozhi et al. have reported the antiviral activity of andrographolide on SARS-CoV-2 main protease enzyme through in silico investigation. This investigation includes molecular docking, target analysis, toxicity prediction, and ADME prediction for andrographolide. Andrographolide, a bioactive component of A. paniculata, showed an inhibitory effect with binding efficiency (−3.094357 kcal/mol) against the main protease enzyme [13, 71]. Moreover, other lead compounds from plants (epicatechin gallate −7.24 kcal/mol, catechin −7.05 kcal/mol, kaempferol −9.41 kcal/mol, gingerol −6.67 kcal/mol, quercetin −8.58 kcal/mol, curcumin −7.31 kcal/mol, and demethoxycurcumin −8.17 kcal/mol) have shown inhibitory effects against the main protease enzyme [72]. In addition, phytoconstituents like boldine, ginkgetin, isoboldine, laurolitsine, lauroscholtzine, laurotetanine, norisoboldine, pseudolycorine, secoboldine, and syringic acid have shown good inhibitory potential against the main protease enzyme through in silico studies. The best binding affinity has been produced by laurolitsine as compared to other molecules [73]. Further, efficacy, drug/herb-drug interaction, safety, and toxicity profile of taraxerol need to be validated against SARS-CoV-2 through preclinical and clinical investigation. Hence, an in silico computational approach could be promising for getting new herbal leads for the management of COVID-19. Other alternate treatment options such as stem cells, monoclonal antibodies, plasma, and polypeptides are also under investigation for COVID-19.

Conclusion

The in silico virtual screening technique is a highly adequate, prudent, and quick approach for identifying a potent herbal lead having therapeutic activity against the main protease enzyme of COVID-19. Taraxerol was found to be a potential lead against the main protease enzyme of SARS-CoV-2. The binding site analysis has revealed that the residues Cys145, Pro168, Met165, Leu167, Gln192, Thr190, and Met49 play an important interacting role in the ligand binding. The results obtained by molecular docking simulation were further validated by performing 100-ns molecular dynamics simulation, and the pharmacokinetic profiling based on their physicochemical parameters was also performed to strengthen the candidature of taraxerol as a safe and optimised inhibitor of the viral main protease enzyme. Further, this herbal molecule needs to be validated through preclinical and clinical investigation for its therapeutic applicability.

Below is the link to the electronic supplementary material.

Supplementary file1 (DOCX 1597 KB)