2. Materials and methods 2.1. Molecular docking studies 2.1.1. Protein preparations The crystal structure of SARS-CoV-2 RdRp (PDB ID: 6M71) (Yan et al., 2020) was retrieved from the protein databank (www.rcsb.org) (Berman et al., 2000). The crystal structure was prepared individually by adding hydrogen atoms and computing the Gasteiger charge using the AutoDock v4.2 program (Morris et al., 2009). Subsequently, the file was saved as .pdbqt format in preparation for molecular docking. Schematic representation of the work-flow for selecting potential natural polyphenolic inhibitors for the SARS-CoV-2 RdRp is shown in Figure 2. Figure 2. Flow chart of the methodology for shortlisting the best natural polyphenolic inhibitor of the SARS-CoV-2 RdRp. 2.1.2. Ligand preparations The SDF structures of GTP, remdesivir, and selected hundred polyphenols (see Table S1 in Supplementary Information) were retrieved from the PubChem database (https://pubchem.ncbi.nlm.nih.gov/) (Kim et al., 2019). The compounds were converted into PDB format, and conformational energies of all the compounds were minimized by using UCSF Chimera (Pettersen et al., 2004). Table 1. Binding energy (kcal/mol) of the natural polyphenols along with the control compounds (GTP and remdesivir) against RdRp of the SARS-CoV-2 (PDB ID: 6M71) by molecular docking study. S. No. Compound Name Binding energy (kcal/mol) S. No. Compound Name Binding energy (kcal/mol) 1 TF3 −9.9 52 Cyanidin −6.3 2 TF2b −9.6 53 Daidzein −6.3 3 TF1 −9.6 54 Glycitein −6.3 4 TF2a −9.3 55 Wogonin −6.3 5 Hesperidin −8.8 56 Phloretin −6.3 6 EGCG −7.3 57 Catechin −6.2 7 Myricetin −7.2 58 Urolithin B −6.2 8 Quercetagetin −7.0 59 Angolensin −6.2 9 Quercetin −6.9 60 Pinosylvin −6.2 10 Curcumin −6.9 61 Formononetin −6.2 11 Dihydrorobinetin −6.8 62 Liquiritigenin −6.2 12 Peonidin −6.8 63 Prunetin −6.2 13 Fisetin −6.8 64 Alpinetin −6.2 14 Robinetin −6.7 65 Biochanin A −6.2 15 5-Deoxygalangin −6.7 66 Rhapontigenin −6.1 16 Kaempferol −6.7 67 Genistein −6.1 17 Scutellarein −6.7 68 Chrysin −6.1 18 (-)-Epicatechin −6.7 69 6-Hydroxyflavone −6.1 19 Purpurin −6.7 70 Equol −6.1 20 Isorhamnetin −6.7 71 Piceatannol −6.1 21 Tricetin −6.6 72 Isorhapontigenin −6.0 22 Gossypetin −6.6 73 Resveratrol −5.8 23 Norathyriol −6.6 74 Danshensu −5.7 24 Coumestrol −6.6 75 Eugenin −5.6 25 Isosakuranetin −6.6 76 Sinapic acid −5.5 26 Pectolinarigenin −6.6 77 Pterostilbene −5.5 27 Tangeritin −6.6 78 Ferulic acid −5.4 28 Nobiletin −6.6 79 Caffeic acid −5.4 29 Pratensein −6.6 80 Isoferulic acid −5.4 30 Hispidulin −6.6 81 Dihydrocaffeic acid −5.4 31 Baicalein −6.5 82 Gentisic acid −5.3 32 Apigenin −6.5 83 Pyrogallol −5.3 33 Morin −6.5 84 4-Hydroxycinnamic acid −5.2 34 Urolithin A −6.5 85 Resacetophenone −5.2 35 Acacetin −6.5 86 Salicyclic acid −5.1 36 Pelargonidin −6.5 87 Syringic acid −5.1 37 Irilone −6.5 88 2-Hydroxybenzoic acid −5.1 38 Naringenin −6.5 89 Gallic acid −5.0 39 Pinocembrin −6.5 90 3-Hydroxybenzoic acid −5.0 40 Kaempferide −6.5 91 4-Hydroxybenzoic acid −5.0 41 Malvidin −6.5 92 Vanillin −5.0 42 Luteolin −6.4 93 p-Coumeric acid −4.9 43 Dalbergin −6.4 94 Vanillic acid −4.8 44 Butein −6.4 95 Paeonol −4.8 45 Biochanin A (1-) −6.4 96 Cinnamic acid −4.7 46 Fustin −6.4 97 Protocatechuic acid −4.6 47 5-Hydroxyflavone −6.4 98 4-Ethylphenol −4.5 48 Pinostrobin −6.4 99 Catechol −4.5 49 Pinobanksin −6.4 100 Tyrosol −4.5 50 Datiscetin −6.3 101 GTP −7.9 51 Galangin −6.3 102 Remdesivir −7.7 2.1.3. Docking studies using AutoDock Vina The energy-minimized structure of all the natural polyphenols, remdesivir, and GTP were docked with the receptor (RdRp of SARS-CoV-2) using AutoDock Vina 1.1.2 (Trott & Olson, 2010). The ligand files were further saved in PDBQT file format, a modified PDB format containing atomic charges, atom type definitions for ligands, and topological information (rotatable bonds). A grid box (30 Å × 30 Å × 30 Å) centered at (121, 120, 125) Å for the SARS-CoV-2 RdRp, was used in the docking experiments. After the receptor-ligand preparation, docking runs were started from the command prompt. The lowest binding energy and best-docked conformation were considered as the ligand molecule with maximum binding affinity. 2.1.4. Protein-ligand interactions LigPlot+ was used to investigate protein-ligand interactions for a given .pdb file containing the docked conformation and also the final simulated conformation (Wallace et al., 1995). The LigPlot+ program self-generated schematic 2D representations of protein-ligand interaction. The output file represents the intermolecular interactions and their strengths, including hydrogen bonds, hydrophobic contacts, and atom accessibilities. H-bonds are shown in green dotted lines, whereas residues involved in hydrophobic interaction are represented in the red semicircle. 2.2. Molecular dynamics simulations All-atom molecular dynamics (MD) simulations were performed on the best eight selected plant-derived natural polyphenols obtained from the molecular docking study along with remdesivir, a well-known RdRp inhibitor, for studying thermodynamic stability of the docked structure. The pmemd.cuda module in AMBER18 (Case et al. 2018) was used for conducting MD simulations, and all simulations were performed utilizing the graphics card, RTX 2080Ti. We adopted the same protocol that was used in our earlier studies (Jonniya et al., 2019; Sk, Roy, Jonniya, et al., 2020; Sk, Roy, & Kar, 2020). The receptor and small molecules were described by the Amber ff14SB (Maier et al., 2015) and GAFF2 (Wang et al., 2004) force field, respectively. Ligands were assigned the AM1-BCC (Jakalian et al., 2002) atomic charges calculated using the antechamber (Wang et al., 2001) module. The complexes were then solvated using an explicit TIP3P (Price & Brooks III, 2004) water model, and nearly 38124 water molecules were needed to solvate each system. Subsequently, all solvated systems were neutralized by adding an appropriate number of Na+ ion. All bond lengths, including hydrogen atoms, were constrained by the SHAKE algorithm (Kräutler et al., 2001). This allows the usage of a 2 fs time-step. The non-bonded cut-off was set to 8 Å and the long range electrostatic interactions were evaluated using the particle-mesh Ewald (PME) (Darden et al., 1993) method. The temperature was kept at 300 K using the Langevin thermostat (Loncharich et al., 1992) with a collision frequency of 2 ps−1. The system pressure was controlled by Berendsen’s Barostat (Berendsen et al., 1984) and fixed at 1.0 bar. We used a time-step of 2.0 fs for all simulations. Briefly, we used two stages of minimization. Firstly, each complex was optimized by 500 steps of steepest descent followed by another 500 steps of conjugate gradient minimization, keeping all atoms of the complex restrained to their initial coordinate with a weak harmonic potential (force constant 2.0 kcal mol−1 Å−2). The second stage of minimization was carried out without any restraints by performing 100 steps of steepest descent, followed by another 900 steps of conjugate gradient minimization to remove any residual steric clashes. Subsequently, all systems were gradually heated from 0 K to 300 K at the NVT ensemble with a force constant of 2.0 kcal mol−1Å−2 acting on all solute atoms. Next, 1.0 ns equilibration MD phase was carried out without any restraint. Finally, we performed 150 ns production simulations for all four systems at the NPT ensemble. Overall, we accumulated 15000 conformations for each simulation, and we used 500 snapshots from the last 50 ns trajectories for binding affinity estimation using the MM-PBSA scheme. The trajectory analysis was done by the AmberTools19 CPPTRAJ (Roe & Cheatham III, 2013) module of Amber18. 2.3. Molecular mechanics Poisson-Boltzmann surface area calculations The binding affinity of EGCG, TF1, TF3, TF2b, TF2a, hesperidin, myricetin, quercetagetin, and remdesivir against the SARS-CoV-2 RdRp, were calculated by the molecular mechanics Poisson-Boltzmann surface area (MM-PBSA) methodology (Jonniya & Kar, 2020; Kar, Seel, et al., 2007; Kar, Wei, et al., 2007; Kar et al., 2011; 2013; Kar & Knecht, 2012a, 2012b, 2012c, 2012d; Kollman et al., 2000; Roy et al., 2020; Wang et al., 2006). The MM-PBSA scheme can briefly be described as follows: (1) ΔGbind=ΔH−TΔS≈ΔEMM+ΔGsolv−TΔS (2) ΔEMM=ΔEinternal+ΔEelec+ΔEvdW (3) ΔGsolv=ΔGpol+ΔGnp where ΔEMM, ΔGsolv, TΔS are the changes in molecular mechanical energy, solvation free energy, and conformational entropy, respectively. Further, molecular mechanical energy is composed of ΔEinternal (bond, dihedral, and angle), ΔEelec (electrostatic) and ΔEvdW (van der Waals) and the change in desolvation free energy is composed of polar solvation (ΔGpol) and non-polar solvation free energy (ΔGnp). The polar solvation free energy, ΔGpol, was calculated by the pbsa module of AMBER18. Due to the high computational cost, we neglect the configurational entropy calculations. Further, to understand the polyphenol-protein interaction more closely, the interaction energy was decomposed into the contributions from each residue of the protein by using the molecular mechanics generalized Born surface area (MM-GBSA) scheme (Gohlke et al., 2003). 2.4. ADMET studies The in-silico pharmacological studies of EGCG, TF2a, TF2b, TF3, and remdesivir were predicted based on their ADMET profile. The ADMET studies (absorption, distribution, metabolism, elimination, and toxicity) were predicted using the pkCSM tool (http://biosig.unimelb.edu.au/pkcsm/prediction) (Pires et al., 2015). The canonical SMILE molecular structures of the above-mentioned compounds were retrieved from the PubChem database (www.pubchem.ncbi.nlm.nih.gov). 2.5. Molecular target prediction Natural compounds interact with a large number of proteins, enzymes, lipids. This interaction plays a crucial role in elucidating the molecular mechanism of the small molecules. So, it is important to identify the molecular targets for new molecules (Gfeller et al., 2014). Swiss Target Prediction website (http://www.swisstargetprediction.ch/index.php) was logged on, and canonical SMILE molecular structures of remdesivir, EGCG, TF2a, TF2b, and TF3 were entered in the search bar option, and results were analyzed.