key: cord-0999065-xpdpssip
authors: Sawant, Sanjay; Patil, Rajesh; Khawate, Manoj; Zambre, Vishal; Shilimkar, Vaibhav; Jagtap, Suresh
title: Computational assessment of select antiviral phytochemicals as potential SARS-Cov-2 main protease inhibitors: molecular dynamics guided ensemble docking and extended molecular dynamics
date: 2021-07-19
journal: In Silico Pharmacol
DOI: 10.1007/s40203-021-00107-9
sha: 329228135ecbea01f00fae2ea8895c0f9e47ed45
doc_id: 999065
cord_uid: xpdpssip

Covid-19 caused by novel coronavirus, 2019-nCoV or SARS-CoV-2 has become most severe pandemic of this century. No specific therapies are available to treat Covid-19 so far. Recently, main protease (M(pro)), a potential drug target from SARS-CoV-2 has been successfully crystallised. The present study is aimed at assessment of bioactive antiviral phytochemicals as potential SARS-COV-2 M(pro) inhibitors, using ensemble docking, molecular dynamics and MM-PBSA calculations. Ensemble docking studies were performed with Autodock vina program. The top 5 compounds having highest binding free energy were subjected to 100 ns molecular dynamics simulations with Gromacs. The resulting trajectories of converged period of MD were further exploited in MM-PBSA calculations to derive accurate estimates of binding free energies. The MD results were analysed with respect to RMSD, RMSF and hydrogen bond formation and occupancy parameters. The drugs remdesivir and nelfinavir were used as standard drugs for comparative studies. In the docking studies five phytochemicals, dalpanitin, amentoflavone, naringin, hinokiflavone, and rutin were found having lowest binding free energies (< − 10 kcal mol(−1)) which is lower than standard drugs. MD studies suggested that the complexes of these five phytochemicals with M(pro) stabilize with well accepted RMSD. Amongst these phytochemicals, hinokiflavove, amentoflavone and naringin were found having better binding affinity with ΔG(binging) than the standard drug remdesivir. Investigations and validation of these inhibitors against SARS-CoV-2 would be helpful in bring these molecules at the clinical settings. [Image: see text] SUPPLEMENTARY INFORMATION: The online version contains supplementary material available at 10.1007/s40203-021-00107-9.

The entire world is facing unprecedented corona virus pandemic and is distressing public health and economies globally. Corona virus emerged in late December 2019 which causes severe acute respiratory syndrome corona virus-2 (SARS-CoV-2) . As per the WHO report, there are ~ 159,613,203 confirmed cases and ~ 3,317,912 deaths due to the SARS-CoV-2 infection (Website) . Within a short span of merely two years the world has seen a dozen concerning mutations in the virus (Li et al. 2020b) . Earlier, the infection was seen to be severe among the elderly patients and patients with co-morbid conditions such as type 2 diabetes, hypertension, obesity, pre-existing heart and lung diseases, cancer; but currently it is seen infecting people in all age groups (Felsenstein and Hedrich 2020; Richardson et al. 2020) . There is no concrete cure by the drugs, however recently remdesivir (Beigel et al. 2020) , dexamethasone (Schultze et al. 2020) , favipiravir (Agrawal et al. 2020) , nelfinavir (Ohashi et al. 2021 ) were found to reduce risk of mortality in patients with severe COVID-19.

Corona is RNA virus containing RNA-dependent RNA polymerase (Pachetti et al. 2020; Sevajol et al. 2014 ). The full-length genome sequence of SARS-CoV-2 is reported to have 79.6% sequence similarity with SARS CoV-1 ) and both of these seemingly similar viruses uses the cell entry receptor angiotensin converting enzyme II (ACE2) of human host cells. The genomic and proteomic studies of SARS-CoV-2 identified drug design feasible seven major targets which include spike protein, envelop protein, membrane protein, protease, nucleocapsid protein, hemagglutinin esterase, and helicase (Prajapat et al. 2020) . The encouraging results of protease inhibitors lopinavir and ritonavir in the management of COVID-19 (Baden and Rubin 2020) suggests the potential SARS-CoV-2 proteases in drug design. The 26 to 32 kb genome of SARS-Cov-2 contains around 14 open reading frames (ORFs). Amongst these the major reading frame ORF 1ab encodes two overlapping polyproteins (pp1a and pp1ab) (Ullrich and Nitsche 2020) . Two cysteine proteases, chymotrypsin-like cysteine protease (main protease [M pro ] or 3C-like protease [3CL pro ]) and papain-like protease [PL pro ] cleaves these polyproteins into 16 non-structural proteins (NSP1-16) (Prajapat et al. 2020) . The M pro is highly conserved in its sequence throughout coronaviruses (Mirza and Froeyen 2020) . The role of M pro is established in proteolytic cleavage at 11 sites at C-terminus of viral polyprotein 1ab (Stobart et al. 2012) . Hence, this unique feature makes M pro as high profile and best characterized target for anti-SARS therapy (Anand 2002) . Monomeric form of M pro is inactive, while homodimeric is active in catalytic regulation. The structure of M pro comprises of three domains. Domain I and II (β-barrel) contain amino acid residues 8-101 and 102-184, respectively, while domain III (α-helix) constituting amino acid residues 185-200 forms a linker loop (Shi et al. 2004 ). The amino acid residue Cys145 and His41 at the active site forms a catalytic dyad which is essential for the catalytic function (Anand 2002) .

