Computational assessment of select antiviral phytochemicals as potential SARS-Cov-2 main protease inhibitors: molecular dynamics guided ensemble docking and extended molecular dynamics.

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 (Mpro), 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 Mpro 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 Mpro stabilize with well accepted RMSD. Amongst these phytochemicals, hinokiflavove, amentoflavone and naringin were found having better binding affinity with ΔGbinging 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. SUPPLEMENTARY INFORMATION: The online version contains supplementary material available at 10.1007/s40203-021-00107-9.

Introduction

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 (COVID-19) (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 (Zhou et al. 2020) 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 [Mpro] or 3C-like protease [3CLpro]) and papain-like protease [PLpro] cleaves these polyproteins into 16 non-structural proteins (NSP1-16) (Prajapat et al. 2020). The Mpro is highly conserved in its sequence throughout coronaviruses (Mirza and Froeyen 2020). The role of Mpro is established in proteolytic cleavage at 11 sites at C-terminus of viral polyprotein 1ab (Stobart et al. 2012). Hence, this unique feature makes Mpro as high profile and best characterized target for anti-SARS therapy (Anand 2002). Monomeric form of Mpro is inactive, while homodimeric is active in catalytic regulation. The structure of Mpro 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 Mpro 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 Mpro. The molecular dynamics (MD) guided cluster of Mpro 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 Mpro in COVID-19.

Methods

Molecular docking with ensemble docking approach

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 Mpro. The crystal structure of Mpro (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 Mpro 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://pubchem.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.

Molecular dynamics simulation studies and MM-PBSA calculations

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).

Results and discussion

Molecular docking studies

SARS-COV-2 Mpro 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 (Xu et al. 2020). 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 Mpro, the bound peptide inhibitor is covalently bound to Cys145 residue (Jin et al. 2020) 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 Mpro 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 Mpro 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.

Fig. 1

Binding pocket of SARS-CoV-2 Mpro. A Binding pose of co-crystallized inhibitor showing different pockets, B conformations of bound co-crystallized inhibitor (red) and docked inhibitor (yellow) at the binding site of SARS-CoV-2 Mpro, C 2D interaction diagram for co-crystallized bound ligand and, D 2D interaction diagram for docked inhibitor

All the optimized structures of the phytochemicals were docked at the binding site of each conformation of Mpro with the same docking protocol for which we employed the in-house bash script. The docking results were subsequently analysed for each Mpro conformation. Lower the binding free energy higher is the binding affinity for Mpro conformation. The docking scores (binding free energies) of only top 10 phytochemicals against each Mpro 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 site conformation of Mpro 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 Mpro conformation 7, while its binding is slightly unfavourable in conformation 10 of Mpro (binding energy − 7.6 kcal mol−1). Each Mpro conformation represents the flexibility in the binding site. Figure 2 illustrates the flexibility in the Mpro conformations, especially at the binding site.

Fig. 2

Aligned conformations of Mpro showing the binding site and orientation of side chains of few residues

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 Mpro and key interactions are given in Table 1.

Table 1

Ensemble docking results for top 10 phytochemicals

Sr. no.	Phytochemical (Pubchem CID)	The Mpro conformation number	Binding free energy (ΔG) kcal mol−1	Interactions at the binding site
				H-bond interactions	Hydrophobic interactions
1	Rutin (5280805)	Conf_7	− 10.4	Thr24, Ser144, Asn142	Ser46, Cys145, Glu166, His41, Met165, Gly143
2	Naringin (442428)	Conf_9	− 10.3	Thr24, Cys145, Thr26, His41	Thr25, Leu27, Ser144, Met49, Gln189, Leu50, His164, Thr190, Met165
3	Dalpanitin (44257304)	Conf_0	− 10.2	Gly143, Leu141, Ser144, His163, Glu166	Cys145, Asn142, His172, Met165, Thr25, Met49, Leu27, His41, Gln192
4	Amentoflavone (5281600)	Conf_6	− 10.2	His41, Thr190	Cys145, Pro168, Leu27, Thr25, Thr26, Ser144, His163, Glu166, Gln189
5	Hinokiflavone (5281627)	Conf_9	− 10.0	Asn119, Gly143, Asn28, Cys145	His41, Met49, Leu27, Thr26, Thr25, Gln189, Met165, Glu166, His164
6	Agathisflavone (5281599)	Conf_21	− 10.0	Thr190, Cys145, Glu166	Met165, His41, Leu141, Thr25, Leu27, His163, Gln189, Met49, His164
7	Azadirachtin (52813032)	Conf_20	− 9.9	Gln189, Gly143, Cys145, Ser144, Glu166	Cys44, His41, Leu27, Thr26, Thr25, His163, Met165, Leu167, Met49
8	Morelloflavone (5464454)	Conf_22	− 9.4	Glu166, His164, Tyr54, Arg188	Met165, His163, Cys145, His41, Asp187, Thr190, Leu167, Pro52, Gln189
9	Robustaflavone (5281694)	Conf_22	− 9.2	His163	Met49, His41, His172, Glu168, Cys145, His164, Tyr54, Thr25, Gly143
10	Marsdenoside (11343196)	Conf_12	− 8.8	Asn142, Thr24, Gly143	His41, Pro168, Met165, Leu27, Thr26, Cys44, Ser46, Cys145, Met49
11	Remdesivir (121304016)	Conf_12 Conf_15	− 9.0	Glu166, Cys145, Ser144	His41, Cys44, Met49, Leu27, Thr25, Gln189, Thr190, His163, Asn142
12	Nelfinavir (64143)	Conf_21	− 8.5	Cys145, Ser144, Gly143, His172, Phe140	His41, Glu166, Met165, Pro168, Thr190, Tyr54, His164, His163, Leu141

