Identification of potential plant-based inhibitor against viral proteases of SARS-CoV-2 through molecular docking, MM-PBSA binding energy calculations and molecular dynamics simulation.

The Coronavirus disease 2019 (COVID-19), caused by the novel coronavirus, SARS-CoV-2, has recently emerged as a pandemic. Here, an attempt has been made through in-silico high throughput screening to explore the antiviral compounds from traditionally used plants for antiviral treatments in India namely, Tea, Neem and Turmeric, as potential inhibitors of two widely studied viral proteases, main protease (Mpro) and papain-like protease (PLpro) of the SARS-CoV-2. Molecular docking study using BIOVIA Discovery Studio 2018 revealed, (-)-epicatechin-3-O-gallate (ECG), a tea polyphenol has a binding affinity toward both the selected receptors, with the lowest CDocker energy - 46.22 kcal mol-1 for SARS-CoV-2 Mpro and CDocker energy - 44.72 kcal mol-1 for SARS-CoV-2 PLpro, respectively. The SARS-CoV-2 Mpro complexed with (-)-epicatechin-3-O-gallate, which had shown the best binding affinity was subjected to molecular dynamics simulations to validate its binding affinity, during which, the root-mean-square-deviation values of SARS-CoV-2 Mpro-Co-crystal ligand (N3) and SARS-CoV-2 Mpro- (-)-epicatechin-3-O-gallate systems were found to be more stable than SARS-CoV-2 Mpro system. Further, (-)-epicatechin-3-O-gallate was subjected to QSAR analysis which predicted IC50 of 0.3281 nM against SARS-CoV-2 Mpro. Overall, (-)-epicatechin-3-O-gallate showed a potential binding affinity with SARS-CoV-2 Mpro and could be proposed as a potential natural compound for COVID-19 treatment.

Introduction

Since initially being reported as pneumonia of unknown etiology in Wuhan, China, COVID-2019 has almost crippled the health care system globally and the world has been brought to a standstill. The outbreak caused by a novel coronavirus (CoV), SARS-CoV-2 was declared a Public Health Emergency of International Concern (PHEIC) on January 30, 2020 and a pandemic on March 11, 2020 [1]. With 216 countries, areas or territories affected by the pandemic, cases are escalating with every passing day [2]. According to the World Health Organization (WHO) statistics, as of December 11, 2020, 68,845,368 confirmed cases and 1,570,304 deaths have been reported worldwide [3]. Currently, there is yet no effective treatment or vaccine approved or developed against this SARS-CoV-2 and hence, there is an urgent need for development for effective drugs against this virus [4, 5].

Coronaviruses are positive single-stranded RNA viruses of approximately 30 kb genome size [6]. For effective replication of these viruses, the large viral polyprotein translated using the host cell machinery needs to be cleaved by main protease (Mpro) along with another papain-like protease(s) (PLpro) of virus origin [7]. As inhibition of main protease (Mpro) and PLpro would stop viral replication, these enzymes have become the best-characterized drug targets for the coronaviruses [8]. Presence of significant antagonistic effect of host cell immune system, PLpro has become a popular target for its inhibitors [9]. Besides, the Mpro of SARS-CoV-2 shows 96% similarity with SARS-CoV-1 and also remains as a promising target for the development of new antiviral drugs (Supplementary Fig. 1) [2, 10, 11].

Plant products have been used for generations in Chinese traditional and Indian Ayurvedic medicines for treating viral diseases and the safety of their usage has already been established [12, 13]. Herbal treatment was also found to be effective during the outbreak of severe acute respiratory syndromes (SARS) in 2003 [14]. Therefore, readily accessible natural compounds with minimum side effects are expected to impede the current outbreak and provide the first line of cure and prevention options [13]. In Indian Ayurveda, plants like Tea (Camellia sinensis), Neem (Azadirachta indica) and Turmeric (Curcuma longa) are being widely used for the treatment of viral diseases since ancient era, thereby indicating the presence of active antiviral entities [15-20]. Additionally, these plants are also being used widely for the treatment of viral diseases by the people of different communities of North Eastern part of India traditionally.

The urgency of novel drugs for controlling SARS-CoV-2 generates high demand for computational drug discovery approaches. Virtual screening of natural compounds may provide an alternative approach to identify the potential drug candidate to combat pandemic infectious diseases like novel coronavirus disease, particularly where specific FDA-approved drugs are not available [1, 21, 22]. These approaches offer a great opportunity of testing the hypothesis of potential drug effect of the natural compounds like the other synthetic counterparts. Based on the traditional use as well as on Indian Ayurveda, this study intends to identify putative candidate compound(s) as a potential therapeutic agent for COVID-19 using computational approaches by screening a large number of natural compounds from Tea, Neem and Turmeric functioning as potential inhibitors of the Mpro of SARS-CoV-2.