Currently, several synthetically derived antiviral drugs such as remdesivir, nelfinavir, favipiravir, lopinavir and ritonavir are being prominently used to curb the viral replication in SARS-CoV-2 (Mothay and Ramesh 2020; de Oliveira et al. 2020) . Moreover, the combination of these one of these antivirals with anti-malarial drugs such as ivermectin, hydroxychloroquine (Rolain et al. 2007; Yao et al. 2020 ) and few corticosteroids (Schultze et al. 2020 ) have shown promising results in the treatment of severe COVID-19 patients. Despite of these promising results these antiviral drugs have sever adverse effects (Agrawal et al. 2020; Horby et al. 2020; Wang et al. 2020) accounting to rising fatalities in COVID-19. On these grounds there is an urgent need to develop safer and sole drug against most promising target M pro of SARS-CoV-2.

Historically, natural products have made a major contributions in pharmacotherapy, especially in infectious diseases (Atanasov et al. 2021 ). In the present study, an attempt was made to identify potential of reported antiviral phytochemicals against SARS-CoV-2 M pro . The molecular dynamics (MD) guided cluster of M pro conformations were used to investigate the potential of 206 antiviral phytochemicals database through ensemble docking. The phytochemicals with better binding interactions in comparison to standard remdesivir were further subjected to extended period MD simulations to gain deeper insights of binding insights, energetic and binding affinity predictions. The results of this study may be useful in further design and development of potential candidate molecules targeted at M pro in COVID-19.

