MD simulation and MM/PBSA identifies phytochemicals as bifunctional inhibitors of SARS-CoV-2

Abstract The global spread of SARS-CoV-2 has resulted in millions of fatalities worldwide, making it crucial to identify potent antiviral therapeutics to combat this virus. We employed structure-assisted virtual screening to identify phytochemicals that can target the two proteases which are essential for SARS-CoV-2 replication and transcription, the main protease and papain-like protease. Using virtual screening and molecular dynamics, we discovered new phytochemicals with inhibitory activity against the two proteases. Isoginkgetin, kaempferol-3-robinobioside, methyl amentoflavone, bianthraquinone, podocarpusflavone A, and albanin F were shown to have the best affinity and inhibitory potential among the compounds, and can be explored clinically for use as inhibitors of novel coronavirus SARS-CoV-2. Communicated by Ramaswamy H. Sarma


Introduction
A novel strain of human virus was reported during the latter part of December 2019 from the province of Wuhan in China. Its infection in humans was traced to zoonotic transmission from bat (York, 2020;Zhou et al., 2020). Due to the rapid spread of SARS-CoV-2 across multiple continents over a short span of a few months, the WHO declared SARS-CoV-2 as an International Public Health Emergency (PHEIC) on January 30, 2020 (Statement on the Second Meeting of the International Health Regulations (2005) Emergency Committee Regarding the Outbreak of Novel Coronavirus (2019-NCoV), 2020). From boosting PPE production and emphasizing the significance of social distancing, using mask, the Emergency Use Authorization (EUA) of remdesivir/ therapeutic antibodies and the use of the well-known glucocorticoid dexamethasone, the globe has taken important steps to combat this disease (Li et al., 2021).
In the absence of any specific drug to combat this virus, clinical treatment has been limited to the synergistic use of antibiotics and existing antiviral agents to treat associated symptoms (Li & De Clercq, 2020;Morse et al., 2020). Patients have been prescribed chloroquine, remdesivir, hydroxychloroquine, and combinations of other antimicrobial drugs based on their symptoms (Colson et al., 2020;Wang et al., 2020). Furthermore, research on SARS-CoV-2 genome and proteome is expanding our understanding of the virus almost every other day. Since the commencement of the COVID-19 outbreak, many COVID-19 vaccine candidates have entered clinical trials in less than 6 months and have been conditionally approved in less than 10 months, setting a new record for vaccine development speed. In the WHO EUL/PQ evaluation procedure as of July 2021, ten vaccines had been finalized.
Coronaviruses (CoV) are the largest RNA viruses in the coronaviridae family, with genome sizes ranging from 27 to 32 kb (Marra, 2003). It is a single-stranded positive enveloped virus that infects humans as well as animals (Kim et al., 2020, p. 1). The viral RNA is translated into two large overlapping polyproteins (pp1a and pp1b) (Chan et al., 2015;Irigoyen et al., 2016). This polyprotein complex is cleaved by a series of proteolytic enzymes, papain-like protease (PLpro) and mainprotease (Mpro), which releases vital proteins that assist in the viral maturation process. Mpro aids in the polyprotein processing while PLpro proteolytically cleaves the polyprotein from its amino-terminal (Mengist et al., 2020). Since Mpro and PLpro are involved in the proteolytic processing of the virion polyprotein, a process vital to the maturation and production of infectious virions, therefore these proteases can serve as potential therapeutic targets. We hypothesize that inhibiting these proteases can halt process of maturation and production of infectious virions and thus stop the spread of the disease. Several crystal structures of Mpro and PLpro complexed with covalent and noncovalent inhibitors have been reported, which has aided in the identification of potential antiviral compounds (Kirchdoerfer & Ward, 2019;Lan et al., 2020;Mengist et al., 2020;Osipiuk et al., 2020; RCSB PDB -6W63: Structure of COVID-19 Main Protease Bound to Potent Broad-Spectrum Non-Covalent Inhibitor X77, n.d.; Tian et al., 2020;Zhang et al., 2020). So far, many natural molecules have been reported with antiviral potential against Mpro (Aanouz et al., 2021; and PLpro (Kumar, Kashyap, et al., 2021;Kumar, Sharma, et al., 2021). This study focused on screening phytochemicals with potential to inhibit SARS-CoV-2 Mpro and PLpro. Structure-assisted virtual screening was performed using the crystal structure of Mpro (PDB ID: 6W63) (RCSB PDB -6W63: Structure of COVID-19 Main Protease Bound to Potent Broad-Spectrum Non-Covalent Inhibitor X77, n.d.) and PLpro (PDB ID: 6W9C) (Osipiuk et al., 2020) to identify phytochemicals for novel Coronavirus inhibition. Candidate inhibitors were restricted to natural compounds that were retrieved from sweetlead and foodb databases. Further, Molecular Dynamics (MD) simulations were conducted on the obtained complexes to rationalize the affinity and stability of the protein-ligand complexes.

Target preparation
The structures for Mpro (PDB ID: 6W63) and PLpro (PDB ID: 6W9C) used in this study were retrieved from RCSB PDB. The protonation state of Mpro and PLpro was examined using Hþþ server, and all missing hydrogen atoms were added (Gordon et al., 2005). The structures were processed using AutoDock4 (Morris et al., 2009) and converted into pdbqt for use in virtual screening with PyRx (Dallakyan & Olson, 2015).

Ligand preparation
The 3 D coordinates of molecules were retrieved from foodb (www.foodb.ca) and sweetlead (Novick et al., 2013) database. These molecules were processed by adding the Gasteiger charge in Open Babel. The molecules were minimized with mmff94 forcefield and converted into pdbqt format in the PyRx module.

Virtual screening
For this study, we used the AutoDockVina module of PyRx for rigid docking in which the ligands were kept flexible. An exhaust parameter of 8 was selected, and a total of nine poses were generated for each docking. A docking grid with dimension X: 22.764 Y: 25.00 Z: 25.00 and center X: À11.61 Y: 14.19 Z: 70.29 were used for Mpro. Likewise, a grid of dimension X: 22.76 Y: 25.00 Z: 25.00 with center X: À26.6363 Y: 45.4429 Z: À11.09 was used for the PLpro. The grid was selected in such a manner that it encircles the protein's active site residues. Each ligand was docked in the particular protein's binding site using the AutoDock Vina module of PyRx and assessed based on binding affinity. PyMOL was used to visualize the top-ranked molecules, and the best conformations with high affinity were selected.

Docking of the co-crystallized ligand
The ligands X77 and P85 are the inhibitors reported against Mpro and PLpro, respectively, in their co-crystallized forms. These ligands were included in the same protocol as described above for the virtual screening to reproduce the same conformation as present in its co-crystallized form to validate our virtual screening protocol.

Admet prediction
We sought to predict the ADMET properties (Absorption, Distribution, Metabolism, Excretion, and Toxicity) of the compounds that were screened from the virtual screening process. The pkCSM server (Pires et al., 2015) was used to predict various parameters such as water solubility (log mol/ L), intestinal absorption in human (% Absorbed), total clearance (log ml/min/kg), AMES toxicity, maximum tolerated dose in human (log mg/kg/day), BBB (blood-brain barrier) permeability, and CYP450 (Cytochrome P450) substrate.

MD simulation using GROMACS
The complex obtained after virtual screening were further subjected to MD-simulation using GROningenMAchine for Chemical Simulations (GROMACS) (Abraham et al., 2015). GROMOS43a1ff in GROMACS was used for the generation of protein topology files and for assigning the partial atomic charges. Ligand topologies were generated using the PRODRG tool (van Aalten et al., 1996). The protein was solvated in a cubic box of radius 1.5 nm using a simple point charge (PC) water module (Mark & Nilsson, 2001). The protein-ligand in a box was further neutralized by the addition of Na þ /Clion using the genion tool in gromacs. This complex was minimized in 50000 steps by the steepest descent algorithm using the energy cutoff of 10.0 kJmol À1 . Followed by this, equilibration was done in two steps: NVT (constant number of particles, Volume, and Temperature) equilibration at 300 K using the Parrinello-Rahman barostat pressure coupling method for 1 ns (Parrinello & Rahman, 1981); NPT (Number of particles, Pressure, and Temperature) equilibration at a pressure of 1.0 bar for 1 ns using Berendsen thermostat (Darden et al., 1993). The Linear Constraint Solver (LINCS) algorithm was used to limit the lengths of heavy atom (Hess, 2008). The short-range forces were measured using the verlet cutoff scheme with a minimum cutoff of 1.2 nm (Grubm€ uller et al., 1991). Finally, three independent simulations of each 100 ns with different initial velocities were carried out using the leap-frog algorithm (Van Gunsteren & Berendsen, 1988) and time step of 2 fs was used. The analyses of simulations were done using snapshots collected every 10 ps and analyzed for the stability of the polypeptide structure by calculating Root Mean Square Deviation (RMSD), Root Mean Square Fluctuation (RMSF), hydrogen-bond formation, Solvent Accessible Surface Area (SASA), Radius of gyration (Rg), and Principal Component Analysis (PCA) using GROMACS inbuilt tools.

MM/PBSA calculations
Calculation of the binding energy of the complex was performed using the g_mmpbsa package (Kumari et al., 2014). Electrostatic energy, Van der Waals energy, and free solvation energy (including polar and non-polar solvation energies) were calculated to estimate the binding energy of the complex, neglecting the entropy calculations.

Principal component analysis
The principal component analysis (PCA) approach was used to separate global and local motions in a macromolecule (Amadei et al., 1993) and to investigate the relative conformational dynamics and atomic variations of substructures in their native and ligand-bound forms. PCA analysis was performed on the trajectories produced during the 100 ns  molecular dynamics simulation. By extracting the translational and rotational motions and tracking the relative Ca backbone atomic variations, a cross-correlation matrix was developed using the GROMACS gmx covar method. The eigenvalues depicted the energetic contribution of the corresponding principal component (PC), whereas the vector represented the direction of motion. The overall versatility of protein was determined by analyzing the predicted eigenvalues and eigenvectors. The GROMACS gmx anaeig tool was used to compute the 2 D projections with respect to the first two eigenvectors.

Compounds identified from virtual screening
The virtual screening protocol followed in this study was validated by docking of bonafide inhibitors X77 and P85 against Mpro (6W63) and PLpro (6W9C), respectively. The conformation with the lowest docking energy score was selected and compared with the pose of the standard inhibitors in the crystal structure ( Figure 1C). The RMSD values obtained for the docked conformations of X77 and P85 over their respective binding pose in the crystal structure of Mpro and PLpro was 0.706 Å and 1.02 Å, respectively. Therefore, the docking was considered reasonable for validating the reliability of the PyRx parameters for screening and docking of compounds with proteins. The binding affinity of the respective ligands was used as a criterion for selection of potential compounds during virtual screening and a total of 14 compounds were shortlisted.
The 14 hits that were obtained using Mpro as target were then screened for their ability to inhibit PLpro. PLpro in complex with the ligands exhibited a unique mechanism of inhibition, in which the inhibitor occupied a site between the two subsites, one formed by residue Glu161-Leu162 and the other formed by residues Gly163, Pro247, Glu263, Asn267, His272, and Ile300. This induced conformational     interactions and are illustrated in Figure 1A and Supplementary Figure 1B.

ADMET prediction
The values obtained from pkCSM server, given in Table 1 indicated that all the compounds were fairly soluble in water. All the compounds but kaempferol-3-robinobioside and quercetin 3-O-xylosyl glucuronide can be absorbed in humans through the intestinal route. None of the compounds were predicted to act as CYP 450 substrate. All the compounds except isoquercitrin were predicted to be negative in AMES toxicity test.

Analysis of molecular dynamics simulation
The shortlisted hits were then subjected to MD simulations using the Gromos54a7 forcefield of GROMACS package in triplicate of 100 ns with a different initial random speed. The MD trajectories were analyzed using various GROMACS in-built tools.

RMSD analysis
RMSD analysis was performed to assess overall structural divergence from the crystal structure in response to ligand binding. The root mean square deviation of all the structures were in the range of 0.26 nm to 0.37 nm with reference to the initial structure. The apo-Mpro and apo-PLpro retained stability until 100 ns. It was observed that the apo-Mpro had higher deviation compared to the docked complexes thus implying that the docked ligands stabilized the Mpro and PLpro structures. All the structures exhibited typical rise in RMSD up to 25 ns due to the kinetic force shock induced by the initiation of molecular dynamics. After 25 ns, all the complexes attained stable plateau except for the complexes of Mpro with bianthraquinone and isoquercitrin. Figure 2 (and supplementary Figure 2) illustrates that the RMSD values  reached an equilibrium and remained stable thereafter and attained equilibrium. The average RMSD values obtained for apo-Mpro, apo-PLpro, and all the complexes are summarized in Table 2 (and supplementary Table 2). In the case of PLpro, the structure of complex of PLpro with withastramonolide, 3, 8 0 -biapigenin, theasinensin B, and isoginkgetin showed higher average RMSD values as compared to apo-PLpro. All the other compounds recorded a thoroughly stable RMSD profile in complex with PLpro. The RMSD profile of the complexes of most of these phytochemicals with Mpro and PLpro were found similar to the apo-Mpro and apo-PLpro. As evident from Figure  2, the binding of each phytochemical to Mpro and PLpro was stable in its respective ligand-binding pocket.

RMSF analysis
Further, the fluctuation of the backbone coordinates from their average position for each protein-ligand complex as well as their apo structures were also calculated. The structural divergence values were plotted against residues of the protein to reflect the flexibility and the motion of the proteins throughout the simulation. The RMSF plot of the apo-Mpro, apo-PLpro and their complexes with ligands are shown in Figure 3 (and supplementary Figure  3). As anticipated, the terminal residues had higher RMSF values, which suggested their mobile nature. The secondary structure elements in the protein-ligand complexes recorded minimal deviation as compared to the apo-protein. The RMSF of loops in the ligand-binding site was significantly higher in ligand complexes than in apo structures. A minor increase in the average RMSF value (Table 2 and supplementary Table 2) was observed for methyl amentoflavone (0.17 nm), albanin F (0.16 nm), quercetin 3-O-xylosyl glucuronide (0.17 nm), and withastramonolide (0.17 nm) when compared to the apo-PLpro (0.15 nm).

Radius of gyration (Rg) analysis
Rg analysis was carried out to monitor the changes in protein compactness in response to ligand binding. The average values of Rg for apo-Mpro and apo-PLpro were observed to be 2.12 nm and 2.24 nm, respectively. The average Rg values for different Mpro and PLpro complexes indicated only subtle changes in compactness of proteins in response to ligand binding (Table 2 and supplementary Table 2). However, the average Rg values for the complexes of PLpro with 3, 8 0 -biapigenin and neotheaflavin 3-gallic acid with 2.25 nm, and podocarpusflavone A and theasinensin B with 2.6 nm, were slightly higher as compared to the apo-PLpro. As shown in Figure 4 (and supplementary Figure 4), all complex systems exhibited relatively similar and consistent Rg values during the course of the simulation thus suggesting the compact packing of the secondary structure elements.

SASA analysis
SASA analysis was performed to understand the behavior of the protein (hydrophilic or hydrophobic) towards the solvent and was calculated by considering all the polar and nonpolar interactions between the atoms. The average SASA values obtained for apo-Mpro and PLpro were found to be 135.8 nm 2 and 146.3 nm 2 , respectively. Moreover, the average SASA values obtained for various Mpro and PLpro-complexes exhibited no significant changes when compared to their apo structures (Table 2 and supplementary Table 2). The results of SASA analysis for apo-Mpro, PLpro, and their complexes ( Figure 5 and supplementary Figure 5) indicated their stability during the simulation.

Hydrogen bond analysis and center-ofmass distance
The H-bond utility of GROMACS was utilized to determine the number and apportioning of hydrogen bonds between the protein and ligand during the simulation. The average number of hydrogen bonds was recorded in the range of 1-8 for the Mpro-ligand complexes and 1-4 for the PLpro-ligand complexes. As illustrated in Figure 6 (and supplementary  Figure 6), the optimum number of hydrogen bonds persisted throughout the simulations, thus implying the stability of the complexes. Over the 100 ns MD simulation, the center-ofmass (CoM) distance between the compounds and the residue HIS41 in Mpro and GLN167 in PLpro was measured and represented in Figure 10 (and supplementary Figure 10).

MM/PBSA binding free energy calculation
The binding energy of the complexes obtained from the g_mmpbsa package are shown in Tables 3A and 3B (and  supplementary Table 3 A and 3B). The estimated MM/PBSA binding free energy with respect to time is shown in Figure  9 (and Supplementary Figure 9) which suggest the efficient and steady binding of the ligands to the proteases under consideration. Isoginkgetin, kaempferol-3-robinobioside, methyl amentoflavone, bianthraquinone, podocarpusflavone A, and albanin F recorded higher binding energy than the co-crystallized inhibitor X77 with Mpro. All the ligands also showed considerably higher binding energy with PLpro than the reported inhibitor P85.

Principal component analysis
The Essential dynamics (ED) or Principal Component Analysis (PCA) was used to monitor protein expansion during simulations. The large-scale average motility of the apo-Mpro, apo-PLpro, and their complexes with phytochemicals was identified by PCA and indicated structures that underlie the atomic fluctuation. The overall movement of the system could be cumulatively expressed by the summation of eigen values (Figure 7 and Supplementary Figure 7). Multiple eigen values are an indication of a broad spread of the protein's conformational space. The module gmx anaeig reads the set of vectors and eigenvalues as input files to project an MD trajectory with selected vector values. Although the total motion of the atoms was distributed over many eigenvectors, the top 10 eigenvectors, sorted by their corresponding eigenvalues, accounted for more than 70% of the total collective motion. The traces of the covariance matrix obtained after diagonalization is given in Table 4 (and Supplementary  Table 4). The complexes of Mpro with all the phytochemicals cover a lesser traceable area in comparison to the apo-Mpro. On the other hand, the complexes of PLpro with bianthraquinone, isoginkgetin, neotheaflavin 3-gallic acid, quercetin 3-Oxylosyl glucuronide, and theasinensin B occupied a greater traceable area than apo-PLpro. These observations suggest that the above-mentioned protein-complex structures had higher flexibility in comparison to the apo structure. The 2 D trajectory projections of the vectors in Figure 8 (and Supplementary Figure 8) show a superposition of apo-  protein with its complexes. Briefly, this analysis provided us with wide-scale motility in the protein in response to ligand binding.

Discussion
Nature offers a diverse set of molecules for study and development of drugs for a variety of ailments, including viral infections. Dietary supplements were designated as a food category in 1994 by the Dietary Supplement Health and Education Act, which means they do not require FDA premarket approval (Dietary Supplement Health and Education Act of 1994, n.d.). Natural compounds, in addition to the currently available drugs, could be an interesting alternative in managing the COVID-19 pandemic. Due to insufficient oxygen/nutrient availability and abnormal metabolism of endogenous/exogenous components to generate ATP in their cells, individuals with contagious and non-infectious illnesses have a weaker immune system (Sattler, 2017). The risk of mortality can be higher in COVID-19 patients due to the compromised natural immunity which can be mitigated by routine consumption of immune boosters prescribed in traditional medicine. Phytochemicals have been shown to synergize with traditional drugs and, interact with the cellular targets and intercept several pathways without any systemic toxicity (Paraiso et al., 2020). Based on this, we attempted to probe phytochemicals as immune boosters and assessed their antiviral-activity against SARS CoV-2. We screened the phytochemicals for their ability to inhibit the proteases that are essential for the viral replication machinery. Several studies have reported on the application of virtual screening in anti-coronavirus drug discovery and development (Choudhary et al., 2020;Dhankhar et al., 2020;Ibrahim, Abdeljawaad, et al., 2020). Prior to actually carrying out animal and clinical studies, the cell-based assays can be used to assess the efficacy and toxicity of compounds that have been shortlisted using virtual screening. In our study, we utilised virtual screening to screen 35405 compounds to obtain 14 hits that could simultaneously inhibit the two proteases in SARS-CoV-2. These were kaempferol-3-robinobioside, isoquercitrin, isoginkgetin, podocarpusflavone A, 3, 8 0 -biapigenin, aurasperone D, quercetin 3-O-xylosyl glucuronide, withastramonolide, epigallocatechin 3,5-digallate, vitamin D2, theasinensin B, albanin F, 3,5 digalloylepicatechin, and neotheaflavin 3-gallic acid. Each of these compounds exhibited better binding energy than the already reported inhibitors of Mpro and PLpro (X77 and P85 respectively).
Herbs and phytochemicals have a long history of use in the treatment of a variety of respiratory ailments. Flavonoids, a type of phytochemical, have been shown to play a role in regulating cellular function and eliminate free radicals, which cause oxidative stress in the body. Quercetine 3-o-xylosylglucuronide, isoquercitrin, albanin F, and kaempferol-3-robinobioside, which were shortlisted as potential hits in this study, belong to this class of phytochemicals. Quercetine 3-O-xylosyl-glucuronide (found in apache fruit, blackberry, and raspberry), and kaempferol-3-robinobioside (found in Piper nigrum (black pepper berries) and Annona coriacea (sugar apple fruit) have been reported to exhibit antioxidant activity (Cho et al., 2004;Gu et al., 2018;Rommel & Wrolstad, 1993;Zhu et al., 2018). A variety of foods, including Mangifera indica (mango) and Rheum rhabarbarum (stalk of rhubarb), Annona reticulata (leaves of custard apple), and Camellia sinensis (leaves of tea), contain isoquercitrine which has been shown to be effective against the Zika virus by preventing the internalization of virus particles into the host cell     (Cojocaru et al., 2020;Gaudry et al., 2018;Panda & Kar, 2007;Sakakibara et al., 2003). Albanin F, a flavonoid present in Morus alba (root barks of white mulberries), has been reported to exhibit anti-inflammatory effects (Lim et al., 2013;Papa et al., 2019). Isoginkgetin, a biflavonoid extracted from flowers of Ginkgo biloba has been reported to exert antiviral effects (Haruyama & Nagata, 2013;Isah, 2015). Amentoflavone, a class of bioflavonoids has been demonstrated for several bioactivities such as antioxidant (Saroni Arwa et al., 2015), antiinflammatory (Abdallah et al., 2015), antiviral (Coulerie et al., 2013), and antifungal (Hwang et al., 2012). Most of these compounds are obtained from plants of the family Calophyllaceae, Clusiaceae, and Selaginellaceae. Methyl amentoflavone, podocarpusflavone A, and 3,8 0 -biapigenin are the derivatives of amentoflavone. Methyl amentoflavone has been isolated from Selaginella sinensis (Ma et al., 2001), Ginkgo biloba (Lobstein-Guth et al., 1988), and many species of Cupressaceae. Podocarpusflavone A has been reported from Podocarpus macrophyllus, and in several plants that belong to the species Garcinia. 3,8 0 -biapigenin has been commonly isolated from Hypericum perforatum (Bergh€ ofer & H€ olzl, 1989). All these plants have been used since ancient times as phytomedicines.
Another significant category of nutraceuticals is the bianthraquinone, that has been isolated from lichens, fungi, and other higher medicinal plants (e.g. Polygonaceae, Rhamnaceae, Rubiaceae, Fabaceae, and Xanthorrhoeaceae) (Breimer & Baars, 1976;Cirillo & Capasso, 2015;Singh et al., 2020;Yun-Choi et al., 1990). One of the well-studied medicinal plants, Withania somnifera, popularly known as the winter cherry (ashwagandha), has been regarded as an immunity booster plant known to Indian Ayurveda for over 3000 years. The decoction of its root has been used as a food and health restorative (Dar et al., 2015;Kim et al., 2005). The extract from this herb is reported to contain various phenolic and flavonoid compounds that are beneficial for immune system. Withastramonolide, one of the hits in the virtual screening during this study, is one of the natural phytochemicals present in ashwagandha.
The hits from virtual screening also included catechins (theasinensin B, neotheaflavin3-gallate, and 3,5 digalloylepicatechin) that are present in tea, one of the world's most popular beverages prepared from the leaves of Camellia sinensis. These phytochemicals have been proven effective against many viruses, bacteria, fungi, and prions. The major mechanisms through which these exert antiviral effects include inhibition of virus from entry into the host cell (adenovirus, enterovirus, herpes simplex virus, hepatitis B, hepatitis C, human immunodeficiency virus, influenza, and rotavirus); inhibition of synthesis and transcription of viral RNA and DNA (enterovirus, hepatitis B, hepatitis C, and human immunodeficiency virus); and destruction and functional alteration of different virus (adenovirus, herpes simplex virus, and influenza) (Higdon & Frei, 2003;Reygaert, 2018). Another crucial phytochemical hit obtained from the virtual screening was vitamin D2, which is a type of vitamin D found in food that serves as a dietary supplement, also referred to as ergocalciferol or calciferol. Vitamin D is a fatsoluble vitamin that supports calcium absorption, regulates the growth of the bone, and plays a crucial role in immune system (Cardwell et al., 2018).
The stability of the protein-ligand complexes was assessed by MD-simulation of 100 ns in triplicate making a total of 300 ns. The trajectories were extracted and analysed to assess the stability and binding affinity of the structure under consideration. The RMSD of the backbone was evaluated which assessed the stability of the complex during the analysis. The backbone RMSF analysis showed the deviation in the protein as a result of the binding of the ligand, especially in the loop regions that encompass the ligand in the active site, whereas, the secondary structure elements exhibited minimum deviation. The decreased fluctuation in Rg indicated that protein was compactly packed and binding of the inhibitor did not have any adverse effect on the stiffness of the protein. The SASA analysis validated the rigidity and compactness of the complex by assessing the hydrophobic interactions formed between the amino acids. These hydrophobic interactions ensured stabilization of global proteins by shielding non-polar amino acids from the aqueous environment. Also, all of the shortlisted phytochemicals formed a consistent number of hydrogen bonds with the residues in the active site. The binding energy of the complex obtained from the g_mmpbsa package justified the efficient and steady binding of the ligands to both the proteases under consideration.
Isoginkgetin, kaempferol-3-robinobioside, methyl amentoflavone, bianthraquinone, podocarpusflavone A and, albanin F recorded higher binding energy than the co-crystallized ligand X77 with MPro. In the case of PLPro, all the ligands exhibited considerably higher binding energy than the reported inhibitor, P85. Eventually, the wide-scale motility in the protein due to the binding of the ligand was established by PCA analysis. Overall analysis of the MD-run showed that all of the predicted small molecule inhibitors had a high affinity for both the proteases. These phytochemicals have the potential to negatively regulate viral proliferation in host cells, with a stable complex of complex that is stable.

Conclusion
This study adopted a holistic approach while aiming at the two proteases crucial for viral replication and pathogenesis to address the global emergency caused by SARS-CoV-2. Inhibitors of these two proteases hold the potential to arrest the viral replication machinery and hence to combat the infection. The process of virtual screening provided 14 candidate inhibitors for both the proteases in SARS-CoV-2. Among these isoginkgetin, kaempferol-3-robinobioside, methyl amentoflavone, bianthraquinone, podocarpusflavone A, and albanin F exhibited significant higher binding energy in comparison to the already reported inhibitors of Mpro and PLpro (X77 and P85, respectively). The analysis of the results of MD simulations established that all the predicted small molecular inhibitors were located in the protein cavity without hampering the stability of the complex as a whole. In the absence of effective antiviral drugs, these natural phytochemicals may be used as potent therapeutics for a complementary therapy during the treatment of various symptoms of SARS-CoV-2associated infections and other coronaviruses. Due to the ease of their availability in daily consumables, these compounds have an incredibly high potential in providing immunity against coronavirus in humans. While, further research and clinical trials are required to validate the efficacy and toxicity of these compounds, there is the scope to chemically modify them to achieve desired properties for clinical applications.
Author's contribution MS and PK designed the study. MS and PD performed the virtual screening studies and analysis. MS and JKM performed the MD simulation studies. MS, JKM, NN and PK wrote the manuscript. All the authors have approved the final version of the manuscript.

Disclosure statement
No potential conflict of interest was reported by the authors.