Materials and methods

Processing of the crystal structures of the selected receptors

In the present study, the three-dimensional structure of PLpro was retrieved from RCSB PDB (PDB ID: 6WX4) which was co-crystallized with peptide ligand Vir251. On the other hand, two published structures of SARS-CoV-2 Mpro, co-crystallized with ligands were also downloaded from RCSB Protein Data Bank (PDB), [10, 23]. The PDB IDs of the structures for SARS-CoV-2 Mpro were: (1) 6LU7, SARS-CoV-2 Mpro co-crystallized with N3 (a Michael-acceptor inhibitor) and (2) 6Y2G, SARS-CoV-2 Mpro co-crystallized with 13b (α-ketoamide inhibitor) (Supplementary Table 1). To understand the binding pattern of SARS-CoV-2 Mpro with its ligands, the binding site of N3 in 6LU7 was analyzed and compared with the binding site of 13b in 6Y2G. Both these ligands were found to bind at a similar binding site of the protein. The same approach was applied for identification of the binding site in case of the SARS-CoV-2 PLpro also, in the crystal structure having PDB ID 6WX4. The structures of SARS-CoV-2 proteases, 6LU7 and 6WX4 were further processed for structure-based analysis. Prior to molecular docking analysis, the target proteins were cleaned by deleting alternate conformations, adjusting the terminal residues and correcting the bond orders. The cleaning process was followed by preparation of the protein for docking analysis. Proteins were prepared by “Prepare Protein” module of “Macromolecules” tool of BIOVIA Discovery Studio (DS) 2018. During the course of protein preparation, the water molecules were removed keeping the target enzymes intact together with their co-crystal ligands. Energy minimization of the protein was performed using CHARMm-based smart energy minimizer of BIOVIA DS 2018 at maximum of 200 steps and 0.1 kcal mol−1 RMSD gradient [24, 25] for making the structure ready for further analysis.

Processing of the plant-based and synthetic small molecules

A list of compounds from Tea, Neem and Turmeric were collated from published literatures and phytochemical databases. The 2D structures of the compounds were retrieved from the PubChem Compound Database of NCBI and a library was prepared. The compound library was passed through Lipinski’s Drug Discovery filter [26, 27] and structures of 1295 shortlisted phytochemicals were converted to 3D counterparts using the ChemOffice package [28, 29].

The two ligands, namely, N3 and 13b bound to SARS-CoV-2 Mpro in the crystal structures 6LU7 and 6Y2G, respectively, were taken as controls for in silico binding experiment. These ligands have been previously reported as inhibitors of SARS-CoV-2 Mpro [10, 23]. Similarly, in the case of 6WX4, the inhibitor Vir251 was used as control during the experiment.

The 3D structures of the ligands were further prepared for docking analysis using “Small Molecule” tool of BIOVIA DS 2018 followed by energy minimization of the compounds with the help of “Full Minimization” module of same tool. The minimization of energy was done using CHARMm-based energy minimizer with 2000 steps of Steepest Descent and Conjugate Gradient algorithm with 0.01 kcal mol−1 RMSD gradient.

Molecular docking

CHARMm-based smart minimizer method was used to minimize the energy of the target proteins as well as ligands and prepared for docking study using CDocker of BIOVIA Discovery Studio (DS) 2018 [25]. The binding site of N3 and 13b (α-ketoamide) to SARS-CoV-2 Mpro was used as the active site for molecular docking study. The coordinates used for binding site of Mpro were X: − 10.897; Y: 13.066; Z: 68.557 with radius 12.738 Å. Similarly, the binding site of Vir251 in 6WX4 was considered as the binding site in case of SARS-CoV-2 PLpro. The coordinates for PLpro were X: 8.904; Y: − 27.443; Z: − 37.926 with radius 10.323 Å. After molecular docking, the binding energies for ligand-receptor interactions were calculated with the help of ‘Calculate Binding Energy’ protocol of DS 2018. Lower value of CDocker energy gives the best binding affinity of the ligand to the receptor protein [25]. The best binding complexes were then processed for molecular dynamics (MD) simulations to study the stability of the protein–ligand complexes which could mimic the conditions of in vitro and in vivo experiments [31].

Molecular dynamics (MD) simulation and pharmacokinetic analysis

All MD simulations were conducted using GROMACS 2018.2 package with a united-atom force field (GROMOS96 43a1) [32-38] for a time scale of 100 ns each at 310 K (37 °C) on a DELL PowerEdge R740 Rack server machine comprising 40 physical processor cores within two second Generation Intel® Xeon®Gold Scalable processors accelerated by an NVIDIA®Tesla®V100 Tensor Core graphic processor unit (GPU).