The results of top five phytochemicals which showed most favourable binding free energy in one of the conformations of Mpro 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 Mpro. 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 pi-cation 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.

Fig. 3

Docked poses of ligands in the corresponding most favourable conformation of Mpro. A Rutin, B remdesivir, and C nelfinavir

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 Mpro. 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 5th 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.

Fig. 4

Docked poses of A naringin and B dalpanitin

Another flavonoid, dalpanitin, also showed favourable binding at the binding site of Mpro 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 Mpro 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 Mpro. 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 non-bonded 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).

Fig. 5

Docked poses of A amentoflavone and B hinokiflavone

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 Mpro 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.

Fig. 6

2D-interaction diagrams for docked poses of other phytochemicals

Molecular dynamics studies and MM-PBSA calculations

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 Mpro 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 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 Mpro substantially compared to remdesivir.

Fig. 7

Root mean square deviations in Mpro backbone atoms and phytochemical’s atoms (DAL dalpanitin, HIN hinokiflavone, AME amentoflavone, RUT rutin, REM remdesivir, NAR naringin)

Table 2

Average, minimum and maximum values of few MDS evaluations

Parameter	Dalpanitin	Hinoki-flavone	Amento-flavone	Rutin	Naringin	Remdesivir
Protein backbone atoms RMSD (nm)
Average	0.217	0.280	0.240	0.232	0.253	0.301
Maximum	0.342	0.391	0.368	0.344	0.409	0.430
Minimum	0.045	0.055	0.045	0.038	0.041	0.040
Atoms of ligand RMSD (nm)
Average	0.164	0.197	0.109	0.222	0.144	0.253
Maximum	0.263	0.307	0.219	0.322	0.228	0.362
Minimum	0.018	0.022	0.022	0.041	0.027	0.090
RMSF (nm)
Average	0.134	0.144	0.164	0.145	0.186	0.191
Maximum	0.551	0.587	0.565	0.339	0.438	0.537
Minimum	0.045	0.053	0.059	0.053	0.059	0.060
Radius of gyration (nm)
Average	2.223	2.206	2.265	2.228	2.207	2.300
Maximum	2.278	2.240	2.310	2.272	2.264	2.370
Minimum	2.179	2.165	2.217	2.190	2.152	2.228

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.

Fig. 8

Root mean square fluctuations in the side chain atoms of residues when complexed with phytochemicals. (S1’, S2, S1, S4 represents the respective pockets in Mpro)

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 Mpro 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 Mpro 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.

Fig. 9

Radius of gyration representing the compactness of 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.

Fig. 10

Number of hydrogen bonds formed between the molecules and the residues at the binding site of SARS-CoV-2 Mpro protein

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.

Fig. 11

Hydrogen bond occupancy percentage for the molecules. The respective atoms of each molecule involved in hydrogen bond formation are shown in the 3D structures of the molecules

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).

Table 3

Estimates of different interaction energies of MM-PBSA calculations

Compound ID	van der Waal energy (kJ mol−1)	Electrostattic energy (kJ mol−1)	Polar solvation energy (kJ mol−1)	SASA energy (kJ mol−1)	Binding energy (kJ mol−1)
Dalpanitin	− 116.383 (12.852)	− 27.793 (20.341)	114.910 (28.772)	− 15.018 (2.140)	− 44.284 (9.925)
Hinokiflavone	− 190.126 (13.180)	− 39.192 (11.177)	169.929 (21.726)	− 21.479 (1.705)	− 80.868 (16.646)
Amentoflavone	− 183.032 (12.200)	− 52.456 (9.619)	182.168 (17.116)	− 20.471 (1.210)	− 73.791 (8.534)
Rutin	− 197.796 (15.983)	− 66.631 (23.870)	236.312 (21.347)	− 20.981 (1.238)	− 49.095 (14.757)
Naringin	− 190.726 (7.356)	− 18.878 (8.448)	161.540 (12.986)	− 20.914 (0.861)	− 68.978 (8.328)
Remdesivir;	− 161.493 (17.564)	− 64.115 (6.902)	189.853 (22.297)	− 20.018 (1.379)	− 55.772 (12.978)

Standard deviations are given in parentheses

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 (ΔGbinding) 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 ΔGbinding is slightly higher than other phytochemicals. The MM-PBSA calculations suggests that the phytochemicals hinokiflavone, amentoflavone and naringin with ΔGbinding of − 80.8, − 73.7, and − 68.9 kJ.mol−1 respectively have the best binding affinity for Mpro 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.

Conclusion

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 Mproinhibitors 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 Mpro. 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 Mpro. The MM-PBSA calculations, which were carried out to accurately estimate the ΔGbinding which revealed that the phytochemicals hinokiflavone, amentoflavone and naringin with ΔGbinding 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 Mpro inhibitors.

Below is the link to the electronic supplementary material.

Supplementary file1 (DOCX 32 kb)