Withasomniferol C, a new potential SARS-CoV-2 main protease inhibitor from the Withania somnifera plant proposed by in silico approaches

Exploring potent herbal medicine candidates is a promising strategy for combating a pandemic in the present global health crisis. In Ayurveda (a traditional medicine system in India), Withania somnifera (WS) is one of the most important herbs and it has been used for millennia as Rasayana (a type of juice) for its wide-ranging health benefits. WS phytocompounds display a broad spectrum of biological activities (such as antioxidant, anticancer and antimicrobial) modulate detoxifying enzymes, and enhance immunity. Inspired by the numerous biological actions of WS phytocompounds, the present investigation explored the potential of the WS phytocompounds against the SARS-CoV-2 main protease (3CLpro). We selected 11 specific withanolide compounds, such as withaphysalin, withasomniferol, and withafastuosin, through manual literature curation against 3CLpro. A molecular similarity analysis showed their similarity with compounds that have an established inhibitory activity against the SARS-CoV-2. In silico molecular docking and molecular dynamics simulations elucidated withasomniferol C (WS11) as a potential candidate against SARS-CoV-2 3CLpro. Additionally, the present work also presents a new method of validating docking poses using the AlteQ method.


INTRODUCTION
The coronavirus disease 2019  caused by the new coronavirus represents a significant global health crisis and has posed unprecedented challenges since December 2019, when the first case was reported in Wuhan, China. Later, the new coronavirus was named SARS-CoV-2 (Severe Acute Respiratory Syndrome Coronavirus 2), as it shares 79.5% genome similarity with the other severe acute respiratory syndrome (SARS) virus (Zhou et al., 2020;Lu et al., 2020). Still, the number of COVID-19 patients is increasing, (residues 1-7) is not involved in the dimerization (Zhong et al., 2008) but is required for the activity of the enzyme. In the dimer, the N-finger of the inactive protomer interacts with domains II and III of the active protomer (Yang et al., 2003;Chen et al., 2005;Chen et al., 2006). Recently, different variants of the SARS-CoV-2 main protease were identified from clinical samples, and their structural variations were explored through in silico studies, with the conclusion that mutations in these variants do not alter the active site conformation (Martin et al., 2020). Hence, the 3CL pro is a good target for drug discovery. So far, peptidomimetic alpha ketoamide inhibitors, the Michael acceptor N3 inhibitor, carmofur, ebselen, aldehyde-based compounds, lopinavir/ritonavir, the antiplatelet drug dipyridamole, boceprevir, GC-376, calpain inhibitors (II, XII), and GC-373 are reported as the most promising drugs against the SARS-CoV-2 3CL pro , and they exert their inhibitory effect by binding to the substrate-binding cleft. However, their clinical outcomes in humans are still not determined (Mengist et al., 2021).
For decades, natural products and herbal medicines have been used for combating numerous viral infections, with their favorable efficacy and low toxicity making them a promising resource for drug discovery. An acceptable toxicity of natural products and herbal medicines make them prospective candidates against COVID-19 (Komolafe et al., 2021). Recent reviews show the significance of natural products against COVID-19 by examining their randomized controlled trial (RCT) reports in COVID-19 patients (Feng et al., 2021;Di et al., 2021). In this regard, exploring potential herbal medicine candidates is a promising strategy for combating a pandemic in the present global health crisis. Withania somnifera (Solanaceae, WS), popularly known as 'Ashwagandha' and 'Indian Ginseng' is a medicinal plant used as a herbal tonic to treat various kinds of diseases, such as cancer, arthritis, asthma, aging, inflammation, and neurological disorders in Indian traditional medicine (Dar, Hamid & Ahmad, 2015). Its pharmacological activities are mainly due to the presence of diverse secondary metabolites, such as alkaloids, flavanol glycosides, glycowithanolides, steroidal lactones (withanolides), sterols, and phenolic acids. Recently, a clinical trial has been initiated by the AYUSH ministry of India for the use of WS along with other plants to evaluate its efficacy against COVID-19 (CTRI (Clinical Trial Registry -India), registration number: CTRI/2021/06/034496, date of registration: June 30, 2021) (Chopra et al., 2021). In addition, the Indian Government has collaborated with the U.K's London School of Hygiene and Tropical Medicine (LSHTM) to conduct a study on ''Ashwagandha'' for promoting COVID-19 recovery (India, 2021). All these reports clearly show WS significance. Hence, evaluating binding interactions of phytocompounds isolated from WS against SARS-CoV-2 proteins may help to understand its mechanism of action against COVID-19. In this concern, in the present investigation, we selected withanolides with a wide range of therapeutic applications (e.g., antimicrobial, anti-tumour, antiinflammatory, anti-oxidant, anti-stress) for the analysis against the SARS-CoV-2 3CL pro . Initially, all withanolides reported in the WS plant were collected through manual literature search, and their binding efficiency against SARS-CoV-2 3CL pro was evaluated using in silico molecular docking and MD simulation studies. The docking results were additionally validated using the complementarity principle implemented via the AlteQ method.