For performing the MD Simulations, topology files of the small molecules N3 (co-crystal ligand) and (−)-epicatechin-3-O-gallate (ECG) were obtained from PRODRG [39, 40]. All defined systems were solvated with extended-SPC explicit solvent water model. Four Na+ ions were added to neutralize the system before energy minimization and position restraint MDs (NVT and NPT for 100 ps each). A water box of 5 Å from the surface of the protein was created for all three systems. The systems were neutralized with counter-ions and energy minimization was performed using steepest descent for 50,000 steps. For all three systems (SARS-CoV-2 Mpro, SARS-CoV-2 Mpro-N3, and SARS-CoV-2 Mpro-ECG), the protein backbone was frozen and solvent molecules with counter-ions were allowed to move for two 100 ps position restrained equilibration MD runs. All simulations were performed under periodic boundary conditions with NVT followed by NPT ensemble. During the position restraint MD runs, V-rescale and Berendsen's coupling algorithms were used to keep the temperature (310 K) and pressure (1 bar) constant, respectively. Finally, 100 ns of production MD runs were performed allowing all molecules to move in all directions according to a classical Newtonian leap-frog MD integrator. For all the systems, the pressure was maintained at 1 bar by isotropic pressure coupling in x, y and z components to a Parrinello–Rahman barostat with the time constant τ = 2.0 ps and compressibility of 4.5 × 10 − 5 bar−1 in all three dimensions. The electrostatic interactions were calculated by the PME algorithm, with coulomb cutoff of 12 Å and interpolation order of 4 within a grid spacing of 0.16 nm. The time steps for the simulations were 2 fs, and the coordinates were stored every 5 ps. The Van der Waals (vdw) forces were treated using a cutoff of 12 Å with force-switch vdw-modifier and rvdw-switch 1. LINCS algorithm was used to constrain all bond lengths. GNU image manipulation tool 2.8 (GIMP) was used to label and reconstruct the assembled diagrams. Protein–ligand interactions were evaluated by LigPlot+ [41]. Further, the pharmacokinetic analysis of the best binding molecule was done using pkCSM [42, 43].

MM-PBSA calculations

Binding free energies (ΔGBind) of the protein–ligand complexes were calculated from molecular mechanics Poisson–Boltzmann surface area (MM-PBSA) evaluated by adaptive Poisson–Boltzmann solver 3.0 (APBS 3.0) using g_mmpbsa package [44, 45]. The MM-PBSA approach is one of the most widely used methods to compute interaction energies among the biomolecular complexes. Together with MD simulation, MM-PBSA can decode significant conformational fluctuations and entropic contributions to the binding energy. In general, the binding free energy (GBind) between a protein and a ligand in solvent can be expressed as:where GComplex is the total free energy of the protein − ligand complex and GProtein and GLigand are total free energies of the separated protein and ligand in solvent, respectively. The same equation can be used to calculate other types of biomolecular interactions such as protein − protein, protein − DNA complexes, etc.

HOMO–LUMO energy calculation

Density functional theory (DFT) analysis was performed in order to determine the electrostatic properties of N3 and ECG. The frontier orbital energies, i.e., highest occupied molecular orbital (HOMO) and lowest unoccupied molecular energy (LUMO) were calculated using Becke’s three-parameter exchange functionals (B3LYP) module of Dmol3 package of BIOVIA DS 2018. The HOMO and LUMO energy profile describe the electron donor and acceptor properties of any compounds which ultimately gives the idea on the reactivity of the compounds to their receptors [46]. In this study, by analyzing the HOMO and LUMO energy profile, the energy gap of N3 and ECG was determined.

Quantitative structure–activity relationship (QSAR) analysis

To study the inhibitory potentiality of ECG, the best binding molecule, a quantitative structure–activity relationship (QSAR) analysis was performed. For the QSAR Study, 21 reported inhibitors were taken from ChEMBL BioAssay (CHEMBL854506) database [47, 48]. The molecules were loaded in BIOVIA Discovery Studio 2019 and their energy minimization was done using CHARMm. Total Energy was predicted using the DMol3 module of the Discovery Studio and other descriptors were also calculated [49]. Activities were predicted for each of the molecule using the generated QSAR equation, validated using the test data set and plotted against their experimental activities. As the R2 value for the QSAR Equation was reasonably high, the equation was used to predict the activity of ECG. Similarly, QSAR analysis of N3 was also performed to compare inhibitory potentiality of ECG with N3.

Results and discussion