Molecular docking study was carried out in order to gain the deeper insights into the possible binding modes and estimate the binding affinity in terms of binding free energies of the library of 197 antiviral natural compounds at the binding site of SARS-CoV-2 M pro . The crystal structure of M pro (PDB ID: 6LU7) with resolution 2.16 Å was retrieved from protein databank (http:// www. rcsb. com). The crystal structures often suffer inaccuracies mostly due to steric clashes, missing atoms or bonds and the lack of information of hydrogen atom positions. In this situation, optimization of protein structure with proper positioning of hydrogen bonds is necessary (Jain 2009 ). The optimization of M pro protein structure was accomplished in following manner. Initially, the water molecules were removed from the protein structure without removing the bound peptide like ligand from the binding site. Hydrogen atoms were added in protein as well as in bound ligand. The positions of hydrogen atoms in the resulting structure were optimized with a gradient norm of 0.05 in Tinker 8 program (Rackers et al. 2018 ). The bound ligand was removed from the structure and 25 ns MD simulation was performed on the well equilibrated and solvated system. The details of MD simulations are given in Sect. 2.2. Around 2500 trajectories obtained post MD studies were clustered into 22 conformations on the basis of root mean square deviation cutoff of 2 Å.

The 2D/3D structures of phytochemicals and remdesivir were retrieved from PubChem database (https:// pubch em. ncbi. nlm. nih. gov). These structures were energy minimized in UCSF Chimera 1.8 (Pettersen et al. 2004 ) program after assigning Gasteiger charges and combination of steepest descent and conjugate gradient geometry search criteria until gradient converges to 0.05 and 0.01 respectively. Autodock vina program (Trott and Olson 2010) was used to perform the ensemble docking, with multiple conformations of protein and multiple ligands. A grid box encompassing the entire binding sites with dimension of 22 Å 3 along the x, y and z axis was set and in order to achieve the exhaustive search of chemical space the exhaustiveness 64 was selected. In-house bash script was employed to perform ensemble docking. The molecular docking results were analyzed in terms of the estimates of the binding free energy in Kcal/mol and the interactions of produced with residues at binding site.

The best five phytochemicals with lowest binding free energy and a standard drug remdesivir were subjected to the 100 ns MD simulations using Gromacs 4.5.6 (Berendsen et al. 1995). The MD simulations were performed on remote server of the Bioinformatics Resources and Applications Facility (BRAF), C-DAC, Pune. The topology of the protein was constructed with the parameters implemented in CHARMM-36 (Best et al. 2012) while the ligand topologies were generated in the CGenFF server (Vanommeslaeghe et al. 2009 ). In order to properly solvate the system a simple point charge (SPC216) water model (Nguyen et al. 2014) was used and the system was neutralized with appropriate ions such as sodium and chloride. This solvated system was subjected to unrestrained energy minimization with steepest descent criteria in order to remove the steric clashes. The system was equilibrated at constant pressure, volume and temperature of 300 K conditions for 100 picoseconds. This equilibrated system was subjected to 100 ns MD simulation with retrain on covalent bonds using LINCS algorithm (Hess et al. 1997 ) and the cut-off value of 12 Å for the long range electrostatics such as Coulomb and Lennard Jones with the Particle Mesh Ewald (PME) (Petersen 1995) method. Post MD analysis was carried out through measurement of deviations in the protein and ligand atoms as root mean square deviations (RMSD). The extent of fluctuations in the atoms of residues in protein in terms of root mean square fluctuations (RMSF) measurement was also taken in to account in estimating the stability of protein-ligand complexes. Radius of gyration (Rg) which is a root mean square distance of a collection of atoms from their common centre of mass and which provides the information of compactness of protein structure was analysed. Hydrogen bond formation is one of the most important non-bonded interactions. The hydrogen bond formation during the progress of MDS was also analysed and the number of hydrogen bonds formed and the residues important in hydrogen bond formation were also evaluated. Binding energy of protein-ligand complex formation was estimated through the combination of molecular mechanics energies with the Poisson Boltzmann surface area continuum solvation (MM-PBSA). The MM-PBSA calculations were performed with gmmpbasa program (Baker et al. 2001; Kumari et al. 2014 ).

SARS-COV-2 M pro is involved in cleaving the polyproteins pp1a and pp1ab and its inhibition could prevent the viral maturation (Liu and Wang 2020 ) and viral replication (Chang et al. 2019) . The binding site residues Thr24, Thr26, and Asn119 have been reported to play a crucial role in making key interactions (Liu and Wang 2020) . The protease inhibitors nelfinavir and lopinavir have been reported clinically beneficial in SARS-CoV infections (Li and De Clercq 2020; Li et al. 2020a) . Nelfinavir is also reported effective in in-vitro studies in SARS-COV-2 . Several natural compounds have been reported to possess antiviral activity (Hsieh et al. 2010; Islam et al. 2020; Lin et al. 2014; Orhan and Senol Deniz 2020; Perez 2003; Xian et al. 2020 ). We prepared a library of 201 potential antiviral phytochemicals. The PubChem CID's of all the selected phytochemicals are given in Table S1 (supporting information).

In the crystal structure of M pro , the bound peptide inhibitor is covalently bound to Cys145 residue which is one of the amino acid of catalytic dyad. In order to validate the docking protocol, initially this bound inhibitor was re-docked at the catalytic site of M pro with a grid box setting (see "Methods" section) large enough to accommodate the entire catalytic site as well as to hold the ligand structure within its internal dimensions. The objective of this re-docking was to only confirm whether the co-crystallized pose is accurately reproduced. The covalent interaction could not be achieved by such docking protocol and thus the covalent interaction with Cys145 residue was ignored. It was found that the docked conformer adopted almost similar conformation to the bound co-crystallized inhibitor. The conformations of bound co-crystallized inhibitor and docked inhibitor at the catalytic site are shown in Fig. 1 . The structure of SARS-CoV-2 M pro has three domains. The binding site is located in a cleft between domain I and II. The binding site has four pockets S1, S1', S2 and S4. The S1 pocket has the residues Phe140, Asn142, Glu166, His163, Leu141 and His172, while the adjacent S1' pocket which holds hydrophobic benzyl group of inhibitors is surrounded by hydrophobic residues Tyr24, Thr25, Thr26 and Leu27. The residue Cys145 in this pocket is covalently bonded to the inhibitor. The S2 pocket is surrounded by residues His41, Met49, Tyr54, Met165 and Asp187, while the S4 pocket is surrounded by residues Leu167, Phe185 and Gln192. When the optimized structure of this ligand was docked at the binding site, all these hydrogen bonds were formed except the hydrogen bonds with His164 and Gln189. Thus, the docking protocol was found adequate to regenerate the nearly similar binding pose of co-crystallized ligand and it was also validated that the grid box setting, docking and scoring function is good enough to search entire binding site encompassing all the pockets effectively.

All the optimized structures of the phytochemicals were docked at the binding site of each conformation of M pro with the same docking protocol for which we employed the inhouse bash script. The docking results were subsequently analysed for each M pro conformation. Lower the binding free energy higher is the binding affinity for M pro conformation. The docking scores (binding free energies) of only top 10 phytochemicals against each M pro conformation was taken into account. This is because, it is expected that the phytochemical with higher binding free energy possibly would have poor binding affinity. It is hypothesized that there is a complementarity between the specific binding at the binding site of SARS-CoV-2 M pro , C 2D interaction diagram for co-crystallized bound ligand and, D 2D interaction diagram for docked inhibitor site conformation of M pro and the conformation of phytochemical having the least binding free energy. For an example, rutin has the most favourable binding affinity (binding energy − 10.4 kcal mol −1 ) for the M pro conformation 7, while its binding is slightly unfavourable in conformation 10 of M pro (binding energy − 7.6 kcal mol −1 ). Each M pro conformation represents the flexibility in the binding site. Figure 2 illustrates the flexibility in the M pro conformations, especially at the binding site.

The results of ensemble docking for top 10 phytochemicals are given in Table S2 (supporting information). The top 10 phytochemicals with their corresponding best binding free energy estimates for specific conformation of M pro and key interactions are given in Table 1 .

The results of top five phytochemicals which showed most favourable binding free energy in one of the conformations of M pro are elaborated in following discussion. Rutin (3,3′,4′,5,7-pentahydroxyflavone-3-rhamnoglucoside), also known as vitamin P or rutoside or quercetin-3-O-rutinoside, is a flavonols abundantly found in plants, such as passion flower, buckwheat, tea, and apple (Ganeshpurkar and Saluja 2017) . It showed the most favourable binding affinity (binding energy − 10.4 kcal mol −1 ) on conformation 7 of M pro . It also shows favourable binding in few other conformations (conformation 5, 6 and 8). The disaccharide rutinose seems to accommodate in a polar binding pocket forming a key hydrogen bond interactions with Thr24 residue from S1' pocket and Asn142 residue from S1 pocket (Fig. 3A) . One of the hydroxyl substituent on core flavone moiety forms a hydrogen bond with Ser144, a residue present adjacent to S1 pocket. The 4H-chromen-4-one core structure of rutin (quercetin core) accommodates in S2 pocket where the pication and pi-pi interactions were observed with protonated His41 and Met165 residues respectively. The 3,4-dihydroxyl phenyl ring on core chromen-4-one was also found forming pi-pi stacking interaction with Cys145 residue. Few other hydrophobic interactions with Gly143, Glu166 and Ser46 were also formed with hydrophobic part of rutin.

Evidently the docked pose of remdesivir in the conformation 12 of Mpro (binding energy − 9.0 kcal mol −1 ) also showed the hydrogen bond with Ser144, pi-stacking interactions with His41 (Fig. 3B) . The pyrolotriazine ring of remdesivir occupies the S2 pocket, while the polar part of structure accommodates in polar S1 and S4 pockets. The decahydroisoquinoline ring of nelfinavir accommodates in the hydrophobic S4 pocket, while the polar hydroxyl group on phenyl ring forms a hydrogen bond with His172 and Phe140 residues in S1 pocket (Fig. 3C) . However, the docked pose of nelfinavir (conformation 21, binding energy − 8.5) was found slightly less favourable as compared to remdesivir.

Naringin, a flavonoid, containing a disaccharide neohesperidose, also seems to bind more favourably with the binding energy − 10.3 kcal mol −1 on the conformation 9 of M pro . The aglycone of naringin called naringenin is reported to have potential in SARS-CoV-2 (Tutunchi et al. 2020) , while naringin is also reported effective in some viruses (Zandi et al. 2011) . The docking studies revealed that the hydroxyl groups on sugar moieties forms a hydrogen bond with Thr24, Cys145 and Thr26 residues in S1' pocket (Fig. 4A) . The benzopyran ring orients in S2 pocket so as to form the hydrogen bond with His41 with hydroxyl group at 5 th position on benzopyran ring. The aromatic ring of benzopyran also forms a pi-stacking interaction with hydrophobic His41 and Met49 residues. The hydroxyl phenyl substituent occupies the S4 pocket subtly with hydrophobic pi-stacking with Pro52 residue. Naringin compared to remdesivir produces the hydrogen bond interactions with Cys145 and hydrophobic interactions with Met49. The lower value of binding free energy may be due to the favourable hydrogen bonds formed in the case of naringin. Whether these hydrogen bonds remain intact during conformational changes in the binding site could be better evaluated through MD studies.

Another flavonoid, dalpanitin, also showed favourable binding at the binding site of M pro initial crystallographic conformation with binding energy of − 10.2 kcal mol −1 . The core chromen-4-one ring of dalpanitin occupies S1 pocket forming hydrogen bonds with Glu166, His163, Ser144, Leu141 and Gly143, while the sugar substituents occupy the adjacent S2 and S4 pockets (Fig. 4B) . Interestingly, none of the two sugar substituents were found involved in hydrogen bond formation. Not only the initial conformation of M pro but also few other conformations (conformation 1, 2, 10-14, and 16-19) were favourable for binding of dalpanitin (see table SX in supporting information). However, the initial crystallographic conformation may be the best suited for better occupancy of dalpanitin presumably apparent from the binding energy of − 10.2 kcal mol −1 . Few earlier reports have shown the potencial of dalpantin as antimicrobial and antiviral phytochemical (Mohotti et al. 2020 ).

Yet another flavonoid, specifically a biflavonoid, amentoflavone, containing two apigenin units coupled through 8 and 3′ position and which is one of the constituents of famous plant ginkgo biloba, also shows favourable binding with conformation 6 of M pro . Although this phytochemical lacks sugar units, it has multiple hydroxyl groups conferring it the reasonable polarity and also possess antiviral potential (Li et al. 2019) .

It is found that these apigenin units in amentoflavone occupies S1 and S2 pockets in such a way that the hydrogen bonds are formed with the key residue His41 and Thr190; while the phenyl ring central to two chromen-4-one rings forms a pi-stacking interaction with Cys145 (Fig. 5A) . The residues His41 and Cys145 belong to the catalytic dyad and both are found participating in nonbonded interactions. It is worthwhile to see if these interactions remain intact during the dynamic conformational fluctuations during MD studies. However, the binding energy of − 10.2 kcal mol −1 clearly suggests the favourable binding affinity. The conformations 20 and 21 were also found favourable for binding of amentoflavone (binding energy of − 10.0 kcal mol −1 ).

Surprisingly, the congener of amentoflavone, hinokiflavone which is again a biflavonoid where one of the apigenin is substituted by another apigenin unit through ethereal phenoxy group was found having better binding affinity (binding energy − 10.0 kcal mol −1 for conformation 6) almost comparable to amentoflavone. The phenoxy substituted apigenin unit occupies the S1 and S2 pockets and form the hydrogen bonds with residues Cys145, Asn28, Gly143, Asn119 (Fig. 5B) . The pyran-4-one ring of this unit also forms the key pi-staking interaction with acidic anionic Glu166 residue. The residue of catalytic dyad, His41 was observed forming a pi-staking hydrophobic interaction with central phenyl ring and phenyl ring of core apigenin moiety. Few other hydrophic residues such as Met49, Pro52, Leu50, Gln189 forms the hydrophobic interactions with core apigenin moiety. These interactions almost coincide with the interactions found in amentoflavone and the structural resemblance may be responsible for similar type of interactions.

The other phytochemicals most of which are flavonoids were found to bind at the binding site of M pro but with slightly less propensity. Amongst these the phytochemical agathisflavone (binding energy − 10.0 kcal mol −1 ) form key interactions with the catalytic dyad residues and also with Thr190, Ser144, and Glu166 (Fig. 6) . Other pthytochemicals which include azadirachtin, morelloflavone, robustaflavone, and marsdenoside having the binding eneries − 8.8 to − 9.9 kcal mol −1 have complex structures. However, most of these phytochemicals bears a hydroxyl groups and carbonyl functionalities which could form the key hydrogen bonds with the residues at the binding site.

Molecular dynamics simulation (MDS) provides greater insights into binding of ligands at the binding site of the protein. MDS offers the conformational sampling over a sufficiently long period through which the flexibility in the position of ligand atom, protein backbone and side chain atoms (Jorgensen and Tirado-Rives 1996) . Further, it can give an impetus about the resultant energetic of protein, ligand and protein-ligand complex under solvated conditions which in turn provides insights of the binding free energy and binding affinity. This is especially more advantageous than the rigid docking studies where specific conformation of protein, mostly crystallographic conformation, is held fixed or rigid during docking simulation. However, the MD simulation needs to be performed for sufficiently long duration or until the system reaches the convergence or sufficient stability. In present study, 100 ns MDS studies were performed to evaluate the conformational stability of the top ranked phytochemicals viz. rutin, naringin, dalpanitin, amentoflavone, and hinokiflavone molecules and the results were compared with the standard drug remdesivir.

The measurement of protein and ligand RMSD provides good estimates of conformational stability of the system (Sargsyan et al. 2017) . The RMSD measurement involves measurement of the deviations from the starting positions of protein or ligand atoms. The MDS of sufficiently long duration captures the protein folding events, influence of loop flexibility and biding site adaptation. The RMSD analysis results for the protein showed that the complex with dalpanitin is stabilized quickly after 10 ns MDS till 80 ns and thereafter rises slightly till 90 ns and again converge to stable system towards the end of simulation time (Fig. 7) . The average, maximum and minimum RMSD values along with other estimates are given in Table 2 . In the case of dalpanitin the average RMSD of 0.232 nm suggests the stable system when it is bound to dalpanitin. Similarly, the complex of M pro with rutin stabilizes after around 50 ns and remains stable thereafter till the end of simulation with average RMSD of 0.232 nm. Possibly the number of hydrogen bonds formed with these two phytochemicals contributes in stabilization of system. The complexes with hinokiflavone, naringin, amentoflavone also stabilize towards the end of simulation. However, in the case of amentoflavone major deviations were observed during 30-50 ns. These deviations may be occurring due to binding site adaptation. The RMSD deviations in the complex with hinokiflavone are slightly higher than other phytochemicals, but are well within acceptable Fig. 7 Root mean square deviations in M pro backbone atoms and phytochemical's atoms (DAL dalpanitin, HIN hinokiflavone, AME amentoflavone, RUT rutin, REM remdesivir, NAR naringin) limits of around 0.25 nm. The larger deviations in the complex with remdesivir are clearly evident with higher average value of 0.301 nm. When the RMSD deviations in ligand atoms were analysed it was observed that hinokiflavone and remdesivir has the higher RMSD deviations than other phytochemicals. The RMSD in naringin and amentoflavone is almost stable throughout the entire MD simulation. The sharp deviations were observed in dalpanitin and rutin at around 60 and 70 ns respectively. Collectively, the results of RMSD suggest that the phytochemicals under study stabilize the M pro substantially compared to remdesivir.

The stability of complex of protein and ligands can be further evaluated in terms of Root Mean Square Fluctuations (RMSF) measurement. It is a measure of elasticity of the protein residues in terms of the fluctuations in the atoms of protein residues during the MD simulation (Fuglebakk et al. 2012 ). The RMSF for individual amino acids correlates with the trend observed in RMSD of complexes (Fig. 8) . It is evident that there are major fluctuations in the residues ranging from the residues 20-30, 40-60, 140-170, and 175-200 corresponding to S1', S2, S1 and S4 pockets respectively. The cluster of residues 20-30 belongs to the binding site residues forming the S1' pocket, while the cluster of residues 40-60 belongs to the S2 pocket of the binding site. The larger cluster of residues 160-180 represents the S1 and S4 pocket of the binding site. The other two clusters of residues 220-240 and 270-290 are farther away from the binding site. The fluctuations of the magnitude up to 0.5 nm were observed in the residues making the binding pocket. Noticeable fluctuations in the residues farther away from the binding site were observed which may be due to compensating the binding site adaptability. The larger fluctuations are clearly evident in the case of remdesivir and naringin. The least fluctuations are observed with dalpanitin and rutin, and consequently the RMSD for these phytochemicals is also lower than other ligands.

The analysis of radius of gyration (Rg) is another criteria to judge the root mean square distance of a collection of atoms from their common center of mass which can provide an insight of the overall compactness of the system (Monajjemi and Oliaey 2009) . The average Rg values for the SARS-CoV-2 M pro protein complexed with hinokiflavone, naringin, rutin, amentoflavone, dalpanitin is around 2.2 nm, while the average Rg value for remdesivir is 2.3 nm (Fig. 9) . The results suggest that the complex of M pro with phytochemicals remains compact throughout the entire simulation compared to remdesivir. The subtle differences in the Rg values suggest that all the systems are almost equally adopts the compact protein structures and this may be in part due to major contributions from the α-helices in the protein structure. The binding affinity of ligands majorly depends upon the non-bonded interactions such as hydrogen bonds, hydrophobic interactions and ionic interactions at the binding site (Du et al. 2016) . More the number of hydrogen bonds formed and longer the life time of hydrogen bonds are indicative of stronger binding affinity. The maximum eight numbers of hydrogen bonds were formed with rutin of which minimum 4 hydrogen bonds are consistently formed throughout the simulation (Fig. 10) . In the case of dalpanitin, maximum 7 and 5 hydrogen bonds were found consistently formed. In the case of three ligands, naringin, amentoflavone, and remdesivir maximum 6 hydrogen bonds out of which minimum 3 are consistently found forming. Hinokiflavone forms maximum 4 hydrogen bonds out of which only 2 seems to be consistently forming during simulation period.

The hydrogen bond formation phenomenon was further investigated to understand which residues at the binding site are contributing in hydrogen bond formation. For this the hydrogen bond occupancy percentage was analysed. The cut off radius of 0.35 nm and angle cut off of 120º criteria was employed while calculating percent occupancy of hydrogen bond (Torshin et al. 2002) . The hydrogen bond occupancy results are shown in (Fig. 11) . The standard drug remdesivir was found forming major ensemble of hydrogen bonds at S1 pocket with residue Glu166 and at S2 pocket with Asp187 both having the highest percent occupancy suggesting major contribution in stabilizing the system. In case of amentoflavone the hydrogen atoms from the hydroxyl group of central phenyl ring form the hydrogen bond consistently with His164 at the S1 pocket with highest occupancy.

Interestingly, it was found that dalpanitin forms the hydrogen bonds with highest occupancy with S4 pocket residues Arg188 and Gln192. The S1' residue Thr26 and S2 pocket residue from catalytic dyad His41 forms hydrogen bonds with highest occupancy in the case of hinokiflavone. Naringin forms the hydrogen bonds with highest occupencies with S1 pocket residues Cys145 and Ser144 and S2 pocket residues Ser46 and His164. Rutin was found forming hydrogen bonds with S1 pocket residues Ser144 and Leu141 and S4 pocket residue Arg188. The results of hydrogen bond analysis and particularly estimation of maximum number of hydrogen bonds formed and percent occupancy of hydrogen bonds suggest that the phytochemicals hinokiflavone, amentoflavone and naringin may have more favourable binding interactions and consequently may have better binding affinity.

The accurate estimates of binding affinity can be obtained from MM-PBSA calculations. In present work, the g_ mmpbsa program was used to calculate the van der Waal energy, electrostatic energy, polar solvation energy, SASA energy and binding energy (Table 3) .

In MM-PBSA calculations the non-bonded interactions such as van der Waal energy and electrostatic energy measured in terms of Coulomb and Lennard-Jones (LJ) potential functions respectively has the major influence on the binding free energy (ΔG binding ) estimate. The MM-PBSA results showed that the van der Waal energy is lower for phytochemicals hinokiflavone, amentoflavone, rutin and naringin. However, due to larger polar solvation energy for rutin its ΔG binding is slightly higher than other phytochemicals. The MM-PBSA calculations suggests that the phytochemicals hinokiflavone, amentoflavone and naringin with ΔG binding of − 80.8, − 73.7, and − 68.9 kJ.mol −1 respectively have the best binding affinity for M pro compared to remdesivir. These phytochemicals could be further evaluated for their potential use in management of SARS-CoV-2 viral replication devising in-vitro and in-vivo testing. 

Currently, Covid-19 has emerged as a global pandemic with enormous global threat to public health and no approved specific drug therapy exist to treat the disease. This study provides a comprehensive overview of the structural characteristics of the SARS-COV-2 proteases and its potential inhibitors. The libraries of phytochemicals with potential antiviral activities were investigated as possible SARS-COV-2 M pro inhibitors through the molecular modelling studies. The results of docking studies indicated that dalpanitin, hinokiflavone, amentoflavone, rutin, and naringin have the lower binding free energy against SARS-COV-2 M pro . These molecules and the reference drug remdesivir were further analysed through 100 ns molecular dynamics simulations in order to obtain more precise and accurate estimates of binding free energy and binding affinity. The MDS study revealed the importance of key interactions at the binding site. Particularly, His41 residue at S2 pocket and few residues at S1 pocket such as Glu166, His164 in eliciting the prominent binding interaction at the binding site of SARS-COV-2 M pro . The MM-PBSA calculations, which were carried out to accurately estimate the ΔG binding which revealed that the phytochemicals hinokiflavone, amentoflavone and naringin with ΔG binding of − 80.8, − 73.7, and − 68.9 kJ mol −1 respectively have better binding affinities than the standard drug remdesivir. These phytochemical could turn out to be the promising lead compounds in further development of potential SARS-COV-2 M pro inhibitors. 

Favipiravir: a new and emerging antiviral option in COVID-19

Structure of coronavirus main proteinase reveals combination of a chymotrypsin fold with an extra alpha-helical domain

Natural products in drug discovery: advances and opportunities

Covid-19-the search for effective therapy

Electrostatics of nanosystems: Application to microtubules and the ribosome

Remdesivir for the treatment of Covid-19-final report

GROMACS: a message-passing parallel molecular dynamics implementation

Optimization of the additive CHARMM all-atom protein force field targeting improved sampling of the backbone ϕ, ψ and side-chain χ1 and χ2 dihedral angles

Antiviral drug discovery: norovirus proteases and development of inhibitors

Repurposing approved drugs as inhibitors of SARS-CoV-2 S-protein from molecular modeling and virtual screening

Insights into protein-ligand interactions: mechanisms, models, and methods

SARS-CoV-2 infections in children and young people

Measuring and comparing structural fluctuation patterns in large protein datasets

The pharmacological potential of rutin

LINCS: a linear constraint solver for molecular simulations

Lopinavir-ritonavir in patients admitted to hospital with COVID-19 (RECOVERY): a randomised, controlled, open-label, platform trial

Synergistic antiviral effect of Galanthus nivalis agglutinin and nelfinavir against feline coronavirus

Natural products and their derivatives against coronavirus: a review of the non-clinical and pre-clinical data

Effects of protein conformation in docking: improved pose prediction through protein pocket adaptation

Structure of M pro from SARS-CoV-2 and discovery of its inhibitors

Monte Carlo vs molecular dynamics for conformational sampling

Open Source Drug Discovery Consortium (2014) g_mmpbsa-a GROMACS tool for high-throughput MM-PBSA calculations

Therapeutic options for the 2019 novel coronavirus (2019-nCoV)

Amentoflavone inhibits HSV-1 and ACV-resistant strain infection by suppressing viral early infection

The impact of mutations in SARS-CoV-2 spike on viral infectivity and antigenicity

The epidemic of 2019-novel-coronavirus (2019-nCoV) pneumonia and insights for emerging infectious diseases in the future

Antiviral natural products and herbal medicines

Potential inhibitors against 2019-nCoV coronavirus M protease from clinically approved medicines

Structural elucidation of SARS-CoV-2 vital proteins: computational methods reveal potential drug candidates against main protease, Nsp12 polymerase and Nsp13 helicase

Screening for bioactive secondary metabolites in Sri Lankan medicinal plants by microfractionation and targeted isolation of antimicrobial flavonoids from Derris scandens

Energy study at different temperatures for acetylcholine receptor protein in gas phase by Monte Carlo, molecular and Langevin dynamics simulations

Binding site analysis of potential protease inhibitors of COVID-19 using AutoDock

Effects of water models on binding affinity: evidence from all-atom simulation of binding of tamiflu to A/H5N1 neuraminidase

Potential anti-COVID-19 agents, cepharanthine and nelfinavir, and their usage for combination treatment. iScience

Natural products as potential leads against coronaviruses: could they be encouraging structural models against SARS-CoV-2?

Emerging SARS-CoV-2 mutation hot spots include a novel RNA-dependent-RNA polymerase variant

Antiviral activity of compounds isolated from plants

Accuracy and efficiency of the particle mesh Ewald method

UCSF Chimera-a visualization system for exploratory research and analysis

Drug for corona virus: a systematic review

Tinker 8: software tools for molecular design

Presenting characteristics, comorbidities, and outcomes among 5700 patients hospitalized with COVID-19 in the

Recycling of chloroquine and its hydroxyl analogue to face bacterial, fungal and viral infections in the 21st century

How molecular size impacts RMSD applications in molecular dynamics simulations

Risk of COVID-19-related death among patients with chronic obstructive pulmonary disease or asthma prescribed inhaled corticosteroids: an observational cohort study using the OpenSAFELY platform

Insights into RNA synthesis, capping, and proofreading mechanisms of SARS-coronavirus

Dissection study on the severe acute respiratory syndrome 3C-like protease reveals the critical role of the extra domain in dimerization of the enzyme: defining the extra domain as a new target for design of highly specific protease inhibitors

Temperature-sensitive mutants and revertants in the coronavirus nonstructural protein 5 protease (3CLpro) define residues involved in long-distance communication and regulation of protease activity

Geometric criteria of hydrogen bonds in proteins and identification of `bifurcated' hydrogen bonds

AutoDock Vina: improving the speed and accuracy of docking with a new scoring function, efficient optimization, and multithreading

Naringenin, a flavanone with antiviral and anti-inflammatory effects: a promising treatment strategy against COVID-19

The SARS-CoV-2 main protease as drug target

CHARMM general force field: A force field for druglike molecules compatible with the CHARMM all-atom additive biological force fields

Remdesivir in adults with severe COVID-19: a randomised, double-blind, placebo-controlled, multicentre trial

int/? gclid= EAIaI QobCh MI_ 9nyhu yh7AI VlQ4r Ch0Fd APVEA AYASA AEgL5 RfD_ BwE

Bioactive natural compounds against human coronaviruses: a review and perspective

Nelfinavir was predicted to be a potential inhibitor of 2019-nCov main protease by an integrative approach combining homology modelling, molecular docking and binding free energy calculation

Antiviral activity of four types of bioflavonoid against dengue virus type-2

A pneumonia outbreak associated with a new coronavirus of probable bat origin

Publisher's Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations

The author SDS, RBP, and VPZ are thankful to Prof. M. N. Navale, Founder President, Sinhgad Technical Edication Society, Pune for encouragement and support.Funding This research did not receive any specific grant from funding agencies in the public, commercial, or not-for-profit sectors.

The authors declare that they have no competing interests.