Collection of phytocompounds
The collection and identification of withanolides from the Withania somnifera plant was manually done through literature curation. The NCBI PubMed database search engine (https://www.ncbi.nlm.nih.gov/pubmed/) was used to collect peer-reviewed research articles. All research articles with the search term ''withania somnifera'' were collected for the analysis. Each article was reviewed separately, and the phytocompounds reported in the WS plant with valid experimental evidence were selected for the analysis. The structure information of the phytocompounds was retrieved from the PubChem database in the SMILES format.

Molecular fingerprinting
A structure based molecular screening of the collected phytocompounds was performed against compounds with a determined SARS-CoV-2 inhibitory activity. A structure based molecular screening was performed using 2D molecular fingerprints (FP) in RDKit (2020.03.1). In the present study, the Morgan circular FP (i.e., the extended connectivity) (Rogers & Hahn, 2010) was used to generate the FP of the phytocompounds (the FP radius was set to four). The molecular similarity exploration was performed using the Tanimoto coefficient (Tc) (Bajusz et al., 2015). The Tc similarity score ranges from zero to one, with zero representing the minimum and one representing the maximum similarity.
The generated FPs were compared to the calculated FP of 8702 molecules collected from the ChEMBL database (Davies et al., 2015) (release version ChEMBL 27) with the assay ID CHEMBL4495582, and the target ID CHEMBL4523582. The assay ID CHEMBL4495582 corresponds to the SARS-CoV-2 inhibitors and target ID CHEMBL4523582 corresponds to the replicase polyprotein 1ab of the SARS-CoV-2.

Molecular docking
Detailed investigation of phytocompounds against SARS-CoV-2 3CL pro was performed by molecular docking studies. Before docking, 3D conformations of phytocompounds were generated and optimized from their SMILES notations using RDKit and MMFF94 force field. The 3D structures were then converted into the pdbqt file format using AutoDockTools 4 script prepare_ligand4.py (Morris et al., 2009). Similarly, the receptor was prepared as described below using UCSF Chimera 1.14 (Pettersen et al., 2004). In the present study, we used the SARS-CoV-2 3CL pro conformation which was generated in our previous study based on the k-means clustering of the unbound 3CL pro during a course of a 900 ns MD simulations (PDB ID: 6LU7) (Novak et al., 2021a). The k-means clustering was based on all non-hydrogen backbone atoms, and it gave two different conformations of the SARS-CoV-2 3CL pro , the primary A conformation and the secondary B conformation, which were present 86.7% and 13.3% of the simulation time, respectively. The structural differences of these conformations compared to the crystallographic conformation of the SARS-CoV-2 3CL pro (6LU7) are thoroughly discussed in (Novak et al., 2021a), while for the present work, the A conformation was chosen, as it is the more prevalent one. The A conformation of the SARS-CoV-2 3CL pro was prepared by adding the Gasteiger charges to all atoms, followed by merging of the nonpolar hydrogens. Such structure was then saved in the pdbqt file format. Docking was performed using AutoDock Vina (Trott & Olson, 2010a) locally on a personal computer with 8 Intel R CoreTM i7-6700K CPU @ 4.00 GHz, 32 GB RAM, and 64-bit Windows 10 Pro operating system. The center of the grid box was set on the Cys145 CA atom (dimensions: x = 13.3, y = 58.2, z = 45.4), with a box size of 20 × 25 × 25 Å. The grid spacing was set to 1 Å and number of modes and exhaustiveness were both set to 100.

Selection of docked poses using the AlteQ method
Validation of the conformations generated by the docking tools can be assessed using different methods (B-H & Brenk, 2009), among which the RMSD (root-mean-square deviation) based methods are the most widely used. A docked ligand pose with an RMSD score less than 2 Å, compared to the reference ligand (the crystallographically determined ligand pose), is generally considered as a good pose (Trott & Olson, 2010b). However, selection of docked poses based solely on RMSD has several flaws and can lead to misclassification of both the correct and incorrect poses (Cole et al., 2005;Kroemer et al., 2004). In the current work, we used the SARS-CoV-2 3CL pro conformation from a molecular dynamics simulation, and hence measuring the correctness of docking poses based on the RMSD with the reference model is not an appropriate method. Besides, the covalently bound peptidomimetic ligand of SARS-CoV-2 3CL pro (6LU7) was used for the MD simulations and measuring the RMSD between covalently bound peptidomimetic ligand with the non-covalently docked ligand is also not appropriate. Therefore, in the present work distance-based measures (i.e., ligand-receptor contacts) were used to check for the correctness of the docked poses. Many distance-based measures have been developed (Rueda et al., 2010;Hawkins et al., 2008), and in these methods, the cut-off length for the ligand-receptor contacts has to be defined. The complementarity principle coupled with the AlteQ method is a newly developed method for determining ligand-receptor contacts in which no cut-offs have to be introduced (Potemkin & Grishina, 2008); it calculates ligand-receptor contacts based on the electron density overlaps between the ligand and the receptor atoms using the Slater's type atomic contributions (Potemkin, Grishina & Potemkin, 2017;Grishina & Potemkin, 2019). Recently, we used this method to calculate electron density overlaps in EGFR (Kandagalla et al., 2021) and CDK (Rimac, Grishina & Potemkin, 2020) receptor-ligand complexes, where we showed that all interactions determined by overlaps of electron clouds follow the complementarity principle, expressed by the following equation: (1) where ρ ligand represents ligand's contribution to electron density in the m th point of the molecular space, ρ enzyme represents enzyme's contribution to electron density in the same point, and RLRE is defined in a following manner (Eq. 2): where dist ligand (dist enzyme ) represents the distance between the mth point and the ligand's (enzyme's) atom having the highest contribution to ρ ligand (ρ enzyme ) at that point. The above complementarity model (Eq. 1) was used for validation of the docking poses. In this regard, the complementarity model was developed from the experimentally solved crystal structure of the SARS-CoV-2 3CL pro . A total of 29 experimentally solved SARS-CoV-2 3CL pro crystal structures with non-covalently bound inhibitors were selected from the RCSB Protein Data Bank (https://www.rcsb.org/) (Berman et al., 2002) for the analysis (PDB ID's 5R7Y, 5R7Z, 5R80, 5R81, 5R82, 5R83, 5R84, 5RE4, 5RE9, 5REB, 5REH, 5REZ, 5RF1, 5RF2, 5RF3, 5RF6, 5RF7, 5RFE, 5RG1, 5RGH, 5RGI, 5RGK, 5RGU, 5RGV, 5RGW, 5RGX, 5RGY, 5RGZ, 5RH0, 5RH1, 5RH2, 5RH3, 5RH8, 5RHd, 6M2N, 6W79, 7JU7, 7KX5, and 7L5D). On the other hand, the five conformations of each phytocompound showing the lowest binding energy with 3CL pro were collected from docking studies for further analysis. All the collected conformations were prepared using the UCSF Chimera 1.14 (University of California, USA) (Pettersen et al., 2004). The chain A conformations with the ligands were retained for the electron density analysis. The 3D maps of the electron density (ρ) were calculated for all conformations using the in-house developed quantum free-orbital AlteQ method (Potemkin & Grishina, 2008). The AlteQ method represents molecular electron density as a sum of Slater's type atomic increments, and it can be expressed as (Eq . 3): where N is the number of atoms in the molecule, ρA is the A atomic increment in molecular electron density, The AlteQ method separates the molecular density into atomic contributions and allows one to separately consider contributions of the enzyme and the ligand atoms (Eq. 5). The generated electron density maps were used to construct a linear regression model to establish a correlation between the distance and the electron density overlap between the ligand and receptors atoms using Eqs. (1) and (2). The generated 3D electron density maps were statistically processed using the scikit-learn (Pedregosa et al., 2011) and the SciPy library (Virtanen et al., 2020) in Python (v 3.7.6), and the plots were generated using the Matplotlib library (Hunter, 2007). The AlteQ software is freely available online via the ChemoSophia webpage (http://www.chemosophia.com/).

Molecular dynamics (MD) simulations
MD simulations for all ligands (WS1, WS4, WS11, and two conformations of the WS7 ligand) were run in complex with the SARS-CoV-2 3CL pro in the previously acquired conformation (PDB: 6LU7). The best docked ligand positions, which were determined using the AlteQ method, were used as the starting points for the MD simulations. The AMBER ff14SB force field (Maier et al., 2015) was used to model the protein and the GAFF force field (Wang et al., 2004) was used to model the ligand. Other simulation parameters such as periodic boundary conditions, NVT conditions, pressure, temperature were the same as described in our previous article (Novak et al., 2021b;Pathak et al., 2021). In brief, protein-ligand complexes were solvated in a truncated octahedral box of TIP3P water molecules spanning a 12 Åthick buffer, and Na + and Cl − ions were added according to (Machado & Pantano, 2020) to achieve a neutral environment with a salt concentration of 0.15 M. Such structures were then submitted for geometry optimization in the AMBER16 program (Case et al., 2016) employing periodic boundary conditions in all directions. For the first 1,500 cycles, the complex was restrained and only water molecules were optimized, after which another 2,500 cycles of optimization followed where both water molecules and the complex were unrestrained. Optimized systems were gradually heated from 0 to 310 K and equilibrated during 30 ps using NVT conditions, followed by productive and unconstrained MD simulations of 300 ns employing a time step of 2 fs at constant pressure (1 atm) and temperature (310 K), the latter held constant using Langevin thermostat with a collision frequency of 1 ps −1 . Bonds involving hydrogen atoms were constrained using the SHAKE algorithm (Ryckaert, Ciccotti & Berendsen, 1977), while the long-range electrostatic interactions were calculated employing the Particle Mesh Ewald method (Darden, York & Pedersen, 1993). The non-bonded interactions were truncated at 11.0 Å. Analysis of the trajectories was performed using the cpptraj module of AmberTools16 (Roe & Cheatham, 2013).

Binding free energy calculations and decomposition
The binding energies, G BIND , of the simulated complexes were calculated using the MM-GBSA (Molecular Mechanics-Generalized Born Surface Area) and the MM-PBSA (Molecular Mechanics-Poisson-Boltzmann surface area) protocols (Genheden & Ryde, 2015;Hou et al., 2011), available as a part of AmberTools16 (Ferenczy, 2015). G BIND is calculated from snapshots of MD trajectories (Ferenczy, 2015) with an estimated standard error of 1-3 kcal/mol (Genheden & Ryde, 2015). G BIND is calculated in the following manner: where the symbol <> represents the average value over 100 snapshots collected from the last 30 ns part of the corresponding MD trajectories (every 150th frame was taken for the calculation). The calculated MM-GBSA and MM-PBSA binding free energies were decomposed into specific residue contribution on a per-residue basis according to established procedures. This protocol calculates the contributions to G BIND arising from each amino acid side chains and identifies the nature of the energy change in terms of interaction and solvation energies (Gohlke, Kiel & Case, 2003;Rastelli et al., 2010). The entropy term was not calculated.

Collection of Withania somnifera compounds from literature and molecular fingerprinting
A search of the NCBI PubMed database resulted in 1,401 research articles which included the term ''withania somnifera'', with the dates until February 28, 2021. All the articles were reviewed and withanolides were collected. Altogether, a total of 80 compounds were collected. However, several studies have already reported WS withanolides against the SARS-CoV-2 3CL pro andhence these compounds were excluded from the analysis. A total of 11 specific types of withanolides, such as withaphysalins, withasomniferols, and withafastuosins, were considered for a detailed binding interaction since there were no previous studies focused on their binding interactions with the SARS-CoV-2 3CL pro . Details about withanolides selected in the present investigation and their respective article information are provided in File S1. Among the selected withanolides, withasomniferol A, withasomniferol B, and withasomniferol C are newly isolated WS phytoconstituents (Anjaneyulu & Rao, 1997).
Similarity of the collected phytocompounds with reported SARS-CoV-2 inhibitors was checked using a molecular similarity analysis. Tc scores of the phytocompounds were calculated relative to the reported SARS-CoV-2 inhibitors, and they ranged from zero to 0.4 (File S2).

Molecular docking
From all tested phytocompounds, compounds with ID's WS1, WS4, WS7, and WS11 showed the lowest energies of −8.0, −8.2, −7.6, and −7.8 kcal/mol, respectively, when binding to the SARS-CoV-2 3CL pro (File S3). Their molecular interactions with the SARS-CoV-2 3CL pro active pocket residues are shown in Figs. 1 and 2, and 2D structures of phytocompounds are shown in Fig. 3. Their binding energies against the 3CL pro and their ADME properties are shown in Table 1. The analysis of the ADME properties of potential drug candidates is essential in the early stage of drug discovery to reduce failure rates in the clinical phase of drug discovery. According to Lipinski's rule of 5 (RO5), to be orally active, drug-like compounds should have molecular weight below 500, number of hydrogen bond donors below 5, number of hydrogen bond acceptors below 10, and the log P value should not exceed 5. As shown in Table 1, all the selected compounds obey these rules except for the WS1 compound, which has a molecular weight above 500. All four compounds were found to be non-inhibitors of the cytochrome P450 enzymes (CYP3A4, CYP2D6, CYP2C9, CYP2C19, and CYP1A2). The cytochrome P450 enzymes play an essential role in metabolism of various molecules. Their inhibition can cause serious drug-drug interactions that can cause unanticipated adverse effects (Lynch & Neff, 2007). Since all our lead compounds were predicted to be non-inhibitors of the cytochrome P450 enzymes, this suggests that they possess a low potential for drug-drug interactions. Additionally, none of the suggested lead compounds are Pan-assay interference compounds (PAINS), which indicates that they will not give false positive results in high-throughput screens. All PAINS interact nonspecifically with numerous biological targets, as opposed to with specific targets. Furthermore, stability of these compounds with SARS-CoV-2 3CL pro was verified through MD simulations. Additionally, before the MD analysis, poses of the compounds generated by docking were evaluated through the complementarity centered AlteQ method.

Selection of the docked poses with the complementarity-centered AlteQ method
Firstly, a linear regression model was constructed to establish a correlation between the distance between the ligand and receptors atoms and their electron density overlap (Eq. 1) for the experimental SARS-CoV-2 3CL pro complexes and this information was used for evaluating the correctness of the docked conformations. The correlation coefficients (R 2 ) in the 5R7Y, 6M2N, 5RGH, 5RF7, 5RG1, 5RH8, and 5REB complexes were found to be below 0.50 and hence were not included in model generation. The average R 2 of all the other experimental conformations was 0.732 (File S4). In a recent article, we described two ways of using the complementary based AlteQ method for validating docking scores . In the present work, we evaluated the correctness of the docked poses from the intercept-versus-slope graph. The intercept and the slope coefficients for all experimental complexes collected from Eq. 1 are correlated and fall onto the same regression line (R 2 = 0.97). The following model (Eq. 6) was obtained (Fig. 4).

Figure 4 A linear regression model established by correlating the intercept and the slope (
Full-size DOI: 10.7717/peerj.13374/ fig-4 The distance between the receptor and the ligand decreases with the slope, as can be observed in Fig. 4. In the 5RH1 complex, the ligand and the receptor atoms are very close, and the slope value was found to be −1.83. In contrast, the slope value of the 5R81 complex was found to be very high (−0.88), which indicates long distances between the ligand and the receptor atoms. The regression line obtained from the experimental SARS-CoV-2 3CL pro complexes was used a reference line to measure appropriateness of the generated docked conformations.
The top five docking poses of WS1, WS4, WS7, and WS11 were evaluated, as they show the lowest binding energies when binding to the 3CL pro . Their slope and intercept values were collected from Eq. 1 and plotted against the experimental conformations (Fig. 4). As expected, different electron density overlap patterns were observed in the docked conformations; they are slightly away from the regression line of the experimental conformations due to different conformations and their orientation within the active pocket, and they correspond to different interaction patterns. The distance of the docked conformations from the regression line was measured and the complete information is provided in File S5. In the top five conformations, the distance of the conformation 1 (based on the docking results) of WS1(0.062), WS4 (0.045), WS11(0.105) ligands were found to be much lower compared to the other four conformations, which is in accordance with the experimental conformations. The only exception was conformation 1 of the WS7 ligand, whose distance was found to be much greater (0.111). The electron density overlap patterns of a few conformations of the ligands WS4, WS7, and WS11 were more distant from the rest of the points (Fig. 4, dotted circle), and these conformations show a unique pattern in the electron density overlap, which was not observed in the experimental conformations.

Molecular dynamics simulations and free energy calculations
To check for the stability of the conformations, MD simulations were performed. For ligands WS1, WS4, and WS11 the best docked conformations were chosen since their electron density overlap patterns were very similar to the electron density overlaps of the experimental conformations. For the WS7 ligand, two different conformations were considered, namely the first (−7.6 kcal/mol, WS7_v1) and the fifth conformation (−7.3 kcal/mol, WS7_v2), with both conformations having different electron density overlap patterns with different binding energies. The distance of the conformation 1 of the WS7 ligand (0.111) was found to be much higher when compared to the fifth conformation (0.062) from the regression line, but the fifth conformation showed a unique pattern in the electron density overlap, and its points (i.e., slope and intercept values from Eq. 3) are located far away from the rest of the points (Fig. 4, dotted circle).
From the five complexes tested, three complexes were found to be stable for the entire duration of the simulation (WS4, WS7_v1, and WS11), while two complexes showed a brief dissociation (∼8.5 ns in the case of WS7_v2 and ≈ 26.0 ns in the case of WS1), which can also be seen from Fig. 5. The root-mean-square fluctuation (RMSF), radius of gyration and intermolecular hydrogen bond plots are provided in the File S6.
From Table 2, it can also be seen that WS1 and WS7_v2 bind the weakest (calculated using the MM-GBSA approach). From the three stable compounds (WS4, WS7_v1, and WS11), WS7_v1 and WS11 have approximately the same G BIND , which is significantly lower than that of WS4. Additional calculations of G BIND were performed using the MM-PBSA protocol. These results are completely in accordance with the MM-GBSA results and show the same trend; WS7_v1 and WS11 have virtually the same binding constant (∼−26.00 kcal/mol), with the other three ligands having at least 12 kcal/mol higher G BIND (detailed information can be found in File S7). The top ten contributing amino acids for the binding of the two best ligands (i.e., WS7_v1 and WS11) are shown in Table 3 and their positions relative to the catalytic Cys145 residue are shown in Fig. 6. Relative to the figure, the WS7_v1 ligand is located slightly more to the left, and the WS11 ligand slightly more to the right. This results in a fact that these ligands share nine out of ten most contributing residues (Table 3), with exception being Thr45, which interacts only with the WS7_v1 far left hydroxy group, and Met165, which interacts only with the WS11 cyclopentanyl group (Fig. 7, Table 3, bolded). While the hydrogen bond interaction with Thr45 is moderately important in case of WS7_v1 binding, the van der Waals interaction with Met165 seems to be the most important one for WS11_v1 binding. Additionally, the WS7_v1 ligand forms two very strong interactions with partial G BIND lower than −2.0 kcal/mol, with one of them being the Cys145 residue, a crucial residue for the 3CL pro function (Cole et al., 2005), and the other one being Thr25 which forms a very strong hydrogen bond with the keto oxygen atom of WS7. On the other hand, WS11 forms much stronger interactions with the other crucial residue, (His41), as well as with Cys44, with both interactions being hydrophobic in nature, while in the case of WS7_v1 these interactions have a van der Waals character. This indicates that WS7 and WS11 interactions have a preference for different parts of the binding pocket and that they could both be good candidates for 3CL pro inhibitors. This also implies that there is still room for developing or finding a ligand which could interact more strongly with both parts of the binding pocket

DISCUSSION
Natural products play a significant role in the discovery of novel and effective therapeutics to fight the present COVID-19 pandemic. Herbal extracts and spices are natural immune boosters and/or anti-infective agents currently used in many parts of the world (Gbadamosi, 2020). In traditional folk medicine, spices, botanical detoxifiers, antioxidants (Gbadamosi & Afolayan, 2016) and plant hematinics (Gbadamosi, 2012) are used as antiviral mediators to prevent or minimize the impact of various diseases. A lack of targeted treatments has encouraged the exploration of novel drug lead compounds in which computational approaches offer a comparatively fast and cost-effective approach. COVID-19 is characterized by a disrupted immunological balance, hyperinflammation, cytokine storm, and multiorgan failure (Saggam et al., 2021). WS plant is reported to mitigate early disease progression and protect vital organs (Saggam et al., 2021). In addition, Saggam et al. (2021) suggested that the WS plant's antiviral capabilities could disrupt viral entrance and its life cycle. The WS phytocompounds were explored for their antiviral potential against SARS-CoV-2 proteins using a computational molecular docking tools. Also, recent reports established that the phytocompound withanone disrupts the host-virus interaction by destabilizing the ACE2-spike protein receptor-binding domain complex (Balkrishna et al., 2021). Withaferin A and withanone were proven to prevent viral entry and replication by blocking the 3CL pro and TMPRSS2 enzymes (Kumar et al., 2022). Kushwaha et al. (2021) tested non-characteristic phytocompounds (quercetin-3-rutinoside-7-glucoside, rutin, and isochlorogenic acid B) present in the WS plant against SARS-Cov-2 3CL pro through molecular docking studies. All the above reports focused and discussed computational approaches and further concluded the effect of WS phytocompounds in blocking the viral entry and replication. However, the reports of WS compounds against SARS-CoV-2 are lacking and the amount of experimental and clinical data is limited. This highlights the necessity of performing systematic research in COVID-19 to investigate WS phytocompounds as antiviral therapeutics. According to Lurie, Keusch & Dzau (2021), given recent indications that the SARS-CoV-2 pandemic will be a long-term health problem, there will be a considerable demand for SARS-CoV-2 medicines and adjuvants research and development. Therefore, in the present study, a total of 11 reviewed withanolides from WS, belonging to the group of withaphysalins, withasomniferols, and withafastuosins, were collected through manual literature curation. The initial molecular docking analysis revealed that the binding energy of all the selected phytocompounds against the SARS-CoV-2 3CL pro