The plant molecules with the best binding affinity to SARS-CoV-2 Mpro and SARS-CoV-2 PLpro were selected based on lowest CDocker energy values. The result revealed that, among all the plant-based compounds used in the study, a tea polyphenol, ECG showed the lowest CDocker energy for SARS-CoV-2 Mpro, that is, − 46.22 kcal mol−1, followed by Phloretin (− 36.36 kcal mol−1) of tea (C. sinensis). Interestingly, in case of SARS-CoV-2 PLpro, the same polyphenol, i.e., ECG showed the lowest CDocker energy with value − 44.72 kcal mol−1 followed by Nordihydroguaiaretic acid (− 35.37 kcal mol−1) of tea (C. sinensis). Based on the CDocker energy, the best compounds having good binding affinity to both the selected receptors are tabulated in Table 1. Since ECG showed the best binding affinity with SARS-CoV-2 Mpro, this complex was considered for further MD simulation.

Table 1

Best binding plant-based compounds (best four) selected based on lower CDocker energy

SARS-CoV-2 Mpro (PDB ID: 6LU7)			SARS-CoV-2 PLpro (PDB ID: 6WX4)
[Compound Name] and [PubChem CID]	2D Structure	CDocker energy (kcal mol−1)	[Compound Name] and [PubChem CID]		2D Structure	CDocker energy (kcal mol−1)
[(−)-epicatechin-3-O-gallate (ECG)] [107905]		− 46.2268	[(−)-epicatechin-3-O-gallate (ECG)] [107905]			− 44.7288
[Phloretin] [4788]		− 36.3649	[Nordihydroguaiaretic acid] [4534]			− 35.3797
[Myrecetin] [5281672]		− 35.091	[Propyl gallate] [4947]			− 35.1536
[Epicatechin] [72276]		− 33.3683	[Phloretin] [4788]			− 33.0821

To compare the binding of ECG with SARS-CoV-2 Mpro with the control, molecular docking was also performed between control molecules and the target proteins. Because ECG showed the highest binding affinity with SARS-CoV-2 Mpro, co-crystal ligands N3 and 13b were allowed to dock with SARS-CoV-2 Mpro. Between the controls, N3 (Supplementary Fig. 2) showed the lowest CDocker energy (− 91.37 kcal mol−1) after docking with SARS-CoV-2 Mpro. Surprisingly, the other control, that is, the ligand 13b showed CDocker energy − 40.82 kcal mol−1, which was lower than the CDocker energy of ECG. The co-crystal ligand of PLpro, i.e., Vir251 showed CDocker value − 75.038 kcal mol−1 with its target enzyme. Moreover, the numbers of H-bond formation between the ligand and receptor were analyzed and observed that both SARS-CoV-2 Mpro-N3 and SARS-CoV-2 Mpro-ECG complexes were formed with six H-bonds. It is assumed that more numbers of H-bonds give a better tolerance to the mutability of the virus [2]. Because ECG showed the highest binding affinity with Mpro and formed more numbers of H-bonds in comparison to all the test compounds, the complex of SARS-CoV-2 Mpro- ECG was taken for further studies, whereas, the complex SARS-CoV-2 Mpro-N3 was taken as control. Therefore, these two complexes were further analyzed by MM-PBSA binding energy calculation during 100 ns MD simulations to compute the binding behavior of the ligand to the receptor by mimicking in vitro and in vivo conditions [31, 50].

From the results obtained from the MM-PBSA analysis, the average binding free energies (ΔGBind) of SARS-CoV-2 Mpro–N3 and SARS-CoV-2 Mpro–ECG complexes were calculated during the 100 ns simulation. The resulting ΔGBind of SARS-CoV-2 Mpro–N3 and SARS-CoV-2 Mpro–ECG complexes were found to be − 282.96 ± 32.97 kJ mol–1 and − 192.40 ± 27.10 kJ mol–1, respectively. The details of MM-PBSA calculation of the complexes are summarized in Table 2. The ΔGBind plot of the complexes is depicted in Fig. 1a.

Table 2

Average statistics of MM-PBSA calculation

MMPBSA
Protein–ligand complexes	ΔEMM (kJ mol–1)		ΔGSol (kJ mol–1)		ΔGBind (kJ mol–1)
	ΔEVdW (kJ mol–1)	ΔEElec (kJ mol–1)	ΔGPB (kJ mol–1)	ΔGSASA (kJ mol–1)
N3–SARS-CoV-2 Mpro	− 361.67 ± 33.14	− 50.54 ± 21.37	158.36 ± 26.85	− 29.11 ± 1.95	− 282.96 ± 32.97
ECG–SARS-CoV-2 Mpro	− 275.01 ± 18.68	− 76.97 ± 29.45	179.86 ± 30.32	− 20.28 ± 1.11	− 192.40 ± 27.10

ΔE Van der Waals contribution from MM, ΔE electrostatic energy as calculated by the MM force field, ΔG solvation free energy comprising ΔGPB, the polar (ΔGpolar) and non-polar (ΔGnon-polar) solvation free energy contribution calculated by PB equation, ΔG the energy contribution from solvent-accessible surface area (SASA), ΔG binding free energy

Fig. 1

a Free energy of binding (ΔGBind), b molecular mechanics potential energy (ΔEMM), c polar solvation energy (ΔGPolar) under ΔGPB and d Non-polar solvation energy (ΔGNon-polar or ΔGSASA) under ΔGPB between the protein and ligand during the simulation was calculated using g_mmpbsa tool. Calculations were performed on the 100 ns trajectories at the interval of every 1 ns. In all the figures, the black line represents SARS-CoV-2 Mpro -N3 complex and red lone represents SARS-CoV-2 Mpro -ECG

The calculated averages of molecular mechanics potential energy (ΔEMM), polar solvation energy (ΔGPolar) and non-polar solvation energy (ΔGNon-polar or ΔGSASA) under ΔGPB between the protein and ligand during the simulation were evaluated from MMPBSA calculation. Average ΔEMM for SARS-CoV-2 Mpro–N3 and SARS-CoV-2 Mpro–ECG complexes were found to be − 412.21 ± 41.69 and –351.98 ± 32.10 kJ mol–1, respectively, and the average ΔGPolar for SARS-CoV-2 Mpro–N3 complex and SARS-CoV-2 Mpro–ECG complex were found to be 158.36 ± 26.85 and 179.86 ± 30.32 kJ mol–1, respectively. Calculated average of ΔGNon-polar or ΔGSASA for SARS-CoV-2 Mpro–N3 complex and SARS-CoV-2 Mpro–ECG complex (red) were found to be − 29.11 ± 1.95 and − 20.28 ± 1.11 kJ mol–1, respectively. The graphical representations of ΔEMM, ΔGPolar and ΔGNon-polar or ΔGSASA as function of time during the simulation are shown in Fig. 1.

To compare the binding of ECG and N3 to SARS-CoV-2 Mpro, three MD simulations were performed on three different states (SARS-CoV-2 Mpro, SARS-CoV-2 Mpro–N3 complex, and SARS-CoV-2 Mpro–ECG complex) for 100 ns. In the MD simulations, lower RMSD values give the idea of the conformational stability of the protein–ligand complex [51]. It was observed that the root-mean-square deviations (RMSDs) of SARS-CoV-2 Mpro system (Fig. 2a) were stabilized around 0.45 nm from 60 ns onward. Whereas, the RMSDs were stabilized from 30 ns onward for SARS-CoV-2 Mpro–N3 at 0.4 nm and SARS-CoV-2 Mpro–ECG systems at 0.35 nm (Fig. 2a). There were less conformational changes in SARS-CoV-2 Mpro–N3 and SARS-CoV-2 Mpro–ECG systems during the simulations in comparison to the SARS-CoV-2 Mpro system. However, RMSD of SARS-CoV-2 Mpro–ECG system was found to be almost similar to RMSD of SARS-CoV-2 Mpro–N3 system.

Fig. 2

Analytical depiction of root-mean-square deviation (RMSD) of three systems. a RMSDs of the protein with least-square fit to protein in SARS-CoV-2 Mpro (black), SARS-CoV-2 Mpro–N3 (red) and SARS-CoV-2 Mpro–ECG (green) systems. b RMSDs of the heavy atoms of small molecules N3 (black) and ECG (red) with least-square fit to the protein backbone. c Radius of gyration (Rg) of SARS-CoV-2 Mpro (black), SARS-CoV-2 Mpro–N3 (red) and SARS-CoV-2 Mpro–C31 (green) systems. The plot clearly depicts Rg of the three systems by their nature. Ligand-bound systems are seemed to reach their stably folded states earlier than the ligand-free system. d Self-explanatory depiction of root-mean-square fluctuations of SARS-CoV-2 Mpro in SARS-CoV-2 Mpro–only (black), SARS-CoV-2 Mpro–N3 (red) and SARS-CoV-2 Mpro–ECG (green) systems

The RMSDs of the heavy atoms of small molecules N3 and ECG to the backbone of SARS-CoV-2 Mpro indicated that both the ligands were unstable with respect to the protein (Fig. 2b). However, a noticeable relative conformational change was observed in the case of ECG with ~ 0.7 nm higher between 56 and 60 ns than that of the relative conformational change of N3. Additionally, RMSDs of the ligands were calculated by least-square fitting to their simulations and observed that both ligands also were found to fit similarly to the binding pocket (Supplementary Fig. 3).

To study the conformational stability and protein compactness, radius of gyration (Rg) of the SARS-CoV-2 Mpro protein was calculated in all three states. It was observed that Rg of SARS-CoV-2 Mpro without the ligands (Fig. 2c, black) significantly decreased from 2.2 nm to ~ 2.1 nm and from 55 ns onward Rg stabilized throughout the simulation. Whereas in case of SARS-CoV-2 Mpro in the presence of N3 and ECG ligands (Fig. 2c, red and green, respectively), Rgs get stabilized near 2.15 nm from 30 ns onward. Although, Rg reduced drastically in the initial phase, but stable nature of Rg from 30 ns onward suggests the stable folding of the protein along with the ligands.

The residue wise root-mean-square fluctuation (RMSF) of SARS-CoV-2 Mpro was calculated for SARS-CoV-2 Mpro–only, SARS-CoV-2–N3 and SARS-CoV-2 Mpro–ECG systems during the course of MD simulation. From Fig. 2d for SARS-CoV-2 Mpro–only plot (black), it was observed that there were distinctly variable peaks in residue positions 21–24, 30–34, 44–56, 92–110, 117–121, 130–145, 182–197 and 220–305 indicating higher fluctuations in these residue locations than that of the SARS-CoV-2–N3 and SARS-CoV-2 Mpro–ECG systems. Although the RMSF of SARS-CoV-2–N3 and SARS-CoV-2 Mpro–ECG systems revealed almost similar pattern of residue fluctuation, residue positions of 44–50 in SARS-CoV-2–N3 system and residue positions of 271–286 in SARS-CoV-2 Mpro–ECG system were noticeably flexible in the plot (Fig. 2d).

To evaluate the interactions of the ligands with SARS-CoV-2 Mpro during the simulation, another two sets of MD simulations of 100 ns were performed on each of the ligands N3 and ECG with solvent water molecules. It further revealed that the small molecule ECG only exhibits a linear LJ interaction at around − 250 to − 280 kJ mol−1 (Fig. 3a). Short-range coulombic interactions seemed to be at the same point for both N3 and ECG (Fig. 3b). From solvent accessible surface area calculation, ΔG of solvation of SARS-CoV-2 Mpro in SARS-CoV-2 Mpro–N3 and SARS-CoV-2 Mpro–ECG systems were deduced and were represented in Supplementary Fig. 4.

Fig. 3

Protein–ligand interaction analysis. a Short-range coulombic (black) and Lennard–Jones potentials (red) of inhibitor N3 interacting with SARS-CoV-2 Mpro. b Short-range coulombic (black) and Lennard–Jones potentials (red) of inhibitor ECG interacting with SARS-CoV-2 Mpro

Detailed interactions of inhibitors N3 and ECG with SARS-CoV-2 Mpro were studied through molecular visualization using molecular graphics tools like, UCSF Chimera and LigPlot + . The X-ray crystal structure of SARS-CoV- 2 Mpro obtained from PDB (PDB ID: 6LU7) revealed the original H-bonding interactions of the inhibitor N3 with the residues Gly143, HIS163, HIS164, GLN166, GLN189 and THR190 (Fig. 4a). Interestingly, after 100 ns of simulation, it was observed that only three residues of SARS-CoV-2 Mpro namely GLU166, ALA191, and GLN192 were involved in H-bonding interactions with inhibitor N3 (Fig. 4b). Nevertheless, the inhibitor N3 remained in the binding pocket of SARS-CoV-2 Mpro due to the strong Van der Waals (VdW) interactions between the residues, GLN189, MET165, HIS163 164, ALA193 194, THR190, ASN142, LEU167, PHE185, HIS41, MET149, LEU50, GLY143, and CYS44. Molecular surface (MS) representation revealed that inhibitor N3 anyhow occupied the binding pocket of SARS-CoV-2 Mpro from 0 to 100 ns (Fig. 4c, d).

Fig. 4

Interaction of inhibitor N3 (DRG) with the key residues at the binding pocket of SARS-CoV-2 Mpro. a Ligplot depiction of hydrogen bonding and hydrophobic interactions between N3 and SARS-CoV-2 Mpro as found in PDB entry 6LU7 at its native state. b Ligplot depiction of hydrogen bonding and hydrophobic interactions between N3 and SARS-CoV-2 Mpro as at 100 ns. c Inhibitor N3 at the binding pocket of SARS-CoV-2 Mpro as found in PDB entry 6LU7 at its native state (at 0 ns). d Inhibitor N3 at the binding pocket of SARS-CoV-2 Mpro as at 100 ns. e Calculation of distances between the atoms of N3 and SARS-CoV-2 Mpro residues forming h-bonds depicted from a during the simulation. f Calculation of distances between the atoms of N3 and SARS-CoV-2 Mpro residues forming h-bonds depicted from b during the simulation

Distances between the atoms of the amino acid residues of SARS-CoV-2 Mpro and atoms of the small molecule were measured using the gmx distance program of GROMACS to evaluate the H bonding interaction. Calculated distances of the atomic pairs involved in H-bonding are depicted in Fig. 4a, b. Figure 4e, f represents the evolution of H-bonds throughout the simulation with respect to the configuration at 0 ns and 100 ns.

A similar analysis is presented for SARS-CoV-2 Mpro–ECG complex in Fig. 4. Initially, SARS-CoV-2 Mpro–ECG docked complex possessed six H-bonds, two from PHE140 and SER144 each, one from HIS164 and the other from MET165 (Fig. 5a). However, the residues aforementioned for the inhibitor N3 (shown in Fig. 4) were also contributing strong VdW interactions with ECG. Interestingly at the end of the 100 ns simulation, the SARS-CoV-2 Mpro- ECG complex possessed seven H-bonds. The residues involved in H-bonding with ECG were HIS41, LEU50, GLY143, LEU141, GLU166, ASP187 and GLN189 (Fig. 5b).

Fig. 5

Interaction of inhibitor ECG (LIG) with the key residues at the binding pocket of SARS-CoV-2 Mpro. a Ligplot depiction of hydrogen bonding and hydrophobic interactions between ECG and SARS-CoV-2 Mpro at its native state from docking experiment. b Ligplot depiction of hydrogen bonding and hydrophobic interactions between ECG and SARS-CoV-2 Mpro at 100 ns. c Inhibitor ECG at the binding pocket of SARS-CoV-2 Mpro at its native state from docking experiment (at 0 ns). d Inhibitor ECG at the binding pocket of SARS-CoV-2 Mpro at 100 ns. e Calculation of distances between the atoms of ECG and SARS-CoV-2 Mpro residues forming h-bonds depicted from a during the simulation. f Calculation of distances between the atoms of ECG and SARS-CoV-2 Mpro residues forming H-bonds depicted from b during the simulation

On comparing N3-Mpro and ECG-Mpro complexes, it was observed that the inhibitor ECG showed better binding when compared to N3 in terms of number of H-bonds and the fitness of the molecules at the binding pocket. Fitness of ECG at the binding pocket of SARS-CoV-2 Mpro could be visualized in MS representation in Fig. 5c, d and it was observed that the ECG was going deep inside the pocket during the simulation at 100 ns (Fig. 5d). Further, ECG was found to attach to the binding pocket of SARS-CoV-2 Mpro more firmly than N3 with more and consistent numbers of H-bonds with potential residues at SARS-CoV-2 Mpro pocket (Fig. 5e, f).

This analysis indicated that our predicted compound might be more permissive to future mutations. Due to the good binding affinity, (−)-epicatechin-3-O-gallate, an abundant green tea flavonoid [52] might act as a potential inhibitor of SARS-CoV-2 Mpro. pkCSM results also revealed good pharmacokinetic properties of the suggested compound (Supplementary Table 2). This analysis indicated that our predicted compound might be more permissive to future mutations. Due to the good binding affinity, (−)-epicatechin-3-O-gallate, an abundant green tea flavonoid might act as a potential inhibitor of SARS-CoV-2 Mpro. To analyze the drug-like properties of ECG, the pharmacokinetic parameters were studied. For the effectiveness of a drug, satisfying all these parameters are very necessary. In this study, we observed good Caco2 permeability of ECG with 63% of intestinal absorption. According to the algorithm used in pkCSM, the drugs having more than 30% intestinal absorption are considered as ideal. In analyzing the various parameters of drug metabolism, it was observed that the compound neither acts as substrate nor as an inhibitor of the different isoforms of Cytochrome P450 except CYP1A2. ECG shows the inhibitory effect on CYP1A2. To analyze the deposition and clearance of the compound, renal OCT2 substrate activity of the compound was analyzed. ECG did not show the substrate activity for renal OCT2. Apart from that, the toxicity assessment of ECG was done by analyzing AMES toxicity, hepatotoxicity and hERGI/II inhibitory activity. From the results, it was observed that ECG did not show any carcinogenic activity in ADME-Tox study without having a hepatotoxic effect. It has hERGII inhibition activity while it did not show any inhibitory activity for hERGI. The results of pharmacokinetic parameters are given in Supplementary Table 2. Depending on the various parameters mentioned above, the compound can be considered easily for developing into an independent drug molecule.

The HOMO and LUMO frontier orbital energies of N3 and ECG were calculated by DFT analysis. The calculated HOMO and LUMO energies of N3 and ECG obtained from the DFT are given in Table 3. These physical parameters provide the information on the chemical and biological activities of the compounds. The decreased HOMO–LUMO gap or energy gap indicates the higher reactivity of the compound with lesser stability [53]. From the calculated HOMO and LUMO energies, the energy gap of the compounds was calculated and observed that the energy gap of ECG (0.107783 eV) is slightly higher than the co-crystal ligand N3 (0.091682 eV) indicating similar reactive nature of ECG with N3.

Table 3

Comparison of the values of orbital energy descriptors HOMO and LUMO

Molecule	HOMO energy (au)	LUMO energy (au)	Energy gap (eV)
N3	− 0.18692	− 0.095238	0.091682
ECG	− 0.19072	− 0.0829372	0.107783

Since both N3 and ECG showed almost similar reactivity, therefore we further determined the IC50 value of ECG and N3 by performing QSAR analysis. During QSAR analysis, out of 21 inhibitors, the training set was randomly generated with 80% molecules (17 inhibitors) and rest were kept as a test set for validation of the model. A total of 100 QSAR models were generated using genetic function approximation model and the best model having R2 value of 88.02% was selected (adjusted R2 = 84.02%) (Supplementary Table 3).

The generated QSAR equation was: “Log (IC50)−1 =  − 7.0395 + 0.014618 * Molecular_Weight − 0.2866 * Num_Rotatable Bonds − 0.0051847 * CHARMm Energy + 0.00042285 *Total_Energy_DMol3”.

Further, the activity of ECG and N3 against SARS-CoV-2 Mpro was predicted using the QSAR equation by generating the plot presented in Fig. 6.

Fig. 6

QSAR Plot showing Experimental Activity [log(IC50)−1] against Predicted Activity [log(IC50)−1]

For analyzing the activity of ECG and N3, the molecules were minimized using the same CHARMm module and requisite descriptors were calculated in Discovery Studio 2019. The calculated descriptors for ECG were: Molecular_Weight = 442.372, Numbers of _Rotatable Bonds = 4, CHARMm Energy =  − 62.6023, Total_Energy_DMol3 =  − 1588.42 and for N3: Molecular_Weight = 680.791, Numbers of _Rotatable Bonds = 18, CHARMm Energy =  − 117.884, Total_Energy_DMol3 =  − 2273.32. After putting these variables in the QSAR Equation as mentioned above, predicted IC50 value for ECG was found to be 0.3281 nM whereas predicted IC50 value for N3 was 394.982 nM. These results suggest the presence of more inhibitory potential in ECG for SARS-CoV-2 Mpro in comparison to N3.

Conclusion

Docking results revealed that ECG exhibited the lowest CDocker energy with value − 44.72 kcal mol−1 or − 187.108 kJ mol−1 which was further validated by MMPBSA analysis and could be correlated to the average of ΔGBind that was found to be − 192.40 ± 27.10 kJ mol−1 during 100 ns MD simulation. However, the average of ΔGBind of N3 calculated by MMPBSA method during MD simulation was found to be underestimated as − 282.96 ± 32.97 kJ mol−1 as CDocker energy during docking was found to be − 91.37 k kcal mol−1 or − 382.29 kJ mol−1. Nevertheless, N3 possessed stronger binding affinity toward SARS-CoV-2 Mpro, but taking into account of ECG at the binding pocket of the enzyme, the RMSDs profile of the heavy atoms of N3 and ECG to the backbone of SARS-CoV-2 Mpro indicated that both the ligands were unstable with respect to the protein (Fig. 2b). Moreover, the RMSF analysis revealed the most flexible residue positions as 44–50 and 271–286 in SARS-CoV-2–N3 and SARS-CoV-2 Mpro- ECG systems, respectively. From molecular visualization of SARS-CoV-2–N3 and SARS-CoV-2 Mpro–ECG systems, it was revealed that the inhibitor ECG conferred better binding in comparison to N3 in terms of number of H-bonds and the fitness of the molecules at the binding pocket. ECG was found be located inside the binding pocket of SARS-CoV-2 Mpro more firmly than N3 with more and consistent numbers of H-bonds with potential residues at SARS-CoV-2 Mpro pocket (Fig. 5e, f).

Based on the series of computational analysis including molecular docking, MD simulation analysis and subsequent QSAR analysis, we propose ECG as a potential inhibitor of SARS-CoV-2 Mpro with predicted IC50 value of 0.3281 nM. In this current situation, a synergistic approach seems to be the most plausible choice for mitigating the spread of the pandemic. Although this study is solely based on in silico predictions and the results need further in vitro and in vivo validations, the present findings may positively contribute to the current global search for inhibitors of SARS-CoV-2 Mpro.

Below is the link to the electronic supplementary material.

Supplementary file1 (PDF 209 kb)